text:
# Constraints on Galaxy Density Profiles from Strong Gravitational Lensing: The Case of B 1933+503

## 1 Introduction

Strong gravitational lenses can aid in the determination of two astrophysical quantities of great interest: the density profile of massive galaxies and the value of the Hubble constant.

The density profile of massive galaxies is not known. Massive early-type and late-type galaxies appear to have mass distributions consistent with flat rotation curves (total mass $`Mr`$, density $`\rho 1/r^2`$) on scales outside their cores and inside their outer edges. This is most cleanly shown by spiral galaxy rotation curves (e.g. Rubin, Ford & Thonnard rft78 (1978), rft80 (1980)), elliptical galaxy X-ray halos (Fabbiano fabbiano (1989)), and some gravitational lenses (e.g., Kochanek csk1654 (1995)). It is also supported by modern stellar dynamical studies of early-type galaxies (e.g., Rix et al. rixetal (1997)). The central properties of galaxies are characterized by density cusps, rather than finite core radii. In particular, the central luminosity profiles of massive galaxies generically show density cusps with $`\rho r^\gamma `$ and $`1\gamma 2`$ (e.g. Faber et al. faber (1997)). Halos in dark matter simulations, which appear to form with $`\gamma 11.5`$ cusps (e.g. Navarro, Frenk & White nfw (1997), Klypin et al. Kly (1999), Moore et al. Moo98 (1998) ), suggest that the dark matter cores should have similar properties. However, it is expected that the dark matter profile will be significantly modified by the subsequent evolution and interaction with the baryons (e.g. Blumenthal et al. blum86 (1986), Dubinski dubin (1994), Navarro, Eke & Frenk nef96 (1996), Gelato & Sommer-Larsen gelato99 (1999), Binney, Gerhard & Silk binney00 (2000)). Currently popular scenarios for self-interacting dark matter (SIDM, see e.g. Spergel & Steinhardt SIDM1 (2000); Kochanek & White SIDM2 (2000)) further complicate predictions for the central density structure of galaxies. The resulting central density profile of the dark matter in massive galaxies is not known. Even less is known about the angular shapes of galaxy mass distributions, although spiral galaxies are clearly far rounder than their luminous disks (see Sackett sackett (1998) for a review) and ellipticals can be either flatter or rounder than their stellar profile (see Buote & Canizares buca (1997) and references therein) and are modestly triaxial (e.g. Franx, Illingworth & de Zeeuw franx (1991)). See, e.g., Rusin and Tegmark (rusin (2000)) for a recent summary and more references.

Multiple-image gravitational lenses offer a promising, newer probe of galaxy mass distributions. The image geometry usually constrains the mass inside the critical line of the lens to an accuracy exceeding that of any other method111The accuracy is almost always better than $`5\%`$ including systematics such as the cosmological model. As is true for all mass estimation methods, there is a global scaling with the Hubble constant (distance), which for the lenses is simply that the mass is proportional to $`h^1`$.. Unlike other probes of galaxy mass distributions, lenses also measure the shape of the gravitational potential near the critical line accurately. In most systems, however, this well measured quantity includes contributions from the lens galaxy and tidal perturbations (either local to the lens or along the line of sight) which can be difficult to disentangle (see e.g. Bar-Kana barkana (1996), Keeton, Kochanek & Seljak kks (1997)). In the majority of lenses, our ability to measure the radial or angular mass distribution is limited by the paucity of constraints supplied by the small numbers of images (2 or 4). There are, however, many lenses with either a larger number of discrete images or continuous arcs or rings formed from the host galaxy or extended radio emission regions. These are the systems best used to probe the details of galaxy mass distributions using gravitational lensing.

Gravitational lenses are also a promising method for determining the global value of the Hubble constant by measuring the light propagation delays between the images (Refsdal refsdal64 (1964)). Although time delays have been measured for six gravitational lenses, the inferred values for the Hubble constant are limited by the systematic uncertainties in the mass models for the systems (see, e.g. Impey et al. impey (1998), Koopmans & Fassnacht koop (1999), Williams & Saha williams (2000), Keeton et al. keeton2k (2000) for recent examples). The uncertainties can be reduced by finding additional constraints on the systems with delays (see Kochanek, Keeton & McLeod kkm (2000)), including external constraints on the models, or by focusing future time delay studies on lenses with fewer modeling uncertainties.

Thus, the utility of strong gravitational lensing is maximized in systems with large numbers of constraints. The focus of this paper is one such system, the 10 image lens B 1933+503 (Sykes et al. sykes (1998), Marlow et al. marlow (1999), Chapman et al. chapman (1998), Browne et al. brow (1999), Biggs et al. biggs (2000), Norbury et al. norbury (2000)). B 1933+503 consists of two quadruply imaged radio sources and one doubly imaged radio source produced by a $`z_l=0.76`$ lens galaxy. The source redshift is $`z_s=2.62`$ (Biggs et al biggs (2000), Norbury et al. norbury (2000)). The images range from unresolved to a pair of thin arcs (see the schematic in Figure 1).

Nair (nair (1998)) showed that a simple singular isothermal ellipsoid (SIE) could fit the available constraints (see also Chapman et al. chapman (1998) for some discussion of modeling). In this paper we extend her analysis to a wide range of parametric mass distributions to explore the ways in which the complicated image pattern restricts the mass distribution. From this ensemble of models, we identify common features of the acceptable models, such as the shape of the mass profile/rotation curve, and estimate the uncertainties in the image time delays for use in determining the Hubble constant. In §2 we review the necessary mathematics of lensing and describe the available data. In §3 we catalog the models we use and their results. We summarize the results and suggest further directions in §4.

## 2 Models

In this section we review the mathematical background and notation for our models (§2.1), discuss the mass distributions we use (§2.2) and the available data on B 1933+503 (§2.3).

### 2.1 Notation

We consider a thin lens characterized by a projected lens potential $`\varphi (𝐱)`$ which is a function of a Cartesian angular coordinate $`𝐱`$. The potential is related to the surface mass density by $`^2\varphi (𝐱)=2\kappa (𝐱)`$ where $`\kappa (𝐱)`$ is the surface density in units of the critical surface density. In addition to the primary lens, we include an external tidal field defined by an amplitude $`\gamma `$ and a orientation $`\theta _\gamma `$. For a source at angular position $`𝐮`$, the images are found at solutions $`𝐱_i`$ of the lens equation

$$𝐮=𝐱\varphi (𝐱)\mathrm{\Gamma }𝐱$$
(1)
where the shear tensor is

$$\mathrm{\Gamma }=\gamma \left(\begin{array}{cc}\hfill \mathrm{cos}2\theta _\gamma & \hfill \mathrm{sin}2\theta _\gamma \\ \hfill \mathrm{sin}2\theta _\gamma & \hfill \mathrm{cos}2\theta _\gamma \end{array}\right).$$
(2)
The shear term models not only local tidal perturbations to the lens but many of the additional perturbations along the line of sight (see Kovner kovner (1987), Bar-Kana barkana (1996)), and it is a necessary component of any realistic lens model (see Keeton, Kochanek & Seljak kks (1997)). The inverse magnification tensor of the model is

$$M^1=I\mathrm{\Gamma }\left(\begin{array}{cc}\hfill \varphi ,_{xx}& \hfill \varphi ,_{xy}\\ \hfill \varphi ,_{xy}& \hfill \varphi ,_{yy}\end{array}\right).$$
(3)
and the magnification of an unresolved image is given by det $`M`$. Equation (1) can also be written as the condition for an extremum on a time delay surface, and the surface itself determines the time delay (see Schneider, Ehlers & Falco lensbook (1992)).

The procedures for modeling the positions of lensed point sources are simple to describe and relatively easy to implement. Equation (1) can be read as a map from the image plane to the source plane via a lens model (potential $`\varphi `$ and external shear $`\mathrm{\Gamma }`$). The calculations in this paper used the lensmodel code (Keeton 2000a ). For observed image positions $`𝐱_i`$ and model image positions $`𝐱_i^{mod}(𝐮)`$ the goodness of fit is measured with a simple $`\chi ^2`$ statistic

$$\chi _{point}^2=\underset{i}{}\frac{|𝐱_i𝐱_i^{mod}(𝐮)|^2}{\sigma _i^2}$$
(4)
with (in this case) isotropic positional uncertainties $`\sigma _i`$ for each image. The generalization to anisotropic positional uncertainties, such as will be used here for the arcs, as well as fluxes, is immediate. The statistic is minimized with respect to the unobserved source position, $`𝐮`$, and the parameters of the lens model, where the model parameters may be further constrained by measurements of the properties of the lens (position, shape, orientation $`\mathrm{}`$). For non-precision work the problem can be further simplified by simply minimizing the projected separations of the images on the source plane with various flux, magnification or magnification tensor weightings to better mimic the correct image plane statistic (see Kochanek csk-mon (1991)). Most methods using source plane fitting methods are biased to find solutions with high total magnifications and steep radial density profiles.222 Roughly writing $`𝐱_i𝐱_i^{mod}=\delta 𝐱_i`$ for the difference between model and measured image positions $`𝐱_i`$, and $`𝐮𝐮^{mod}=\delta 𝐮`$ one has (from equations 1 and 3) $`\delta 𝐮=M_{ji}^1\delta 𝐱_i`$. So if one minimizes $`|\delta 𝐮|^2/\sigma ^2`$ instead of $`|\delta 𝐱|^2/\sigma ^2`$, increasing $`M`$, the magnification, will tend to improve the fit. The corresponding radial profile will tend to be steeper because steeper profiles generally have higher image magnifications.

### 2.2 Parametrized Models

We model the mass distribution of the lens using a wide variety of parametrized ellipsoidal density distributions. Parametric models of galaxies have the advantages of being simple and physically motivated. The disadvantage of parametric models is that they do not span the complete space of physically realizable mass distributions. We have chosen a set of parametric models broad enough to encompass the range of radial mass distributions considered to represent galaxies. We have, however, restricted the models to ellipsoidal surface densities in an external shear field and thus use a more restrictive range of angular structures than is in principle possible for a real galaxy. The various model families will be discussed in this section and the formulae for their surface densities can be found in the Appendix.

Before discussing the parameterized models, it should be noted that there are non-parametric means of modeling gravitational lenses, principally the linear programming method of Saha & Williams (saha (1997)), where a gridded, positive definite surface density distribution is constructed to fit the data. Such models avoid the overly restrictive assumptions of parametric models and allow elegant solution algorithms for the problem. Coupled with these advantages are two caveats. First, the bulk of the mass in a galaxy, particularly the early-type galaxies that dominate the lens sample, must actually be composed of stars and dark matter with a positive definite, quasi-equilibrium phase space distribution function. This is more restrictive than a positive definite density distribution, which is in turn more restrictive than its integral, a positive definite surface density distribution.333 In an actual phase space distribution function there may be some unoccupied orbits, i.e. there can be some stochasticity due to parts of the galaxy with smaller total particle numbers such as satellites and globular clusters which may be responsible for perturbing the image flux ratios in many lenses, see Mao & Schneider (maosch (1998)). Many density distributions permitted by the Saha & Williams (saha (1997)) method are not physically realizable configurations for galaxies. Thus, even with symmetry and smoothness constraints, requiring a positive definite non-parametric real space density allows too much freedom in the mass model. Second, the current algorithms solve the lens equations on the source plane which can introduce the biases mentioned above. Nonetheless, as this approach captures different features of the matter distribution than the one considered below, applications of this method to B 1933+503 would be very interesting.

The parametric models are characterized by a lens position, $`𝐱_l`$, a mass scale, an axis ratio $`q`$ (or ellipticity $`e=1q`$), and a major axis position angle, $`\theta `$. The major axis density profile is further specified by parameters describing the radial profile of the galaxy such as core radii, break radii and exponents. The models have surface density distributions $`\kappa (m)=\mathrm{\Sigma }(m)/\mathrm{\Sigma }_c=𝑑z\rho (𝐱,z)/\mathrm{\Sigma }_c`$, where $`\rho (𝐱,z)`$ is the density distribution and $`\mathrm{\Sigma }_c`$ is the critical surface density. The surface density $`\kappa `$ is a function of the ellipsoidal coordinate $`m^2=x^2+y^2/q^2`$, where $`x^{}`$ and $`y^{}`$ are Cartesian coordinates aligned with the major and minor axes of the lens. The precise parametric forms $`\kappa (m)`$ of the surface density distributions for the models are given in the Appendix and their properties are discussed below. In order to simplify the comparison between models below, the discussion is restricted to radial properties of a spherical density $`\rho (𝐱,z)\rho (r)`$.

We start with the singular isothermal ellipsoid (SIE) as a standard model. It is both analytically convenient (Kassiola & Kovner kaskov (1993), Kormann et al. kormann (1994), Keeton & Kochanek kkspiral (1998)) and the simplest mass model which is broadly consistent with our general knowledge of the mass distributions in observed galaxies as reviewed above. The singular isothermal sphere has a global $`\rho 1/r^2`$ density profile, which means it has both a steep central density cusp and a globally flat rotation curve. We then explore two general families of ellipsoidal density distributions, varying in radial profile, that contain the SIE as a limit.

The power-law models, whose spherical form is $`\rho (r^2+s^2)^{(\alpha 3)/2}`$, are a simple generalization of the isothermal spheres to include both a central core radius $`s`$ and a variable logarithmic slope $`\alpha `$. The isothermal sphere is the limit with $`s=0`$ and $`\alpha =1`$, and the power-law model sequence includes several “old-fashioned” models of galaxy luminosity profiles such as the modified Hubble profile ($`\alpha =0`$) and the Plummer model ($`\alpha =2`$). This model family was introduced into gravitational lensing because of its relative analytic simplicity and because observational evidence at the time suggested that galaxies possessed finite core radii (see Blandford & Kochanek 1987; Hinshaw & Krauss krauss (1987)). It has remained the most popular “more complicated than isothermal” lens model (e.g. Kochanek csk1654 (1995), Grogin & Narayan grogin (1996), Chae chae (1999), Koopmans & Fassnacht koop (1999)) because of its history and its relative simplicity even as an ellipsoidal model (Barkana bar\_powlaw (1998), Chae, Khersonsky & Turnshek chae\_powlaw (1999)). Unfortunately, we now know that galaxies are characterized by central cusp exponents rather than core radii. Thus the steeper models, $`\alpha <1`$ are unphysical while the shallower models $`\alpha >1`$ are approximations to realistic models with shallow density cusps. We include a sequence of these models with $`\alpha <1`$ only for historical continuity.

The pseudo-Jaffe models are related to the Jaffe (jaffe (1983)) models for early-type galaxies (see e.g., Keeton & Kochanek kkspiral (1998) and references therein). They have $`\rho (r^2+s^2)^1(r^2+a^2)^1`$ for a spherical model, and include both a central core radius $`s`$ and an outer break radius $`a`$. The model has an approximately flat rotation curve for $`sra`$, truncated by the finite central density inside the core radius and by the finite total mass outside the break radius. As we expect the core radius to be zero, we are mainly interested in the model as a simple probe of the truncation of the mass distribution at the break radius. The pseudo-Jaffe models match the SIE model for $`s=0`$ and $`a\mathrm{}`$, and match the $`\alpha =1`$ power-law model in the limit $`sa`$.

We currently lack an implementation of a family of ellipsoidal models with a variable central cusp exponent (the lensing equivalent of the $`\eta `$ models; see Dehnen dehnen (1993), Tremaine et al. tremaine94 (1994)). The analytic lensing properties of such models were studied by Evans & Wilkinson (ew98 (1998)). In order to explore the regime of shallow density cusps we included the NFW (Navarro, Frenk & White nfw (1997)), de Vaucouleurs and $`\alpha >1`$ models in our survey. We are interested in the consequences of the $`\rho 1/r`$ central density cusp of the NFW models, and not in its origins as the characteristic density profile of dark matter halos before modifications due to baryonic matter. The de Vaucouleurs profile also acts like a shallow, cusped density distribution ($`\rho r^{5/4}`$, see Hernquist hernquist (1990)), as well as being the natural constant mass-to-light ratio model for gravitational lenses. The NFW models are characterized by a break radius $`a`$ between their inner $`\rho 1/r`$ cusp and their asymptotic $`\rho 1/r^3`$ profile. The de Vaucouleurs model effective radius $`R_e`$ plays a physical role similar to the NFW break radius. The $`\alpha >1`$ power law models can model the central regions of galaxies with shallow cusps but cannot be good global mass models as they have steadily rising rotation curves even at large radii.

### 2.3 Data

We constrain the models using the relative positions of the lens galaxy and the lensed images and the flux ratios between the images. Here we summarize the sources for the constraints and their limitations. Table 1 presents the constraints and uncertainties derived from the data in Sykes et al. (sykes (1998)), Nair (nair (1998)), Marlow et al (marlow (1999)) and Biggs et al. (biggs (2000)). Figure 1 illustrates the geometry of the system. The source redshift is $`z_s=2.62`$ (Biggs et al. biggs (2000), Norbury et al norbury (2000)) and the lens redshift is $`z_l=0.76`$ (Sykes et al sykes (1998)). The available infrared and optical HST images are too poor to precisely align the radio maps with the HST images or to perform accurate astrometry or surface photometry of the lens galaxy (see Marlow et al. marlow (1999)).

The lensed images (see Figure 1) are composed of three sets of multiply imaged source components. There are four images of the flat spectrum core (image 1, 3, 4 and 6) which have very precise astrometry from the Sykes et al. (sykes (1998)) VLBI observations. There are two images of one steeper spectrum component (images 1a and 8) whose astrometry is relatively poor because of their low fluxes and proximity to the brighter core images. There are four images of a second steeper spectrum component (images 2a, 2b, 5 and 7) but these cannot be naively treated as four images of the same source point because of their extended structure (see Nair nair (1998)). In particular, the peaks of images 5 and 7 are only doubly imaged and do not have counterparts in images 2a and 2b. There is no gap in the arc formed by images 2a and 2b, so they contain only partial images of the source component producing images 5 and 7. Following Nair (nair (1998)) we treat the peaks of images 5 and 7 as a two-image system. We then consider all four images 2a/2b/5/7 as a four-image system, using error ellipses corresponding to the major and minor axes of the images to compensate for our uncertainty in precisely how the extended structure should be modeled.

The flux ratios of the images depend on both frequency and time. To derive the flux ratios of the images we averaged the flux ratio measurements by Sykes et al. (sykes (1998)) and Biggs et al. (biggs (2000)). The differences are probably dominated by systematic differences in the resolution and sampling, temporal variability, and intrinsic systematic problems with flux ratios (see Mao & Schneider maosch (1998)) rather than simple measurement errors. Thus, the problem with flux ratios is always the precision with which they should be imposed. For images 1, 3, 4 and 6, we used 50% errors on the fluxes because the large flux ratio between image 4 and the other images is unlikely to be due to the effects of lensing – generic lens configurations do not generate such flux ratio patterns for point images. For images 1a and 8 we used the reported flux errors. For images 5 and 7 modeled as a two-image system, we used 10% flux errors (which we view as a lower bound on any flux ratio constraints), and for images 2a, 2b, 5 and 7 modeled as a four-image system we used no flux constraints because we are modeling extended rather than point images. We will explore the consequences of changing these estimates of the uncertainties on the results.

No central radio image is detected in any of the radio maps, which is a powerful constraint on the central surface density of the lens. In our standard models we used a detection limit for central images of 0.072 mJy based on the 3$`\sigma `$ detection limit in the 8 GHz observations by Biggs et al. (biggs (2000)). This is very conservative, both because we use a 3$`\sigma `$ limit as a 1$`\sigma `$ limit and because we apply it to the central image for each image system separately rather than to the sum (i.e. if a model produced an 0.072 mJy central image in its model for each of the 1a/8, 1/3/4/6 and 5/7 image systems we will count it as a 3$`\sigma `$ detection instead of a 9$`\sigma `$ detection). We will explore the consequences of using a tighter 0.024 mJy (1$`\sigma `$) constraint. No position constraint was imposed on this undetected central image.

## 3 Results

We now describe the results of modeling the system with the density distributions discussed in §2.1. These consists of three discrete models (SIE, NFW and deV) and two two-parameter model families (the power law models and the pseudo-Jaffe models). All the models use ellipsoidal density distributions and include an independent external shear. Models without the shear terms fit the data so poorly that we do not discuss their properties. The necessity of an independent external shear is expected (see Keeton, Kochanek & Seljak kks (1997)). The external shear required varied from between 2% and 14 % for the best ($`2\sigma `$ models), and the galaxy positions varied by $`10`$ mas. The images will strongly constrain the mass distribution over the region containing images (from $`1.5h^1`$ kpc to $`4.1h^1`$ kpc) with the mass being most tightly constrained near the critical radius of $`2.6h^1`$ kpc, where $`1\stackrel{}{\mathrm{.}}0=5.1h^1`$ kpc (for $`\mathrm{\Omega }_0=0.3`$, $`\lambda _0=0.7`$, and $`H_0=100h`$ km s-1 Mpc-1). Table 2 summarizes the results.

### 3.1 Statistical Summary

The data supply a total of 26 constraints, and the models have $`N_p=7`$ or 8 parameters, leaving us with $`N_{dof}=19`$ or $`18`$ degrees of freedom. We will explore the results using standard results for Gaussian statistics, which we now summarize.

1. Comparing different models: For $`N_{dof}`$ degrees of freedom, the value of our fit statistic should follow a $`\chi ^2`$ distribution with $`N_{dof}`$ degrees of freedom with a mean value of $`N_{dof}`$ and a variance of $`(2N_{dof})^{1/2}`$. For large $`N_{dof}`$ the distribution is asymptotically Gaussian, a limit we can use relatively safely. Thus, for $`N_{dof}=18`$ the mean and variance of the $`\chi ^2`$ distribution are $`18`$ and $`6`$ respectively.
2. Rescaling the statistics: The best models have $`\chi ^211`$ for $`N_{dof}=18`$ which is somewhat low although statistically possible (23% of models with $`N_{dof}=18`$ would have $`\chi ^2<11`$ or $`>25`$ for Gaussian statistics). However, we may have overestimated the uncertainties or degrees of freedom in the model. We can adjust for this by rescaling the models to make the best fit model have $`\chi _{min}^2=N_{dof}`$, which corresponds to shrinking the typical constraint uncertainties by 20% (an entirely plausible correction). Since we do not rescale the fit statistics, we may be overestimating the range of acceptable models and the uncertainties on their parameters.
3. Parameter significance: For a multi-parameter model like the power-law models, the introduction of a parameter (e.g. adding a core radius to the SIE to fit the more general isothermal ellipsoid) is statistically significant if the change in the fit statistic with the addition of the parameter, $`|\mathrm{\Delta }\chi ^2|>f_{lim}\chi ^2/N_{dof}`$, exceeds the results without that parameter by a factor $`f_{lim}=0.5`$, $`3.0`$ or $`8.1`$ for significance levels of 50%, 90% and 99% compared to the improvement expected for a random variable. These limits are based on the F-test for adding one variable to a model with $`N_{dof}=18`$ or $`19`$ degrees of freedom.
4. Parameter estimation: For a multi-parameter model, parameter ranges are estimated from the change in the fit statistics relative to the best model, $`\mathrm{\Delta }\chi ^2=\chi ^2\chi _{min}^2`$ where the parameter range for one variable is acceptable at the 1$`\sigma `$, 90%, 2$`\sigma `$ and 99% confidence level for $`\mathrm{\Delta }\chi ^2=1.0`$, $`2.71`$, $`4.00`$ and $`6.63`$ respectively. Upper bounds and parameter ranges quoted below will be $`1\sigma `$ bounds ($`\mathrm{\Delta }\chi ^2=1`$) unless otherwise noted.

### 3.2 Model Properties

We start with the standard SIE model with external shear, which has $`\chi ^2/N_{dof}=12.0/19`$. The model results are compared to the constraints in Figure 1 and are remarkably good, given the incredible simplicity of the model. The SIE model is a particular example of both the power law and pseudo-Jaffe two-parameter models, which we consider next. Adding a core radius to the SIE, making a general isothermal ellipsoid, slightly improves the fit to $`\chi ^2/N_{dof}=11.4/18`$ but the improvement is not statistically significant (66% for the F-test). The formal core radius is $`140h^1`$ pc, but the more useful result is an upper limit of $`260h^1`$ pc. The fit is not improved by adding an outer break radius $`a`$ to the isothermal density profile (the pseudo-Jaffe models). The value of $`\chi ^2`$ monotonically rises as the break radius is reduced from the SIE limit ($`a\mathrm{}`$), leading to a lower bound of $`a>1\stackrel{}{\mathrm{.}}47.1h^1`$ kpc. The core radius of the pseudo-Jaffe model plays the same role as that of the isothermal ellipsoid – a small core radius slightly, but insignificantly, improves the fit.

The traditionally favored power-law models, with $`\alpha <1`$, a finite core radius, and a more rapidly declining asymptotic density distribution than the isothermal models, have increasingly large values for $`\chi ^2`$ as the mass distribution becomes more centrally concentrated (see Figure 2). These models require large core radii to fit the data (e.g. $`(653_{97}^{+92})h^1`$ pc for the modified Hubble profile, see Figure 3). In addition, for $`\alpha <0`$ models, external shears of $`>20\%`$ are required, increasing as $`\alpha `$ decreases. These large external shears make the interpretation of the parameters questionable–consequently we have cut off our analysis at the $`\alpha =1`$ model (shear $`23\%`$); the $`\alpha <1`$ models have $`\chi ^2>35`$ and external shears of $`>27\%`$. The power law models are consequently unacceptable unless nearly isothermal ($`\alpha 1`$), as we might expect from our current understanding of the structure of galaxies (see §1). The poor fit of the models is not solely due to the restriction on the visibility of the central image, as these models remain unacceptable when we weaken or eliminate the constraints on the image flux ratios (see below).

The fit does improve significantly for power-law models with $`\alpha >1`$, where the models are mimicking galaxies with softer central density cusps than the SIE model. For these models the core radius must be small, and for simplicity we only show the results with the core radius fixed.444For example, the $`\alpha =1.1`$ power-law model has the same $`\chi ^2`$ for an optimized core radius ($`0\stackrel{}{\mathrm{.}}0002`$) and for the core radius fixed as in Table 2. The best fit is found for $`\alpha =1.14_{0.17}^{+0.11}`$ which corresponds to a density cusp with $`\rho r^{1.86_{0.17}^{+0.11}}`$ The exponent $`\alpha `$, unlike the core radius $`s`$, is a statistically significant new parameter (82% for the F-test). Unfortunately, the $`\alpha >1`$ models are poor global models because of their steadily rising rotation curves.

The NFW and de Vaucouleurs (deV) models provide examples of models with shallow central cusps (near $`\rho 1/r`$) and are reasonable global models for the mass distribution. (The $`\alpha >1`$ models cannot produce such shallow cusps because as $`\alpha 2`$ they become constant surface density sheets rather than $`\rho 1/r`$ cusps because of how the limits are ordered.) Figure 3 shows the variation in $`\chi ^2`$ with the de Vaucouleurs effective radius $`R_e`$, the NFW break radius $`a`$, and, for comparison, the core radius $`s`$ of the isothermal ellipsoid ($`\alpha =1`$) and the pseudo-Jaffe model break radius $`a`$. The NFW model is barely acceptable at 2$`\sigma `$ significance, with a break radius of $`a=(1.53_{0.26}^{+0.36})h^1`$ kpc. The more cuspy de Vaucouleurs models produce acceptable fits for an effective radius of $`R_e=(14_4^{+9})h^1`$ kpc. This is much larger than is realistic for a constant mass-to-light ratio model of the lens galaxy.555While the HST images are not adequate to fit a photometric model, the optical and IR emission is confined to a region comparable in size to the lens, so a real constant mass-to-light ratio model would have $`R_e<1\stackrel{}{\mathrm{.}}0`$ ($`5h^1`$ kpc).

The uncertainties we used for the image flux ratios are somewhat arbitrary (see §2.3), so we explored the effects of changing these limits. We recalculated all the models with 50% flux uncertainties for all images ($`\chi _{min}^2/N_{dof}=9.0/19`$), no flux constraints at all ($`\chi _{min}^2/N_{dof}=6.7/10`$), and using the larger of the quoted errors in the flux measurements or 10% of the flux ($`\chi _{min}^2/N_{dof}=61.5/19`$). In all three of these cases the best-fit model was the SIE, or equivalently the $`\alpha =1`$ isothermal ellipsoid with the statistically insignificant addition of a core radius. The $`\alpha =1.14`$ best-fit model for our standard flux ratio uncertainties remains within 1$`\sigma `$ of the best fitting models for two cases, and had a $`\mathrm{\Delta }\chi ^2=1.4`$ for the case where 10% flux limits are used. In the case where no flux limits are imposed the $`\alpha <1`$ models are still disfavored, with $`\alpha >0.5(1\sigma ),0.4(2\sigma )`$. The overall statistical ordering of the models is essentially unaffected by changes in the choice for the flux ratio uncertainties. We also explored effects of changing the limits on the flux of any central images, discussed below.

The central cusp exponent or core radius of the distribution determines the visibility of the central or “odd” image (Narasimha, Subramanian & Chitre nara86 (1986)), whose flux decreases as the central surface density or concentration of the lens increases. Thus, the limits on the core radius depend on the treatment of the central image. Figure 4 illustrates the consequences of tightening the limit on the flux of central images from 0.072 mJy (3$`\sigma `$ in the 8 GHz map) to 0.024 mJy (1$`\sigma `$) for the isothermal ellipsoid ($`\alpha =1`$), the $`\alpha =1`$ power-law, and the de Vaucouleurs models. The scale radii are forced to decrease, and for the $`\alpha =1`$ power-law model and the de Vaucouleurs model the $`\chi ^2`$ statistics rise significantly. The change in $`\chi ^2`$ for models with low $`\chi ^2`$ is generally small ($`\mathrm{\Delta }\chi ^2<1`$). Mao, Witt & Koopmans (2000) explored the effects of adding a central black hole to the lens galaxy as an alternate means of suppressing the flux of central images. While we do not expect the addition of a central black hole to modify the global trends discussed below (as model comparisons were essentially unchanged when the flux constraints were completely removed), it would be interesting to include central black holes in future analysis.

Finally, we explored the effects of further generalizing the models. In particular we were concerned that the limitation to ellipsoidal models might force an undesirable correlation between the monopole and quadrupole structures of the lens galaxy. We tested this by modeling the galaxy as a combination of two pseudo-Jaffe models where the inner cutoff radius of one equaled the outer core radius of the other (see Keeton et al. keeton2k (2000)) and then optimized the relative masses, ellipticities and orientations of the two components. The best fit models had $`\chi ^21011`$ and 2 or 3 fewer degrees of freedom, so the improvements in the fit were not statistically significant given the reduction in the number of degrees of freedom. Modeling this system as a combination of two other models (Keeton 2000b ) produced similar behavior in the monopole and quadrupole structures discussed below.

### 3.3 Global trends

Parametric models make it difficult to see the global trends implied by the range of acceptable models or to compare the physical properties of the individual models. We can illustrate the similarities and differences of the models by looking at the monopole and quadrupole deflection profiles of the models (see e.g. Kochanek csk-mon (1991) for definitions). For a potential such that $`^2\varphi =2\kappa (r,\theta )=2\mathrm{\Sigma }(r,\theta )/\mathrm{\Sigma }_c`$, the monopole deflection of the lens is

$$\alpha _0(r)=\frac{1}{r}_0^rr^{}𝑑r^{}_0^{2\pi }𝑑\theta ^2\varphi (r)=\frac{2}{r}_0^rr^{}𝑑r^{}_0^{2\pi }𝑑\theta \frac{\mathrm{\Sigma }(r,\theta )}{\mathrm{\Sigma }_c}=\frac{2M(<r)}{\mathrm{\Sigma }_cr}$$
(5)
where $`M(<r)`$ is the mass enclosed by the circle and the coordinates are centered on the lens galaxy. The rotation curve of the galaxy is roughly proportional to the square root of the deflection monopole (projection effects modify this slightly). Generic deflection or rotation curve profiles have three radial regimes. At small radii there is a core region where the deflection can go to zero either because of a finite core radius or the presence of a shallow density cusp. The size of the region is set either by the core radius or the break radius where the steep outer density cusp is matched to the shallower inner cusp. In the inner region the rotation curve rises, and can become flat for a while. At large radius the distribution must have a break beyond which the rotation curve becomes Keplerian (or near Keplerian for the NFW model).

Figure 5 shows the monopole as a function of radius for the models from Table 2. As is typical of lens models (see Kochanek csk-mon (1991)), all models agree very precisely on the mass enclosed by the mean critical line of the lens. For a critical surface density of $`\mathrm{\Sigma }_c=0.62h`$ g cm-2 ($`\mathrm{\Omega }_0=0.3`$, $`\lambda _0=0.7`$), the mass inside $`r=0\stackrel{}{\mathrm{.}}5`$ is $`6.0\times 10^{11}h^1M_{}`$ independent of the parametric model. The statistically acceptable models have slightly rising or falling deflection profiles, and since $`\alpha _0v_c^2`$ the rotation curve of the acceptable models varies by $`<`$ 8.5% from the average between the inner and outer images ($`1.7h^1`$ kpc to $`4.1h^1`$ kpc). The model with the least flat rotation curve ($`\alpha =1.3`$) has a core of less than $`0\stackrel{}{\mathrm{.}}0001`$, and so projection effects should not alter this limit significantly (for the spherically symmetric case the change is less than $`0.01\%`$). The models with a slightly falling rotation curve are the pseudo-Jaffe models, and the models with the rising rotation curves are the $`\alpha >1`$ shallow-cusp power law models. The worst fitting models from Table 1 tend to have rapidly rising and then falling rotation curves as shown on the right in figure 5.

The second simple physical property of the models we can compare is the quadrupole of the lens. We define the quadrupole by the appropriate multipole moment of the angular deflection,

$$\left(\begin{array}{c}\hfill \gamma \mathrm{cos}2\theta _\gamma \\ \hfill \gamma \mathrm{sin}2\theta _\gamma \end{array}\right)=\frac{1}{\pi r^2}𝑑\theta \left[\frac{}{\theta }\varphi \right]\left(\begin{array}{c}\hfill \mathrm{sin}2\theta \\ \hfill \mathrm{cos}2\theta \end{array}\right),$$
(6)
including both the lens galaxy and the external shear in the potential. We have defined the quadrupole such that the quadrupole of an SIE model is proportional to $`r^1`$ and the quadrupole of an external shear is independent of radius. An external shear $`(\gamma /2)\mathrm{cos}(2\theta 2\theta _\gamma )`$ produces a quadrupole moment $`\gamma `$ at orientation $`\theta _\gamma `$. Figure 6 shows the amplitude and major axis PA of the quadrupole moment as a function of radius for the statistically acceptable models.666Note that the PA is related to $`(x,y)`$ via $`x=r\mathrm{sin}\theta ,y=r\mathrm{cos}\theta `$, while the angles in the equations obey the more standard definition $`x=r\mathrm{cos}\theta ,y=r\mathrm{sin}\theta `$. The measured PA of the major axis of the lens, $`40^{}\pm 5^{}`$ (Sykes et al sykes (1998)) is shown for comparison. All the sucessful models show a similar decline and flattening of the strength of the quadrupole combined with a small $`3^{}`$ swing in orientation. The poorly fitting power law models also have very different quadrupole profiles from the acceptable models. When we tested the composite models, where the properties of the monopole and quadrupole are largely decoupled, the monopole and quadrupole deflection profiles of the best two component models were little changed from the acceptable single profile models. Some examples are included in Figure 6.

Although B 1933+507 has yet to show the variability needed to measure a time delay (Biggs et al. biggs (2000)), the constraints on the mass distribution are significantly better than for the lenses with delay measurements. Figure 7 shows the variation in the predicted longest delay between the compact radio components (between images 1 and 6) as a function of the $`\chi ^2`$ statistic for the various models. The delay depends on the outer cutoff of the mass distribution, with the delay rising as the mass distribution is truncated at smaller radii and the depth of the central potential decreases. This is a common feature of lens time delays (see Keeton & Kochanek kk (1997), Impey et al. impey (1998), Koopmans & Fassnacht koop (1999), Witt, Mao & Keeton wmk00 (2000)). Formally we find a delay range of $`(10.6_{1.1}^{+2.4})h^1`$ days. As seems to be typical of lens systems (see Keeton & Kochanek kk (1997), Impey et al. impey (1998)), the ratios of the delays between the images depend weakly on the model. For example, the ratio of the delay between images 1 and 4 to the delay between images 1 and 6 changes by only $``$ 3% ($`<0.4h^1`$ days) for the range of the statistically acceptable models. Models with central cusps which are more realistic than the $`\alpha >1`$ power law models will have an outer break radius similar to the pseudo-Jaffe models. As the models are varied by decreasing the break radius the time delays will rise.

## 4 Conclusions and Prospects

After fitting a wide range of parameterized, ellipsoidal mass distributions to the B 1933+503 system, we conclude that the mass distribution must have an approximately flat rotation curve with a cuspy central density distribution. Over the range spanned by our acceptable models the rotation curve deviates from the average by less than $``$8.5% between the inner and outer images ($`1.7h^1`$ kpc to $`4.1h^1`$ kpc). The central cusp needs to be quasi-isothermal ($`\rho r^2`$) and models with shallow cusps ($`\rho 1/r`$) fail to fit the data. Traditional power law models with finite central core radii are inconsistent with the data, as we would expect from our current understanding of galaxy mass distributions. Both the monopole and quadrupole deflection profiles of the acceptable models are similar and stable to increasing the complexity of the models further. The historically popular power-law models for gravitational lenses should be abandoned in favor of the lensing equivalents of the $`\eta `$-models (Dehnen dehnen (1993), Tremaine et al. tremaine94 (1994)). These models are parametrized by a central cusp exponent and break radius between an interior cusp and a more rapidly declining outer density profile, rather than a core radius and the slope of the asymptotic density profile.

So far, all lenses which have been analyzed to determine the radial mass distribution of the lens prefer mass distributions corresponding to nearly flat rotation curves. This is true of MG 1654+1346 (Kochanek csk1654 (1995)), MG 1131+0456 (Chen et al. ckh95 (1995)), Q 0957+561 (Keeton et al. keeton2k (2000)), PG 1115+080 (Kochanek et al. kkm (2000)) and B 1933+507 (this paper). No lens with an image distribution suitable for determining the radial mass profile has ever found results which are grossly inconsistent with a nearly flat rotation curve. Similar results are also found for observations of nearby early-type galaxies with X-ray halos (Fabbiano fabbiano (1989)) and modern stellar dynamical studies (e.g. Rix et al. rixetal (1997)).

The predicted time delays of B 1933+503 are relatively well constrained at $`(10.6_{1.1}^{+2.4})h^1`$ days between images 1 and 6 over the wide range of possible mass distributions we explored. A $``$22% modeling uncertainty is encouraging because there are excellent prospects for further improving the model constraints. Improved radio maps with a better quantitative understanding of the mapping between images 2a, 2b, 5 and 7, astrometry of images 1a and 8, and tighter constraints on the flux of any central images would all be useful. The time delay ratios between the images can constrain the models if measured to 3% or better. Deeper infrared images of the system, particularly with the Hubble Space Telescope, could enormously improve the relative astrometry between the lens galaxy and the radio sources and accurately determine the structure of the lens galaxy (ellipticity, orientation $`\mathrm{}`$). Most importantly, the existing HST images strongly suggest the presence of lensed images of the radio source host galaxy. Deep imaging to measure the shape of the Einstein ring image of the host would dramatically improve the constraints (see Kochanek, Keeton & McLeod kkm (2000)).

Acknowledgements: CSK was supported by NASA Astrophysics Theory Program grants NAG5-8831 and NAG5-9265. Support for the CASTLES project was provided by NASA through grant numbers GO-7495, GO-7887, and GO-8175 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. JDC thanks M. White for conversations and the ITP for hospitality. JDC was supported in part by the National Science Foundation under Grant No. AST-0074728 and while at the ITP under Grant No. PHY94-07194. We thank the referee for helpful suggestions and J. Munoz for discussions on the input constraints and comments on the draft.

## Appendix A Summary of Models

Here we present the analytic expressions for the model surface densities in terms of the ellipsoidal coordinate $`m`$, defined as $`m^2=x^2+y^2/q^2`$, where $`x^{}`$ and $`y^{}`$ are Cartesian coordinates aligned with the major and minor axes of the lens. The singular isothermal ellipsoid (SIE) is defined by

$$\kappa _{SIE}(m)=\frac{1}{2}\frac{b}{|m|}\mathrm{SIE},$$
(A1)
and the power-law models are defined by

$$\kappa _\alpha (m)=\frac{1}{2}b^{2\alpha }\left(s^2+m^2\right)^{\alpha /21}\mathrm{power}\mathrm{law}.$$
(A2)
The power-law model with $`\alpha =1`$ and $`s=0`$ is the SIE model. The pseudo-Jaffe model is the difference of two softened isothermal ellipsoids ($`\alpha =1`$ with a finite core radius), where

$$\kappa _{PJ}(m)=\frac{1}{2}b\left[\left(s^2+m^2\right)^{1/2}\left(a^2+m^2\right)^{1/2}\right]\mathrm{pseudo}\mathrm{Jaffe}$$
(A3)
with $`sa`$. The core radius of the profile is $`s`$ and the break radius is $`a`$. It matches the SIE model for $`s0`$ and $`a\mathrm{}`$, the general $`\alpha =1`$ power-law model for $`a\mathrm{}`$ and the $`\alpha =1`$ power law model in the limit that $`sa`$ with $`b(as)`$ kept constant. The de Vaucouleurs model is defined by

$$\kappa _{dV}(m)=\frac{b}{2𝒩R_e}\mathrm{exp}\left[k\left(\frac{m}{R_e}\right)^{1/4}\right]\mathrm{deVaucouleurs}$$
(A4)
where

$$k=7.67\mathrm{and}𝒩=_0^{\mathrm{}}ve^{kv^{1/4}}𝑑v=1.683\times 10^3$$
(A5)
and the NFW model by

$$\kappa _{NFW}(m)=2\kappa _s\frac{1(m/a)}{(m/a)^21}\mathrm{NFW}$$
(A6)
where

$$(m/a)=\{\begin{array}{cc}\frac{1}{\sqrt{m^2/a^21}}\mathrm{tan}^1\sqrt{m^2/a^21}\hfill & m/a>1\hfill \\ \frac{1}{\sqrt{1m^2/a^2}}\mathrm{tanh}^1\sqrt{1m^2/a^2}\hfill & m/a<1\hfill \\ 1\hfill & m/a=1\hfill \end{array}$$
(A7)