DISENTANGLING BARYONS AND DARK MATTER IN THE SPIRAL GRAVITATIONAL LENS B1933+503

We observed B1933+503 with the global VLBI array on 1998 November 27 at 1.7 GHz with a bandwidth of 16 MHz. We used 17 telescopes with 10 from the Very Long Baseline Array (VLBA) and 7 from the European VLBI Network (EVN). We adopt the center of component 4 of the lensed images as the phase center for the observations. Observations were conducted on a cycle of 6.5 minutes, with 1.5 minutes on a phase calibrator (B1954+51 = J1955+5131) followed by 5 minutes on the target source B1933+503. The exception was the Lovell telescope, which has a slower slew rate and for which every other observation of the phase calibrator was omitted, yielding a 1.5 minute + 11.5 minute cycle on the phase reference and target. The total observing time was nine hours, providing the good u–v coverage required for high dynamic range imaging and image fidelity. A single observation of 3C345 was obtained for fringe finding and flux calibration.

We reduce the data with the Astronomical Image Processing System¹⁴ package using standard procedures. The images are iteratively CLEANed and self-calibrated (phase-only), before a single amplitude self-calibration solution is performed, to remove residual phase and amplitude errors. The final maps are produced using natural weighting of the visibilities and with a 10 Mλ taper to increase the sensitivity to the extended emission. The natural-weighting image has an rms noise level of 50 μJy beam⁻¹ with a beam size of 5.7 × 2.7 mas at a position angle (P.A.) of −322. The tapered image has an rms of 55 μJy beam⁻¹ with a beam size of 20.1 × 16.5 mas at a P.A. of 6055. Figure 1 shows the final maps for each of the components in B1933+503, except for component 6, which is not detected.

Zoom In Zoom Out Reset image size

The radio source being lensed has a compact core with two extended radio lobes on opposite sides of the core (Sykes et al. 1998; Nair 1998). The radio spectral indices make the identification of the images of the compact core unambiguous (Sykes et al. 1998), and the flux densities of the core components showed little time variability (Biggs et al. 2000). Components 1, 3, 4, and 6 correspond to the images of the compact core, components 1a and 8 are the images of one of the lobes, and components 2, 5, and 7 correspond to the images of the other radio lobe. With component 2 counting as two (merging) images, B1933+503 is a spectacular 10 image radio lens.

We use components 1 and 4 to align the global VLBI data to the previous radio observations (with the Very Large Array, MERLIN, and VLBA) in Sykes et al. (1998) and Marlow et al. (1999). Table 1 lists the image positions for each of the components in the previous observations compiled by Cohn et al. (2001) and in the global VLBI data. The merging pair of images are denoted by components 2a and 2b. In the global VLBI data, several components (1, 4, 5, and 2a) have multiple intensity peaks, but some of their corresponding images (3 and 7) have only single peaks due to, e.g., scatter-broadening during propagation through the disk of the lens galaxy (e.g., Marlow et al. 1999; Norbury 2002). We adopt the flux-weighted average of the peak positions as the image positions and estimate the uncertainties from the separation of the peaks. We also set the positional uncertainty of component 7 and 2b to that of their counter image (component 5). For components that are spatially extended (1a and 8), the uncertainty is set to the geometric mean of the major and minor axis of the beam.

Cohn et al. (2001) also compiled flux ratios from Sykes et al. (1998), Nair (1998), and Biggs et al. (2000). We do not list or use the flux ratios because they could be systematically biased due to scattering and substructures in the lens. In fact, Kochanek & Dalal (2004) showed that the flux ratios in B1933+503 are anomalous due to substructures, and the constraints in B1933+503 lead to a relatively smooth, elliptical macro mass model. Furthermore, Cohn et al. (2001) found that the exclusion of flux ratios had little effect on the main results of their mass models.

We observed B1933+503 with the adaptive-optics (AO) system at the Keck II telescope on 2005 July 31. We used the Kp filter (at 2.2 μm) on the Near Infrared Camera 2 (NIRC2) and took 33 dithered frames with a field of view of 10'' × 10'' that are each six co-adds of 30 s. The dithering allows good sky subtraction in the reduction. The natural tip-tilt reference star of magnitude R ∼ 16.4 located at ∼213 southeast of B1933+503 and a sodium laser guide star were used to correct for atmospheric turbulence.

We use an IRAF-based algorithm to subtract the sky, remove bad pixels and cosmic rays, and produce weight maps for each frame. We use the weight maps and the Drizzle package (Fruchter & Hook 2002) to co-add the images together. We show in Figure 2 the drizzled NIRC2 image of the lens system. By comparing the positions of the three stars in the NIRC2 image to the corresponding stars in previous HST observations of the system (Section 2.6), we determine the pixel scale of the NIRC2 image to be 9.964 mas pixel⁻¹. The closest field star to B1933+503 (∼29 away) in the combined NIRC2 image has an FWHM of ∼57 mas.

Zoom In Zoom Out Reset image size

In addition to the high-resolution NIRC2 image, we observed B1933+503 on 1996 July 31 with the Near Infrared Camera (NIRC; Matthews & Soifer 1994) at the Keck Observatory. The NIRC detector is 256 pixels on a side, with a pixel scale of 015 pixel⁻¹. Thus, the field of view of the camera is 384 on a side, i.e., larger than the NIRC2 narrow camera field of view. The system was observed in both J and K bands, but the image quality of the J-band data was too poor to be useful. The K-band data consist of 59 exposures, each with an exposure time of one minute (five co-adds of 12 s each). The data were processed in a standard fashion, including steps to subtract the dark current, flatten the images, and subtract the sky. For each of the images, the flat field and sky frames were constructed from the frames observed directly before and after the image. The processed images were aligned by measuring the position of a star that appeared in each frame and then were co-added to produce the final image. The photometric zero point of the NIRC data was determined through a comparison to the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), which has K-band magnitudes for the two brightest stars in the image. The uncertainty in the photometric zero point is estimated to be 0.15 by comparing the magnitudes of the two stars. This incorporates the uncertainty due to the difference in throughput of the 2MASS and NIRC K-band filters, which we estimate to be ≲ 0.03 by comparing the magnitude difference between the filters for several stars.

We observed B1933+503 using the Echelle Spectrograph and Imager (ESI; Sheinis et al. 2002) instrument at the Keck Observatory on 2002 June 6. The slit of width 125 was aligned along the major axis of the lens galaxy (P.A. = 138°). The standard echellette mode yielded a spectral resolution with Gaussian width of 30 km s⁻¹. Five exposures of 1800 s each were taken under good seeing conditions of 06.

We reduce the data using the software EASI2D developed by D. J. Sand and T. Treu (Sand et al. 2004). We see prominent emission lines Hβ (λ4861), [O ii] (λ3726, λ3729), and [O iii] (λ4959, λ5007) from the lens galaxy in the spectra. For each line, we set the center of the galaxy to the brightest pixel in the two-dimensional spectrum since measuring the center to much less than a pixel is difficult due to the presence of seeing. We then bin the spectrum in the spatial direction by a factor of three to increase the signal. Thus, the center of the galaxy should be well within the central binned spatial pixel of size ∼05. Any small offsets in the centroid are taken into account in the modeling in Section 7.3. The systemic velocity is difficult to measure from the data directly and is determined in the modeling (Section 7.3). Figure 3 is the rotation curve based on these lines.

Zoom In Zoom Out Reset image size

We carried out infrared spectroscopic observations of the B1933+503 lensed source on 2005 September 9 with the Near Infrared Spectrometer (NIRSPEC; McLean et al. 1998) on the Keck II telescope. The data were taken through the NIRSPEC 6 and 7 filters (N6 and N7, respectively), which gave an approximate wavelength coverage of 1.56 to 2.02 μm and 2.08 to 2.52 μm, respectively. Each observation consisted of 4 × 300 s exposures that used dithering along the slit to improve the removal of the sky background during the reduction stage. The total exposure time was one hour through each filter. The slit was put at a position angle of 2267 to cover the two strong infrared components from the lensed source. The standard star HD162208 was also observed four times each in N6 and N7 for calibration.

We reduce the data with a Python-based pipeline that removes cosmic rays, subtracts the sky, wavelength calibrates using the atmospheric skylines, and extracts a one-dimensional spectrum. Furthermore, we use the standard star HD162208 to correct for the response of the spectrograph. The reduced spectra for both filters are presented in the bottom panel of Figure 4, and the atmospheric transmission is shown in the top panel. The spectra are smoothed using a 7 pixel moving average with each point being weighted by the inverse variance associated with it. We identify a strong emission line in the N6 spectrum (bottom panel, left), but do not see any spectral features in the N7 spectrum (bottom panel, right). Therefore, we believe that the detected emission line is likely Hα blended with two [N ii] emission lines, corresponding to a source redshift of z_s = 1.71 ± 0.01. The uncertainty in the redshift is conservative and accounts for the blending of the lines and also the contamination by atmospheric lines on both sides of the emission line.

Zoom In Zoom Out Reset image size

We note that the strong emission line detected here has been previously reported by Biggs et al. (2000), based on an unpublished spectrum that was taken with the United Kingdom Infrared Telescope. Then the line was interpreted to be Hβ at redshift 2.62 due to a second spectral feature that was believed to be Hα in the K band. Our much better spectral resolution and higher sensitivity data do not detect the second emission line in N7, which would certainly have been detected if it were Hα given the relative flux of the supposed Hβ line in N6 and the noise level. We therefore rule out a redshift of 2.62 for the lensed source.

Archival HST images of B1933+503 in the F160W filter with the Near Infrared Camera and Multi-Object Spectrometer (NICMOS; Proposal ID: 9744; PI: Kochanek) and the F555W and F814W filters with the Wide Field and Planetary Camera 2 (WFPC2; Proposal ID: 9133; PI: Falco) are available. We used MultiDrizzle¹⁵ to combine the exposures in each filter and correct for geometric distortion. The F555W data were discarded due to the low signal-to-noise ratio (S/N). We oversampled the F814W and F160W images with MultiDrizzle to the same pixel scale as the NIRC2 image, and show the color image composed of the three filters in Figure 5. The images of the lensed source in the plane of the lens galaxy (corresponding to radio components 3 and 6) suffer from dust extinction.

Zoom In Zoom Out Reset image size

We use the stars in the WFPC2 F814W images to obtain the pixel scale of the NIRC2 image in Section 2.2 and the alignment of radio and NIRC2 images in Section 3. Furthermore, we use the F814W and F160W photometries in Section 4.3 to estimate the stellar mass of the lens galaxy (that incorporates the effects of dust) in Section 9.4.

A model of the point-spread function (PSF) is needed to extract the intrinsic light profile of the lens galaxy without atmospheric and instrumental blurring. We use the star that is ∼29 northwest of B1933+503 in the observed field as a model of the PSF. Previous works have shown that field stars serve as good PSF models, especially for spatially extended objects (e.g., Marshall et al. 2007; Suyu et al. 2009).

To model the lens galaxy light, we use the Galfit package (Peng et al. 2002). We mask out the lensing arcs and the spiral arm-like features in the southeast region of the galaxy. Optionally, we also mask out the central region of the galaxy. We find that the galaxy is well described by two exponential disk profiles with a common centroid (or a single exponential disk profile, if the central region was masked out). Specifically, we employ Sérsic profiles with index n_sersic ≡ 1 that correspond to exponential disks and also fit a uniform background for the sky. One of the Sérsic profiles corresponds to the disk of the galaxy and the other is centrally localized with a small effective radius (∼005).¹⁶ We identify this latter component as a bulge, although it could also be a bar given the limited resolution. We estimate the uncertainties on the Sérsic parameters based on differences in the best-fit parameter values for different choices of masks. This systematic uncertainty dominates the statistical uncertainty of the fit. For the galaxy centroid, we include the systematic uncertainty of 5 mas (from the alignment of NIRC2 and radio data) which dominates the overall positional uncertainty. Table 2 lists the best-fit values for the mask containing the lensed arcs and spiral arm features.

We see in Figure 6 that the two-component Sérsic model reproduces the overall structure of the lens galaxy light, and the background fit yields uniform sky residuals. The reduced χ² is 1.03 in the fitting region. The residuals show, apart from the three strong lensing arcs (corresponding to radio components 1, 4, and 6 in Figure 2), some small-scale features that could be spiral arms or tidal features in the lens galaxy. In Sections 7.3 and 8.2, we account for the effects of these residual features in the mass modeling by adding systematic uncertainties to the line-of-sight velocities and the radio image positions.

Zoom In Zoom Out Reset image size

The photometry of the lens is required for estimating the stellar mass from SPS models in Section 9.4. To obtain the integrated K-band magnitudes for the disk and the bulge based on the exponential profiles in the previous section, we first determine the photometric zero point of the NIRC2 image by calibrating it with the low-resolution NIRC K-band image. We then integrate the model light profiles and list in Table 3 the K-band magnitudes of the disk and the bulge, where the estimated uncertainty comes from calibration, difference in throughput of NIRC2 Kp and NIRC K filters, and variations between different arc masks.

The HST images have significantly lower S/N than the NIRC2 image. Therefore, we use the structural parameters (axis ratio, position angle, and effective radius) of the exponential profiles from the NIRC2 image (listed in Table 2) to obtain the integrated magnitudes in F814W and F160W with Galfit. We use the nearest field star to approximate the PSF. The separation of the arcs and the lens light is more difficult in these low S/N images, and we conservatively adopt an uncertainty of 0.3 mag from various choices in the arc masks. The integrated magnitudes for the disk and the bulge are listed in Table 3.

We assume that the disk of stars and gas in the galaxy is thin and circularly symmetric. Furthermore, we assume that there is a constant M/L throughout the disk. Hence, the exponential profile for the light in Section 4.2 implies that the profile for the elliptical surface mass density of the projected (inclined) disk is

$$\begin{equation} \Sigma _{\rm d}^{\rm P}(R^{\prime }) = \Sigma _{\rm d,0}\,\exp [-R^{\prime }/R_{\rm d} ], \end{equation} \tag{ 1 }$$

where Σ_{d, 0} is the normalization of the disk (set by the M/L), R_d is the scale radius of the disk, and R' is given in terms of the coordinates x', y' along the major/minor axes centered on the galaxy by

$$\begin{equation} R^{\prime 2} = x^{\prime 2} + (y^{\prime }/q_{\rm d})^2. \end{equation} \tag{ 2 }$$

The axis ratio q_d is related to the inclination angle i of the disk by

$$\begin{equation} q_{\rm d}=\cos i, \end{equation} \tag{ 3 }$$

where i = 0° corresponds to a face-on disk.

In terms of the radial coordinates in the plane of the disk, $R=\sqrt{x^2+y^2}$, the surface mass density is

$$\begin{equation} \Sigma _{\rm d}(R) = \Sigma _{\rm d,0}\,q_{\rm d}\,\exp \left[ -R/R_{\rm d}\right], \end{equation} \tag{ 4 }$$

where the extra factor of q_d in the normalization ensures that an inclined disk and a face-on disk have the same total mass. The scale radius R_d is related to the effective radius (the half-light radius in Table 2, which is also the half-mass radius with a constant M/L) by R_d = R_e/1.678.

For gravitational lensing, the quantity of interest is the dimensionless surface mass density,

$$\begin{eqnarray} \kappa _{\rm d}(R^{\prime })&=&\Sigma _{\rm d}^{\rm P}(R^{\prime })/\Sigma _{\rm crit} \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &=& \kappa _{\rm d,0}\,\exp [-R^{\prime }/R_{\rm d} ], \end{eqnarray} \tag{ 6 }$$

where

$$\begin{equation} \Sigma _{\rm crit} = \frac{c_{\rm l}^2}{4\pi G}\frac{D_{\rm s}}{D_{\rm d} D_{\rm ds}}, \end{equation} \tag{ 7 }$$

κ_{d, 0} is the disk strength, c_l is the speed of light, and G is the gravitational constant. The distances D_d, D_s, and D_ds are angular diameter distances to the lens, to the source, and between the lens and the source, respectively.

The galaxy light fitting suggests that we can model the bulge as a point mass, since the effective radius of the corresponding Sérsic profile is very small compared to the radial range of positions spanned by the lensed images (∼03–08) and most of the points in the rotation curve.

We assume an NFW profile (Navarro et al. 1996) with a triaxial shape (Jing & Suto 2002) for the dark matter halo. The three-dimensional density distribution is given by

$$\begin{equation} \rho _{\rm h}(r) = \frac{\rho _{\rm h,0}}{(r/r_{\rm h,0})(1+r/r_{\rm h,0})^2}, \end{equation} \tag{ 8 }$$

where

$$\begin{equation} r^2 = c^2\left(\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2}\right), \quad a\le b \le c. \end{equation} \tag{ 9 }$$

The parameters a, b, and c describe the triaxial shape of the halo. The orientation of the dark matter halo as seen by a distant observer can be described by two angles ϑ and φ. Following Oguri et al. (2003), we choose (ϑ, φ) as the polar angle of the observer's line-of-sight direction in the halo coordinate system (x, y, z) (see Figure 1 of Oguri et al. 2003). With this definition, Equations (1) and (2) of Oguri et al. (2003) relate the coordinates in the frame of the halo (x, y, z) to the coordinates of the distant observer (x', y', z').

The two-dimensional surface mass density of the projected triaxial halo is elliptical, and the axis ratio (q_h) and position angle (ϕ_h) of the ellipse can be calculated given the values for a/c, b/c, ϑ, and φ (Oguri et al. 2003). We use the Einstein radius for the corresponding spherical mass model, R_{h, E}, to characterize the strength of the halo since lensing can robustly measure this quantity. The Einstein radius is the radius of the ring that is formed by a point source lying perfectly behind a spherical halo with strength R_{h, E}. The normalization ρ_{h, 0} is related to the Einstein radius R_{h, E} by

$$\begin{eqnarray} \rho _{\rm h,0} = A \left\lbrace \begin{array}{@{}ll}[1-\log (2)]^{-1} & {\rm if\ } \tilde{R}_{\rm h,E}=1 \phantom{\frac{\Sigma }{\Sigma }} \\\ms \Big {[} \log \left(\frac{\tilde{R}_{\rm h,E}}{2}\right) + \frac{2}{\sqrt{1-\tilde{R}_{\rm h,E}^2}}& \\\ms \ \ \;\cdot {\rm arctanh} \left(\sqrt{\frac{1-\tilde{R}_{\rm h,E}}{1+\tilde{R}_{\rm h,E}}}\right) \Big {]}^{-1} & {\rm if\ } \tilde{R}_{\rm h,E}<1\\ \Big {[} \log \left(\frac{\tilde{R}_{\rm h,E}}{2}\right) + \frac{2}{\sqrt{\tilde{R}_{\rm h,E}^2-1}} & \\\ms \ \ \;\cdot {\rm arctan} \left(\sqrt{\frac{\tilde{R}_{\rm h,E}-1}{\tilde{R}_{\rm h,E}+1}}\right) \Big {]}^{-1} & {\rm if\ } \tilde{R}_{\rm h,E}>1\\ \end{array}\!{,} \right.\quad\quad \end{eqnarray} \tag{ 10 }$$

where

$$\begin{equation} A=\frac{\Sigma _{\rm crit}}{4 r_{\rm h,0} D_{\rm d}} \tilde{R}_{\rm h,E}^2, \end{equation} \tag{ 11 }$$

and $\tilde{R}_{\rm h,E}=R_{\rm h,E}/r_{\rm h,0}$. Both R_{h, E} and r_{h, 0} are in arcseconds (which can be easily converted to, e.g., kiloparsecs, with the angular diameter distance to the lens, D_d).

For the lensing analysis in Section 8, we also include an external shear component with strength γ_ext and a position angle of ϕ_ext. A shear angle of ϕ_ext = 0 corresponds to an elongation of the images in the east–west direction.

For simplicity, we impose a symmetry condition on the model. To be precise, we require that the centroids of the three components coincide and that the total mass distribution is three-dimensionally axisymmetric. In particular, the NFW profile is either prolate (a = b) or oblate (b = c). Based on the symmetry assumptions and the axis ratio and P.A. of the projected disk (in Table 2), the orientation of the halo is (ϑ, φ) = (50 ± 5°, 145 ± 3°).

Cohn et al. (2001) found that the orientation of the quadrupole moment of their lens model (one-component total mass profile and external shear) agrees with the position angle of the lens galaxy within ∼5°. This implies that the position angle of the projected total mass distribution is aligned with the light distribution of the galaxy. We confirm the alignment by modeling the lens system using a pseudoisothermal elliptic mass distribution (Kassiola & Kovner 1993) in the presence of external shear. The alignment of the projected total mass distribution and the light implies that the NFW halo is oblate (in the inner region probed by lensing) since a prolate halo would lead to a ∼90° difference in the position angles of the projected mass and of the light distributions.

We impose Gaussian priors on (1) the centroid of the total mass distribution, (2) the projected axis ratio and scale radius of the disk, and (3) the orientation of the NFW halo, based on the observed light profile in Table 2. The position angle of the projected disk is set by the orientation of the NFW halo (based on the axisymmetry assumption). We impose uniform priors on the remaining parameters and summarize the priors in Table 4. While all these priors are imposed for the lens modeling, some are not needed for the kinematics, such as the external shear.¹⁷ In total, we have 13 mass parameters: 6 with Gaussian priors, and 7 "free" parameters with uniform priors.

Table 4. Priors on Model Parameters

Description	Parameter	Prior
Centroid in θ1	θ1, c	G(− 0040, 0005)
Centroid in θ2	θ2, c	G(− 0036, 0005)
Disk axis ratio	qd	G(0.63, 0.06)
Disk strength	κd, 0	U(0, ∞)
Disk scale radius	Rd	G(051, 003)
Bulge mass	Mb	U(0, ∞)
Halo flattening	a/c = a/b	U(0.25, 1)
Halo orientation angle	ϑ	G(50°, 5°)
Halo orientation angle	φ	G(145°, 3°)
Halo Einstein radius	Rh, E	U(0, ∞)
Halo scale radius	rh, 0	U(01, 8'') ≡ U(0.75, 60) kpc
External shear strength	γext	U(0, 0.3)
External shear angle	ϕext	U(0, 2π)

Notes. G(μ, σ) is a Gaussian distribution with mean μ and standard deviation σ. U(a, b) is a uniform distribution between boundaries a and b. In cases where b = ∞, the upper boundary is set to a real number that corresponds to masses ≳ 10¹³ M_☉, i.e., larger than galaxy-scale masses.

Download table as:  ASCIITypeset image

For a test mass in the plane of the disk at a radius R from the center, its circular velocity v_c(R) has three contributions:

$$\begin{equation} v_{\rm c}^2 = v_{\rm b}^2 + v_{\rm d}^2 + v_{\rm h}^2, \end{equation} \tag{ 14 }$$

where v_b, v_d, and v_h are the circular velocities of the bulge, disk, and halo mass distributions, respectively.

For the point-mass bulge, the rotational velocity contribution v_b takes on a simple form

$$\begin{equation} v_{\rm b}^2 = \frac{G M_{\rm b}}{R}, \end{equation} \tag{ 15 }$$

where G is the gravitational constant and M_b is the mass of the bulge.

The contribution of the disk also has an analytic expression in terms of Bessel functions I₀, K₀, I₁, and K₁ (e.g., Binney & Tremaine 1987):

$$\begin{equation} v_{\rm d}^2 = 4\pi G\Sigma _{\rm c} R_{\rm d} y^2 (I_0(y)K_0(y) - I_1(y)K_1(y)), \end{equation} \tag{ 16 }$$

where Σ_c is the central surface mass density of the disk (≡ κ_{d, 0}Σ_critq_d), R_d is the scale radius of the disk, y is related to the radial distance R by y = R/2R_e, and R_e is the effective radius of the mass profile of the disk.

For the halo component, we follow Binney & Tremaine (1987) and integrate the oblate spheroidal three-dimensional mass density ρ(m) with m² = x²/q²_h + y² + z² to obtain the circular velocity:

$$\begin{equation} v_{\rm h}^2 = 4\pi G q_{\rm h} \int _0^{R} m^2\frac{\rho (m)}{\sqrt{R^2-m^2\big(1-q_{\rm h}^2\big)}}{\rm d}m, \end{equation} \tag{ 17 }$$

where q_h = a/c is the flattening of the dark matter halo. For the NFW mass density (Equation (8)), we numerically compute the integral in Equation (17).

To model the rotation curve, we follow van der Marel & van Dokkum (2007) to obtain the predicted line-of-sight velocities from the circular velocities which incorporate the effects of the disk inclination, slit, and seeing.

The line-of-sight velocity in the absence of the slit and seeing at a point (x', y') on the sky (with (x', y') = (0, 0) at the center of the galaxy and x' along the major axis of the lens galaxy) is

$$\begin{equation} v_{\rm los}(x^{\prime },y^{\prime }) = \frac{x}{R} \sin (i) v_{\rm c}(R), \end{equation} \tag{ 18 }$$

where $\mbox{\boldmath {${R}$}}=(x, y)$ is the corresponding coordinate in the plane of the disk with x = x', y = y'/cos i, and $R=\vert \mbox{\boldmath {${R}$}} \vert$.

We weight v_los by the modeled galaxy light, and convolve with both a Gaussian PSF with FWHM of 06 and a square top-hat function of size W × P (for slit of width W and pixel scale P) for each (binned) spatial pixel along the slit. The resulting weighted and convolved velocity is the predicted line-of-sight velocity v^pred_los at the corresponding spatial pixel.

The uncertainties of the data points in the rotation curve of Figure 3 are only statistical. We include a systematic uncertainty of 40 km s⁻¹ (which is comparable to the statistical uncertainty) to account for deviations from our assumptions of the smooth axisymmetric mass model, slight offset (if any) of the galaxy centroid in the rotation curve, and the unknown systemic velocity. We explore a range of values for the systemic velocity and choose the value that optimizes the rotation curve fitting. We add the statistical and the systematic uncertainty in quadrature to obtain the final uncertainty on the line-of-sight velocity. The amount of systematic uncertainty is set such that the rotation curve can be modeled with a reduced χ² ∼ 1. The inclusion of the systematic uncertainty is crucial for not underestimating the uncertainty on the resulting mass model parameters.

The likelihood for the rotation curve data is

$$\begin{equation} P(\mbox{\boldmath {${d}$}}_{\rm D} | \mbox{\boldmath {${\eta }$}}) = \frac{1}{Z_{\rm D}} \exp \left[ -\frac{1}{2} \sum _{i=1}^{N_{\rm D}} \frac{\big(v_{{\rm los},i}^{\rm obs} - v_{{\rm los},i}^{\rm pred} \big)^2}{\sigma _i^2} \right], \end{equation} \tag{ 19 }$$

where N_D is the number of data points (=15), v^obs_{los, i} is the observed line-of-sight velocity of data point i, v^pred_{los, i} is the predicted line-of-sight velocity from our model and observational setup, σ_i is the uncertainty in the velocity, and Z_D is the normalization given by

$$\begin{equation} Z_{\rm D} = (2\pi)^{N_{\rm D}/2} \prod _{i=1}^{N_{\rm D}} \sigma _i. \end{equation} \tag{ 20 }$$

We use MultiNest (Feroz & Hobson 2008; Feroz et al. 2009) to sample the posterior PDF of the 13 model parameters. Figure 7 shows the resulting constraints on four of the parameters: halo flattening a/c, halo Einstein radius R_{h, E}, halo scale radius r_{h, 0}, and disk strength κ_{d, 0}. The disk–halo degeneracy appears as the anti-correlation between the disk strength (κ_{d, 0}) and the halo Einstein radius (R_{h, E}): the more massive the disk, the less massive the halo, and vice versa. The halo scale radius is poorly constrained. Furthermore, the kinematics data provide little information on the halo shape.

Zoom In Zoom Out Reset image size

Figure 8 shows the predicted rotation curve of the most probable mass distribution with a reduced χ² = 1.0 based on the kinematics data. The bulge contributes very little to the mass of the system, consistent with the small effective radius observed for the bulge light distribution. At small (large) radii, the disk (dark matter halo) dominates in the enclosed mass.

Zoom In Zoom Out Reset image size

We briefly describe how to obtain the deflection angles for the three mass components in our model. The deflection angles are used to solve for the predicted radio image positions that are needed for constructing the likelihood of the lensing data.

8.1.1. Disk

The lensing deflection angles for an elliptical mass component following an exponential disk profile has no analytical form. Hence, to model the stellar disk component, we employ a "chameleon profile" mimicking the mass profile of an exponential disk whose deflection angles are analytical. This profile is inspired by the chameleon profile used in Maller et al. (2000). We use the following approximation for the dimensionless surface mass density of the disk:

$$\begin{eqnarray} \hspace{-3pc}\kappa _{\rm d}(R^{\prime }) & = & \kappa _{\rm d,0} \exp [-R^{\prime }/R_{\rm d}] \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} & \approx & \kappa _{\rm d,0} s_1 \left[ \frac{1}{\sqrt{R^{\prime 2}+s_2^2}} - \frac{1}{\sqrt{R^{\prime 2}+s_3^2}}\right], \end{eqnarray} \tag{ 22 }$$

where s₁ = 4.6849R_d, s₂ = 1.1720R_d, and s₃ = 1.4518R_d. The projected enclosed mass within radius $\tilde{R^{\prime }}$ for the chameleon profile agrees with that of the exponential profile within 2% for the range of radii spanned by the lensed images ($\tilde{R^{\prime }}=0.6R_{\rm d}$–1.6R_d). Each term in the square brackets in Equation (22) is in the form of a pseudoisothermal elliptic mass distribution whose deflection angles can be computed analytically (Kassiola & Kovner 1993).

8.1.2. Bulge

The bulge is modeled as a point mass, and the deflection angle has a simple closed form (e.g., Schneider et al. 2006)

$$\begin{equation} \alpha (\xi) = \frac{4G}{c_{\rm l}^2}\frac{M_{\rm b}}{\xi }, \end{equation} \tag{ 23 }$$

where G is the gravitational constant, c_l is the speed of light, M_b is the mass of the bulge, and ξ is the impact parameter (i.e., the distance between the light ray and the point mass on the lens plane).

8.1.3. Dark Matter Halo

We use the radio image positions listed in Table 1 to constrain the mass model parameters. We add a positional uncertainty of 10 mas in quadrature to the uncertainties in the table to account for systematic effects such as the residual features in the lens galaxy light fit (Section 4.2), presence of substructure (e.g., Chen et al. 2007), and scatter-broadening through the disk of the lens (e.g., Marlow et al. 1999; Biggs et al. 2004). The lensing arcs in the NIRC2 images could in principle be used in addition to the radio image positions to constrain the lens mass distribution; however, these arcs (especially the one associated with components 3 and 6) in practice suffer from dust extinction. Only the F160W image of the HST data has sufficient S/N to be used for dust correction, and Suyu et al. (2009) showed that dust correction based on two bands are prone to systematic effects. Therefore, we do not include the dust-affected NIRC2 arcs for constraining the mass distribution.

The likelihood of the radio image positions is

$$\begin{equation} P(\mbox{\boldmath {${d}$}}_{\rm L} |\mbox{\boldmath {${\eta }$}}) = \frac{1}{Z_{\rm {L}}} \exp {\left[-\frac{1}{2}\displaystyle \sum _{j=1}^{N_{\rm sys}}\displaystyle \sum _{i=1}^{N_{\rm im}^j} \frac{\big\vert \mbox{\boldmath {${\theta }$}}_{i,j}^{\rm obs}-\mbox{\boldmath {${\theta }$}}_{i,j}^{\rm pred}(\mbox{\boldmath {${\eta }$}})\big\vert ^2}{\sigma _{i,j}^2} \right]}, \end{equation} \tag{ 24 }$$

where N_sys is the number of multiply imaged systems, N^j_im is the number of multiple images in system j, $\mbox{\boldmath {${\theta }$}}_{i,j}^{\rm obs}$ is the observed image position, $\mbox{\boldmath {${\theta }$}}_{i,j}^{\rm pred}(\mbox{\boldmath {${\eta }$}})$ is the modeled image position, σ_{i, j} is the uncertainty in the observed image position, and Z_L is the normalization given by

$$\begin{equation} Z_{\rm L}={(2\pi)^{N_{\rm L}/2} \displaystyle \prod _{j=1}^{N_{\rm sys}} \displaystyle \prod _{i=1}^{N_{\rm im}^j} \sigma _{i,j}} \end{equation} \tag{ 25 }$$

with

$$\begin{equation} N_{\rm L}=\displaystyle \sum _{j=1}^{N_{\rm sys}} {N_{\rm im}^j}. \end{equation} \tag{ 26 }$$

The source position for each system of multiple images is needed to predict the image positions. We model the source position as the weighted average of the mapped source positions from the observed image positions and the deflection angles from the mass distribution. Specifically, for each image system, we take the average of the mapped source position $\mbox{\boldmath {${\beta _i}$}}$ weighted by $\sqrt{\mu _{i}}/\sigma _{i}$, where μ_i and σ_i are, respectively, the modeled magnification and the positional uncertainty of image i.¹⁸ In doing so, we have in effect marginalized the source position parameters by approximating the lensing likelihood as having a delta function at the weighted source position for each image system.

We model the three-component mass distribution of the spiral galaxy based on the lensing data with Glee ("Gravitational Lens Efficient Explorer"), a software developed by A. Halkola & S. H. Suyu.¹⁹ For a set of mass model parameter values, Glee computes the values of the lensing likelihood in Equation (24) and the prior in Equation (13); these are the two ingredients needed to obtain the posterior PDF in Equation (12). As in the kinematics analysis, we use MultiNest (Feroz & Hobson 2008; Feroz et al. 2009) to sample the posterior PDF.

Figure 9 shows the resulting constraints on the same parameters as those in Figure 7. The degeneracy between the disk strength and the Einstein radius of the dark matter halo is visible, though not as strong as in the case of the kinematics-only analysis. The halo is highly flattened (a/c ∼ 0.3). The flattening is degenerate with the disk strength as shown in the top left panel: a flattened halo is less massive and requires a more massive disk. The flattening is also degenerate with the halo Einstein radius: massive halos with high R_{h, E} need to be more flattened to reproduce the overall ellipticity of the projected mass as constrained by the lensing data. The scale radius of the dark matter halo is not constrained, as expected since lensing only probes the distribution in the radial range spanned by the images, i.e., ∼05, or ∼4 kpc. Nonetheless, a small value of ≲ 1'' is rejected by the data at 95% CI.

Zoom In Zoom Out Reset image size

The most probable lensing model (with highest posterior PDF) has a reduced χ² = 0.9. We show in Figure 10 the critical curves (solid) and caustics (dashed) of the most probable lensing model. The open symbols are the observed image positions, and the solid symbols are the modeled source positions. The figure illustrates the configuration of the 10 image system in relation to its three-component source (the source positions are clustered into three groups). The first group of sources is outside the astroid caustics and produces components 1a and 8. The second group is inside the astroid caustics and produces components 1, 3, 4, and 6. The third group is near a fold of the caustics and produces components 2, 5, and 7.

Zoom In Zoom Out Reset image size

Although the kinematics constraints are significantly weaker than the lensing constraints, the κ_{d, 0}–R_{h, E} panel (top middle) in Figure 11 shows that the kinematics and lensing contours are tilted at different angles. Thus, the combination of the two in principle breaks the disk–halo degeneracy. In the case of B1933+503, most of the constraints on the mass distribution come from the lensing data: the 10 radio images, which both span a large range of radii and have positional uncertainties of only a few mas, provide a strong leverage in discerning contributions from the disk, bulge, and halo. This is in contrast to typical lenses with only two or four images which result in stronger disk–halo degeneracies (e.g., Trott et al. 2010). On the other hand, the S/N of the rotation curve data for B1933+503 is only modest given the high lens redshift of z_l = 0.755. Nonetheless, the kinematics data in B1933+503 are informative in excluding the high disk masses allowed by lensing, which has interesting consequences that are discussed in Section 9.4.

A disk is considered to be maximal if the fractional contribution to the circular velocity of the disk at 2.2 R_d is v_d(2.2 R_d)/v_tot(2.2 R_d) = 0.85 ± 0.1 (Sackett 1997). Based on our lensing and kinematics model, B1933+503 has v_d(2.2 R_d) = 248⁺²²₋₂₆ km s⁻¹ and v_tot(2.2 R_d) = 326 ± 8 km s⁻¹. In comparison to the Milky Way's circular velocity at the position of the sun, v_c(R₀) = (219 ± 20) R₀/(8 kpc) km s⁻¹ (Reid et al. 1999), the spiral galaxy in B1933+503 is significantly more massive. The resulting disk contribution in B1933+503 of v_d(2.2 R_d)/v_tot(2.2 R_d) = 0.76^+0.05_−0.06 suggests that the disk is marginally submaximal.

In the first column of Figure 11, we see that the dark matter halo in B1933+503 is oblate with a/c = 0.33^+0.07_−0.05 in order to fit to the lensing observations. The axis ratio of the projected surface mass density of the halo is 0.67 ± 0.04, which is slightly rounder than the axis ratio of projected disk of 0.53 ± 0.03. It appears that only projected mass distributions of the lens with axis ratios of ∼0.6 are consistent with the lensing data. If we consider the subset of the MultiNest sample with a/c ∼ 0.5 (which corresponds to a projected axis ratio of ∼0.75), then the lensing reduced χ² increases from 0.8 to 2.4 due to misfits in systems II, V, VI, and VII in Table 1. The global VLBI data with better positional accuracies and, in particular, the four-image systems II and VI enforce the high ellipticity in the projected surface mass density and hence the highly oblate halo. The high ellipticity is robust against assumptions on the dark matter halo profile. In particular, when we replace the NFW halo with either a singular or cored isothermal profile, the radio data still require a high ellipticity for the halo.

N-body simulations of dark matter halos indicate that the halos are typically triaxial and the axis ratio between the short and long axis (i.e., a/c in our notation) is ∼0.3–0.9 (95% CI) for halo masses of 10¹¹–10¹² M_☉ (e.g., Bett et al. 2007; Macciò et al. 2008). Simulations with baryons find that the central galaxy tends to make the triaxial halo essentially oblate and less flattened (e.g., Abadi et al. 2010). Our study of B1933+503 with the resulting a/c ∼ 0.3 for the oblate NFW halo suggests that baryons are effective at making the halo oblate in the inner regions of the galaxy.

Since the kinematics and lensing data probe the inner ∼1'' = 7.5 kpc of the galaxy, it is not surprising that the scale radius r_{h, 0} of the dark matter halo in our model (where the density transitions from ρ ∼ r⁻¹ to r⁻³) is not constrained but has a lower limit of 21 = 16 kpc (95% CI). The virial radius is difficult to extrapolate from the model given the uncertain scale radius. For galaxies with total mass of ∼10¹² M_☉ (applicable to B1933+503), typical virial radii are ∼250 kpc. This implies a concentration of ≲ 16 that is consistent with concentration–mass relations for galaxy-scale halos (e.g., Macciò et al. 2008; Duffy et al. 2010). To better constrain r_{h, 0}, a rotation curve that extends to larger radii is needed.

From our three-component model of the galaxy, we can determine the fraction of dark matter as a function of radius by integrating the mass enclosed within spherical radii. Since we employ a parameterized model and the parameters of the model are constrained by the lensing and kinematics data, the enclosed mass can be computed for all radii. At small and large radii where there is no (or low S/N) data, the enclosed mass is effectively extrapolated based on the form of the parametrized model and the measured values of the mass parameters. The left panel in Figure 12 shows the mass enclosed for the disk and the dark matter halo. The bulge, which is unresolved and modeled as a point, has mass M_b = (1.6 ± 0.3) × 10¹⁰ M_☉ for all radii. In the inner 10 kpc of the galaxy, we see that the disk dominates in mass, and beyond that, the dark matter halo dominates.

Zoom In Zoom Out Reset image size

The right panel shows the dark matter mass fraction within a sphere of radius r. The rise in the dark matter fraction is similar to the analysis of the spiral lens SDSS J2141−0001 by Dutton et al. (2011). The dark matter mass fraction within 2.2 disk scale lengths for B1933+503 is f_{DM, 2.2} = 0.43^+0.10_−0.09. Within the effective radius, the mass fraction is f_{DM, e} = 0.37^+0.09_−0.08, which is consistent with the ranges of values found in previous lensing and kinematics analyses of early-type galaxies (e.g., Auger et al. 2010a; Barnabè et al. 2011).

The shapes of the low- and high-mass end of the stellar IMF are difficult to determine observationally since low-mass stars are intrinsically faint and high-mass stars are low in abundance. In this section, we describe the estimation of the stellar mass in B1933+503 based on SPS models for Chabrier (2003) and Salpeter (1955) IMFs. We compare these masses to the independently measured disk mass from lensing and kinematics to learn about the IMF.

Using the broadband photometries in Table 3, we follow Auger et al. (2009) to infer the stellar mass for the disk based on composite SPS models from Bruzual & Charlot (2003). These models have been employed in many studies including Sloan Digital Sky Survey (SDSS) galaxies (e.g., Kauffmann et al. 2003; Gallazzi et al. 2005) and spiral lens analyses similar to ours (e.g., Dutton et al. 2011). Dust extinction (intrinsic) is taken into account in these models, and our multiband photometries allow a good handle on dust with only a small broadening of uncertainties (e.g., Bell & de Jong 2001). Our inference of the stellar mass provides uncertainties that incorporate parameter degeneracies in the stellar population models (e.g., between age and metallicity). In Figure 13, we plot the inferred stellar mass in red-dashed (blue-dotted) assuming a Chabrier (Salpeter²⁰) IMF. These two IMFs cover the range applicable to spiral galaxies (Bell & de Jong 2001).

Zoom In Zoom Out Reset image size

The lensing and kinematics analysis provides an independent measurement of the total disk mass of log₁₀(M_disk/M_☉) = 11.17^+0.08_−0.10 that includes stars and gas (thin solid curve in Figure 13). To extract the stellar mass contribution, we assume that the cold gas accounts for 20% ± 10% of the total mass (e.g., Dutton & van den Bosch 2009) and is distributed like the stars. This provides an upper limit on the mass contribution of the gas to the disk (in comparison to scenarios where the gas is more extended than the stars). For each sample in the posterior PDF of the disk mass, we draw a random number f from a Gaussian distribution centered on 0.2 with a standard deviation of 0.1 that is truncated between [0,1], and subtract the fraction f from the total disk mass. The gas-subtracted disk mass from the lensing and kinematics model is log₁₀(M_*/M_☉) = 11.06^+0.09_−0.11, and the distribution is shown by the thick solid curve in Figure 13.

Comparing the thick solid curve to the dashed and dotted curves in the figure, we see that our mass model of B1933+503 favors a Chabrier-like IMF. We can quantitatively rank the two IMFs by computing the Bayesian evidence which is the integral under the product of the PDF from lensing and kinematics (thick solid) and the PDF from the SPS (dashed or dotted). The ratio of the Bayesian evidence for Chabrier to Salpeter is 7.2; in other words, the probability of the IMF being Chabrier is 7.2 times higher than the probability of the IMF being Salpeter for B1933+503.

Using the Chabrier IMF, we obtain the rest-frame stellar mass-to-light ratio of the disk to be M_*/L_V = (0.6 ± 0.3)M_☉/L_{V, ☉}. By passively evolving to z = 0, we obtain the present-day stellar mass-to-light ratio of M_*/L_V = (1.7^+1.3_−0.9)M_☉/L_{V, ☉}, in agreement with typical values found for local galaxies (e.g., Trott et al. 2010; van de Ven et al. 2010; Courteau & Rix 1999; van der Kruit & Freeman 2011, and references therein).

Our finding of a preference toward a Chabrier-like IMF for the spiral galaxy is in agreement with the results in Dutton et al. (2011) and Ferreras et al. (2010). Nonetheless, studies of massive elliptical galaxies favor Salpeter-like IMFs (e.g., Treu et al. 2010; Auger et al. 2010b; van Dokkum & Conroy 2010; Spiniello et al. 2011). This supports a non-universal IMF for galaxies that is dependent on mass and/or Hubble type.