# The Einstein Cross: constraint on dark matter from stellar dynamics and gravitational lensing

## Abstract

We present two-dimensional line-of-sight stellar kinematics of the lens galaxy in the Einstein Cross, obtained with the GEMINI 8m telescope, using the GMOS integral-field spectrograph. The stellar kinematics extent to a radius of ″ (with ″ spaxels), covering about two-thirds of the effective (or half-light) radius ″ of this early-type spiral galaxy at redshift , of which the bulge is lensing a background quasar at redshift . The velocity map shows regular rotation up to km s around the minor axis of the bulge, consistent with axisymmetry. The velocity dispersion map shows a weak gradient increasing towards a central (″) value of km s. We deproject the observed surface brightness from HST imaging to obtain a realistic luminosity density of the lens galaxy, which in turn is used to build axisymmetric dynamical models that fit the observed kinematic maps. We also construct a gravitational lens model that accurately fits the positions and relative fluxes of the four quasar images. We combine these independent constraints from stellar dynamics and gravitational lensing to study the total mass distribution in the inner parts of the lens galaxy.

We find that the resulting luminous and total mass distribution are
nearly identical around the Einstein radius ″,
with a slope that is close to isothermal, but which becomes
shallower towards the center if indeed mass follows light. The
dynamical model fits to the observed kinematic maps result in a
total mass-to-light ratio (in the
-band). This is consistent with the Einstein mass M divided by the (projected) luminosity within
, which yields a total mass-to-light ratio of , with an error of at most a few per cent. We estimate
from stellar populations model fits to colors of the lens galaxy a
stellar mass-to-light ratio from to .
Although a *constant* dark matter fraction of 20 per cent is
not excluded, dark matter may play no significant role in the bulge
of this L early-type spiral galaxy.

###### Subject headings:

gravitational lensing — stellar dynamics — galaxies: photometry — galaxies: kinematics and dynamics — galaxies: structure## 1. Introduction

In the cold dark matter (CDM) paradigm for galaxy formation (e.g. Kauffmann & van den Bosch, 2002), galaxies are embedded in extended dark matter distributions with a specific and (nearly) universal radial profile (e.g. Navarro et al., 1997; Merritt et al., 2005). Measurements of rotation curves from neutral hydrogen (H i) observations in the outer parts of late-type galaxies have provided evidence for the presence of dark matter in these systems already more than two decades ago (e.g. van Albada et al., 1985). In the outer parts of early-type galaxies, however, cold gas is scarce (but see e.g. Franx et al., 1994; Morganti et al., 1997; Weijmans et al., 2008), so evidence for a dark matter halo has to come from other tracers, such as kinematics of stars, planetary nebulae and globular clusters (e.g. Carollo et al., 1995; Gerhard et al., 2001; Romanowsky et al., 2003; Côté et al., 2003) or hot X-ray gas (e.g. Fukazawa et al., 2006; Humphrey & Buote, 2010). However, these tracers are not always (sufficiently) available, the observations are often challenging, and the interpretation modeling dependent. Also in the inner parts of galaxies the amount, shape and profile of dark matter is still poorly known (e.g. Primack, 2004), even though stellar kinematics are readily available.

A fundamental problem in using stellar kinematics (as well as other collisionless kinematic tracers) for this purpose is the mass-anisotropy degeneracy: a change in the measured line-of-sight velocity dispersion can be due to a change in total mass, but also due to a change in velocity anisotropy. Both effects can be disentangled by measuring also the higher-order velocity moments (Dejonghe, 1987; van der Marel & Franx, 1993; Gerhard, 1993), but only the inner parts of nearby galaxies are bright enough to obtain the required high-quality kinematic measurements (e.g. van der Marel, 1991; Gerhard et al., 2001; Cappellari et al., 2006). A unique alternative to break the mass-anisotropy degeneracy is provided by gravitational lensing. In case of strong lensing, the mass of a foreground galaxy bends the light of a distant bright object behind it, resulting in multiple images. From the separation and fluxes of the images the total mass distribution of the lens galaxy can be inferred directly. The velocity anisotropy then can be determined from the observed velocity dispersion, without the need for higher-order velocity moments.

Treu & Koopmans (2004, and references therein) have applied this approach to several strong lensing systems, of which 0047-281 (Koopmans & Treu, 2003) is the best constrained case, with three radially separated velocity dispersion measurements extending to about the effective radius of the lens galaxy. They measure the mass within the Einstein radius by fitting a singular isothermal ellipsoid to the positions of the quasar images. This Einstein mass is used to set the amplitude of the total (stellar and dark) matter distribution, which they assume to be spherical. The constant stellar mass-to-light ratio determines the contribution of the stars with a fixed radial profile, with the remainder due to dark matter with a single power-law profile with slope . They then compare the dispersion profile predicted by the spherical Jeans equations, for an ad-hoc assumption of the velocity anisotropy , with the observed dispersion measurements. Based on a reasonably constrained and an upper limit on , they conclude that a significant amount of dark matter is present in the inner parts of the lens galaxy, with a slope shallower than the nearly isothermal total mass distribution. An additional (external) constraint on (or ) is needed to go beyond this limit on the dark matter distribution. Even so, their results are limited by too few kinematic constraints (which leaves the anisotropy degenerate), and by the use of a simple spherical dynamical model.

Clearly, most lens galaxies are significantly flattened and so cannot be well-described by spherical models. Non-spherical models provide a more realistic description of the lens galaxy, but the increase in freedom requires also (significantly) more spatially resolved kinematic measurements to constrain them. Only a few of the known strong gravitational lens systems are close enough to obtain such kinematic data. Even then, one can make (ad-hoc) assumptions on the velocity distribution, such as velocity isotropy in the meridional plane of an axisymmetric model. The latter restriction to a distribution function of two (instead of three) integrals of motion allows for the recovery of the intrinsic shape and mass distribution in the presence of two-dimensional kinematic data (Barnabè & Koopmans, 2007; Czoske et al., 2008; Barnabè et al., 2009), possibly even without the additional constraint provided by gravitational lensing (e.g. Cappellari et al., 2006).

In this paper, we relax the two-integral assumption and find that one can still derive independent determinations of the intrinsic mass distribution from both gravitational lensing and stellar dynamics with two-dimensional kinematic data. We have observed the gravitational lens system QSO 2237+0305, well-known as the Einstein Cross, with the integral-field spectrograph GMOS on the GEMINI-North Telescope. We fit the extracted stellar velocity and velocity dispersion maps of the inner parts of the lens galaxy at a redshift with axisymmetric dynamical models. The resulting intrinsic mass distribution agrees well with that inferred from the lens model that fits the four quasar image positions and relative fluxes. Furthermore, we compare this total mass distribution with the stellar mass distribution, based on an independent estimate of the stellar mass-to-light ratio , to place constraints on the dark matter distribution in the inner parts of the lens galaxy.

In Section 2, we briefly describe the Einstein Cross, and we present the photometric and spectroscopic observations. In Section 3, we extract the two-dimensional stellar kinematics, describe the dynamical modeling method, and construct lens models. From the latter we infer the total mass distribution, which we compare in Section 4 with the luminous mass distribution inferred from the observed surface brightness. We then build axisymmetric dynamical models, and compare the resulting dynamical mass-to-light ratio with that inferred from the lens models, and in turn with an estimate of the stellar mass-to-light ratio. In Section 5 we discuss our findings and we summarize our conclusions.

Throughout we adopt the WMAP cosmological parameters for the Hubble constant, the matter density and the cosmological constant, of respectively km s Mpc, and (WMAP3; Spergel et al., 2007).

## 2. Observations and data reduction

### 2.1. The Einstein Cross

The Einstein Cross is the well-known gravitational lens system QSO 2237+0305 (RA: 22h 40m 30.3s, Dec: 21′ 31″, J2000). In this system, a distant quasar at redshift is lensed by the bulge of the early-type spiral PGC 069457 at (angular diameter distance Mpc, kpc), resulting in a cross of four bright images separated by about 1.8″.

The Einstein Cross has long been the closest strong gravitational lens system known, and has been very well studied since its discovery by Huchra et al. (1985). There is a wealth of ground- and space-based imaging data at all wavelengths (e.g. Falco et al., 1996; Blanton et al., 1998; Agol et al., 2000; Dai et al., 2003). The resulting precise measurements of the positions and relative fluxes of the quasar images can be used to construct a detailed lens model.

In contrast, kinematic data of the lens galaxy is very scarce, with only one measured central stellar velocity dispersion (Foltz et al., 1992) and two H i rotation curve measurements in the very outer parts (Barnes et al., 1999). There are several previous integral-field studies of the Einstein Cross: TIGER: Fitte & Adam (1994); INTEGRAL: Mediavilla et al. (1998); CIRPASS: Metcalf et al. (2004). However, none of these studies were concerned with the stellar kinematics of the lens galaxy, but instead investigated the quasar spectra.

### 2.2. Imaging

*Left panel*: the contours of the WFPC2/F555W -band image reveal clearly the bulge, spiral arms and bar embedded in the large-scale disk of this early-type spiral galaxy. The ellipticity measured from the MGE fit (contours) is used to estimate the inclination.

*Middle panel*: the central 8″8″ of the WFPC2/F814W -band image, of which the MGE fit is used to construct the (stellar) luminosity density model of the lens galaxy. We use the -band image instead of the longer exposed -band image as it tracers better the old stellar population and is less sensitive to extinction and reddening. The four quasar images are masked out during the MGE fit. In both WFPC2 images the contours are in steps of mag/arcsec. The images are rotated such that North is up and East is to the left.

*Right panel*: The scale-free lens model with slope (dashed contours) fits the positions and relative fluxes of the quasar images, indicated by the filled circles. Superposed is the MGE fit (solid contours).

We use two WFPC2 images retrieved from the HST-archive (F555W -band image, 1600 seconds, PI: Westphal, and F814W -band image, 120 seconds, PI: Kochanek; see left and middle panels of Fig. 1) to determine the luminosity density of the lens galaxy. We correct the -band image for a Galactic extinction of (Schlegel et al., 1998), and we convert to solar units using the WFPC2 calibration of Dolphin (2000), while assuming an absolute -band magnitude for the Sun of 4.08 mag (Table 2 of Binney & Merrifield, 1998). From a de Vaucouleurs profile fit to the -band photometry in the inner bulge-dominated region, we obtain an effective radius ″, which is consistent with previous measurements (e.g. Racine, 1991).

For the construction of the lens model, we use the accurate positions
of the quasar images from the website of the CASTLES survey^{10}

### 2.3. Integral-field spectroscopy

Observations of the Einstein Cross lens system were carried out using the integral-field unit of the GMOS-North spectrograph (Murray et al., 2003; Hook et al., 2004) on July 17 and August 1 2005 as part of the program GN-2005A-DD-7; the seeing was about ″. The data were obtained using the IFU two-slit mode that provides a field-of-view of 5″7″. An array of 1500 hexagonal lenslets, of which 500 are located 1′ away from the main field to be used for sky subtraction, sets the 02 spatial sampling. Eight individual science exposures of 1895 seconds each were obtained during the two nights, resulting in a total on-source integration time of hours. An offset of 03 was introduced between exposures to avoid bad CCD regions or lost fibers. The R400-G5305 grating in combination with the CaT-G0309 filter was used to cover a wavelength range of 7800–9200 Å with a FWHM spectral resolution of 2.8 Å.

For the data reduction we used an updated version of the officially
distributed Gemini IRAF^{11}*L.A. Cosmic* algorithm
by van Dokkum, 2001) to each science exposure, and wavelength
calibration after the extraction of the data. A careful flat-fielding
procedure, crucial for proper removal of fringes in the spectral
direction, was carried out using the Quartz halogen (QH) lamp.
Wavelength calibration was performed with CuAr lamp exposures taken
before each science exposure. At the observed wavelengths, a complex
spectrum of HO absorption features overlaps with the CaT lines we
are interested in. We constructed a correction spectrum from
observations of the white dwarf star Wolf1346, taken with the same
instrumental setup as our science frames. We checked the range of
fiber-to-fiber variations of the spectral resolution of the instrument
by measuring the width of the sky lines in each fiber, resulting in
the nominal FWHM value of Å for most of the data
cubes. In order to combine all the data, we homogenized the science
frames by convolving all spectra to an instrumental FWHM resolution of
3.1 Å (or km s) and
resampling the spectra to the same range and sampling in wavelength.
After interpolating each science frame to a common spatial grid,
taking into account the small spatial offsets applied during the
observations, we sum the spectra sharing the same position in the sky
to produce the final merged data cube.

The significant contribution from the sky lines to the overall spectrum of the lens galaxy made the sky subtraction the most challenging step in the data reduction process. Although the scatter in the instrumental resolution was small within the individual science frames, the combined effect of other steps in the data reduction (e.g. fringing, telluric absorption, resampling of data in wavelength), introduced alterations in the shape of the sky lines (from the sky fibers) that were not equally reproduced in the science fibers. As a consequence any attempt to use the sky from the sky fibers resulted in serious systematic effects (i.e. P-Cygni like residuals) in the galaxy spectra underneath. In order to minimize these effects, we chose to extract the sky spectra from the outer most regions of our fully merged science data cube, where no galaxy light was appreciable. The sky was then subtracted optimally as described below.

## 3. Analysis

### 3.1. Two-dimensional stellar kinematics

To measure reliable stellar kinematics we first co-add spectra using
the adaptive spatial two dimensional binning scheme of
Cappellari & Copin (2003) to obtain in each resulting (Voronoi) bin a
minimum signal-to-noise ratio (S/N). The resulting spectra are
spectrally rebinned to steps of constant , equivalent to a
sampling of 24 km s per pixel. We adopt the single stellar
population (SSP) models of Vazdekis et al. (2003) as spectral
templates, which are convolved to have the same effective instrumental
broadening as the GMOS data. Next, using an updated version (v4.2)
of the penalized pixel-fitting (pPXF) algorithm of
Cappellari & Emsellem (2004), a non-negative linear combination of these
templates is convolved with a Gaussian line-of-sight velocity
distribution (LOSVD) and fit to each observed spectrum, to derive the
mean line-of-sight velocity and velocity dispersion per
bin on the plane of the sky. This version of the code allows for the
inclusion of a sky spectrum into the set of stellar templates, and
determines the optimal scaling of the sky to the overall spectrum
(see also Weijmans et al., 2009). We used this procedure to
*clean* the individual spectra in our merged data cube, before we
derived the final stellar kinematics.

As an example of how this routine performed, we show in the top panel of Fig. 2 the observed spectrum from the center of the lens galaxy (thick black curve), and the cleaned galaxy spectrum (thin black curve) fitted by a composite of stellar population models (blue curve). The sky emission-lines can severely affect the Ca II triplet absorption lines that we use to extract the stellar kinematics, especially at larger radii where the the relative contribution from the sky is more significant. In particular, in our experiments, a poor sky subtraction caused significant systematic effects in the map of at the receding side of the galaxy, where one of the absorption lines is redshifted onto (the residual of) a sky emission line. The map of , however, was (nearly) unaffected, showing regular rotation with a straight zero-velocity curve that coincides with the minor axis of the (bulge) photometry. These tests already show that, at least in the inner parts of the lens galaxy, both the photometry and kinematics are consistent with axisymmetry.

If we now *assume* this symmetry from the start and apply it
during the extracting of the stellar kinematics, we can not only
increase the S/N, but at the same time also suppress or even remove
systematic effects. To do this, we make use of a built-in
functionality in the pPXF algorithm that allows diametrically opposed
spectra to be fitted simultaneously. The symmetry imposed by this
technique can be either mirror-symmetric or point-symmetric,
corresponding to the behavior of the LOSVD in axisymmetric or triaxial
stellar systems, respectively. In practice, this is done following the
method outlined by Rix & White (1992), whereby the two spectra
from opposed positions in the galaxy are laid end-to-end to create a
single vector. The template library is also doubled in length, where
the second half is convolved by the symmetric counterpart of the
broadening function used for the first half (simply reflecting the
broadening kernel around the systemic velocity). The two galaxy
spectra are therefore fitted simultaneously, with the implied
assumption of symmetry.

Performing a symmetric extraction of the kinematics presents a number of advantages, and has been used by numerous authors not just as a quality-check (e.g. Rix & White, 1992), but also when the assumption of symmetry is an integral part of the analysis (e.g. Gebhardt et al., 2003). The two main advantages for the case presented here are in reducing the impact of systematic effects in the data (including sky subtraction) and increasing the effective spatial resolution. The latter is possible since under the assumption of axisymmetry, all moments of the LOSVD are symmetric about the galaxy major axis. We can therefore ‘fold’ the data cube about the major axis, combining spectra from the two hemispheres, which increases the effective S/N within a given spatial element, and therefore reduces the required spatial bin size.

We interpolate the merged GMOS data cube onto a new spatial grid with and coordinates parallel to the major and minor photometric axes respectively. We then fold the data-cube about the major axis by summing together spatial pixels (‘spaxels’) with equal positions. Before spatially binning the spectra to a minimum S/N, the folded data cube is folded again, into a single quadrant by combining spaxels with equal positions. The Voronoi bins are derived in this quadrant, and the resulting bin centers are mirrored around . The spectra are combined according to these Voronoi bins to ensure symmetry around the minor axis. At the field perimeter, it is generally not the case that both sides of the galaxy are sampled at the exact same positions, due to the alignment of the GMOS field with respect to the galaxy’s principle axes. This results in some bins where the S/N of the opposed spectra is far from equal, or in some cases where observations are only present on one side of the galaxy. The former is taken into account by consideration of the error spectrum of each bin. For the latter, the bin is fitted as a single spectrum, giving kinematic values only at that position, without a mirrored counterpart. Fig. 3 shows how the (axi)symmetric kinematic extraction performed at three different position along the major axis.

From Monte-Carlo simulations for data of this spectral range and resolution, a minimum S/N of 20 is deemed optimal. This results in an average error in both and of about km s, while still preserving the spatial resolution of the data. However, towards the edge of the GMOS field, the errors increase to km s in and km s in . The final symmetrically binned data cube consists of bins, with the largest containing individual 0.2″ spaxels. The resulting (original and symmetrized) maps of and are presented in Fig. 4, with superposed contours of the integrated flux from the data cube. The velocity map shows regular rotation with amplitude of up to km s. The dispersion map shows somewhat higher values of at the position of the quasar images. This is caused by imperfect subtraction of the quasar continuum, which acts to strongly dilute the absorption lines and make them appear broader.

In our stellar template fit to the observed spectra, we include additive polynomials up to 10th order to account for remaining systematic effects in the data, e.g. due to imperfect flat fielding. At the same time, we use these additive components to mimic the contaminating contribution of the quasar images, which impart a strong featureless continuum on top of the galaxy’s integrated stellar spectrum. This method nicely removes nearly all contribution from the quasar, except where the quasar images are the brightest, in which case cannot be reliable measured. The middle panel of Fig. 3 provides an example of such a case in which the spectra come from the opposite positions 1.0″ from the center along the major axis, right in the middle of two quasar images. Even though the composite stellar population model fits to the Ca II triplet still look reasonable, we exclude these regions when we fit dynamical models to the and maps. However, between the quasar images, the spectra are dominated by the stellar light of the galaxy, giving a reliable central (″) dispersion km s. Likewise, at the edges of the field, the quasar contamination is small, providing a reliable gradient in the velocity dispersion.

Comparing in Fig. 4 the symmetrized stellar kinematics with what is obtained without imposing axisymmetry, there is a marked reduction of noise and systematic effects, especially in the velocity dispersion map. However, the main features of regular rotation, decreasing at larger radii, as well as corresponding velocity and dispersion amplitudes, are present in both approaches, albeit with larger uncertainty and systematics in the non-symmetrized case. Given the various advantages, we consider hereafter the symmetrically extracted kinematics.

### 3.2. Dynamical models

As mentioned above, both the photometry and (unsymmetrized) kinematics
in the inner parts of the lens galaxy are consistent with axisymmetry.
Hence, in constructing dynamical models we consider an axisymmetric
stellar system in which both the potential and
distribution function (DF) are independent of azimuth and time.
By Jeans’ (1915) theorem the DF only depends
on the isolating integrals of motion: , with energy , angular momentum parallel to the symmetry -axis, and a third integral
for which in general no explicit expression is known. However,
usually^{12}

Schwarzschild (1979) introduced a method that sidesteps our ignorance about the non-classical integrals of motion. It finds the set of weights of orbits computed in an arbitrary gravitational potential that best reproduces all available photometric and kinematic data at the same time. The method has proved to be powerful in building detailed spherical and axisymmetric models of nearby galaxies (e.g. Rix et al., 1997; van der Marel et al., 1998; Gebhardt et al., 2003; Valluri et al., 2004; Cappellari et al., 2006; Thomas et al., 2007) as well as globular clusters (van de Ven et al., 2006; van den Bosch et al., 2006), and since recently also triaxial models (van de Ven et al., 2008; van den Bosch et al., 2008). However, in all cases the stellar systems are significant closer than the lens galaxy in the Einstein Cross, allowing for (even) more and higher-order stellar kinematic measurements necessary to constrain the large freedom in this general modeling method. Instead we construct simpler, but still realistic dynamical models based on the solution of the axisymmetric Jeans equations.

When we multiply the collisionless Boltzmann equation in cylindrical coordinates by respectively and and integrate over all velocities, we obtain the two Jeans equations (see also Binney & Tremaine, 1987, eq. 4-29)

(1) | |||||

(2) |

where is the intrinsic luminosity density. Due to the assumed axisymmetry, all terms in the third Jeans equation, that follows from multiplying by , vanish.

We are thus left with four unknown second order velocity moments , , and and only two equations. This means we have to make assumptions about the velocity anisotropy, or in other words the shape and alignment of the velocity ellipsoid. In case the velocity ellipsoid is aligned with the cylindrical coordinate system , so that we can readily solve equation (2) for . If we next assume a constant flattening of the velocity ellipsoid in the meridional plane, we can write and solve equation (1) for . This assumption provides in general a good description for the kinematics of real disk galaxies (Cappellari, 2008), such as the lens galaxy in the Einstein Cross, which is an early-type spiral galaxy. When , the velocity distribution is isotropic in the meridional plane, corresponding to the well-known case of a two-integral DF (e.g. Lynden-Bell, 1962; Hunter, 1977).

Knowing the intrinsic second-order velocity moments, the line-of-sight second-order velocity moment for a stellar system viewed at an inclination away from the -axis follows as

(3) | |||||

where is the (observed) surface brightness with the -axis along the projected major axis. For each position on the sky-plane, yields a prediction of the second velocity moment , which is a combination of the (observed) mean line-of-sight velocity and velocity dispersion .

Under the above assumptions, besides the anisotropy parameter (and possibly the inclination ), the only unknown quantity is the gravitational potential. In other words, once the gravitational potential is known, the second velocity moment, together with the assumed velocity anisotropy, define which stellar orbits are present, except for their sense of rotation. Almost any velocity field can then be reproduced by arranging the sense of rotation of the individual orbits, without any change to the gravitational potential. The only limit on is when all the orbits rotate in the same direction, but this limit is not very useful since it is generally much larger than the observed velocity in galaxies. This implies that, unless further assumptions are made, the observed velocity is virtually useless for determining the gravitational potential. Henceforth, we use the combined to constrain the gravitational potential.

The gravitational potential is via Poisson’s equation related to the total mass density . We may estimate the latter from the intrinsic luminosity density , derived from deprojecting the observed surface brightness , once we know the total mass-to-light ratio . It is common in dynamical studies of the inner parts of galaxies to consider as an additional parameter and to assume its value to be constant, i.e., mass follows light (e.g. Cappellari et al., 2006). Since may be larger than the stellar mass-to-light ratio , this still allows for possible dark matter contribution, but with constant fraction. As an alternative to these constant- models, we also construct dynamical models in which we use the strong gravitational lensing to constrain the gravitational potential. In § 3.3 below, we use the positions and relative fluxes of the quasar images to estimate the surface mass density , which we deproject to obtain , so that is not anymore a free parameter. For both type of dynamical models, the best-fits to the stellar kinematic data are shown in Fig. 8, and discussed further in § 4.2 below.

We use the Multi-Gaussian Expansion method
(MGE; Monnet et al., 1992; Emsellem et al., 1994) to
parameterize both the observed surface brightness as well
as the estimated surface mass density as a set of
two-dimensional Gaussians. Even though Gaussians do not form a
complete set of functions, in general surface density distributions
are accurately reproduced, including deviations from an elliptical
distribution and ellipticity variations with radius. Representing also
the point-spread function (PSF) by a sum of Gaussians, the convolution
with the PSF becomes straightforward. Moreover, the
MGE-parameterization has the advantage that the deprojection can be
performed analytically once the viewing angle(s) are given. Also many
intrinsic quantities such as the potential can be calculated by means
of simple one-dimensional integrals.
Similarly, the calculation of in
equation (3) reduces from the (numerical) evaluation
of in general a triple integral to a straightforward single integral
(for further details and a comparison with more general
dynamical models see Cappellari, 2008).
The latter integral is given in Appendix A, together
with expressions for and for the
MGE-parameterization of and as given in
Table 2, and further discussed in
§ 4.1 below.
We used the publically available Jeans Anisotropic MGE (JAM)
implementation^{13}

### 3.3. Lens models

0.8 | 0.0829 | -0.0148 | 2.1615 | -0.0588 | 0.0000 | 0.0018 | 0.0589 | 0.0021 | 0.0029 | 0.8858 | 1.5413 | 1.5998 |

0.9 | 0.0738 | -0.0141 | 1.9457 | -0.0493 | 0.0001 | 0.0011 | 0.0499 | 0.0005 | 0.0019 | 0.8862 | 1.5426 | 1.5845 |

1.0 | 0.0651 | -0.0132 | 1.7733 | -0.0415 | 0.0002 | 0.0006 | 0.0423 | -0.0007 | 0.0012 | 0.8866 | 1.5441 | 1.5756 |

1.1 | 0.0569 | -0.0121 | 1.6325 | -0.0350 | 0.0004 | 0.0003 | 0.0357 | -0.0014 | 0.0008 | 0.8872 | 1.5460 | 1.5708 |

1.2 | 0.0490 | -0.0110 | 1.5153 | -0.0293 | 0.0005 | 0.0001 | 0.0300 | -0.0020 | 0.0006 | 0.8878 | 1.5481 | 1.5688 |

Note. – Parameters of the lens models shown in Fig. 6. For each slope , the best-fit quasar source position and Fourier coefficients and () are given, all in arcseconds. The last three columns provide the Einstein radius (in arcsec), the projected mass within this radius and within the critical curve (both in M).

Under the thin-lens approximation, the gravitational lensing properties of a galaxy are characterized by its potential projected along the line-of-sight, also known as the deflection potential . In case of a MGE parameterization of the surface mass density of the galaxy, the projection along the line-of-sight of the corresponding gravitational potential, or, alternatively, solving the two-dimensional Poisson equation becomes straightforward. The calculation of and its derivatives reduce to the (numerical) evaluation of a single integral, as shown in Appendix A.2 for an axisymmetric system. This MGE approach thus provides a powerful technique to build general lens models. However, since in this case the quasar images only provide a few constraints within a limited radial range, we adopt instead a simpler and less general approach.

We use the algorithm of Evans & Witt (2003) to construct a lens model that accurately fits the (optical) positions and relative (radio) fluxes of the four quasar images in the Einstein Cross. The deflection potential is assumed to be a scale-free function of the polar coordinates and in the lens sky-plane, with for realistic models. The angular part is expanded as a Fourier series

(4) |

The positions of the images are related to the position of the source by the lens equation (e.g. Schneider et al., 1992)

(5) |

where and denotes the derivative to . The flux ratios of the images follow from their magnifications , which are given by

(6) |

where is the parity of image .

Since the Fourier coefficients and correspond to a displacement of the source position , we set to remove this degeneracy. For a given slope , we are left with free parameters and eleven constraints (four image positions and three flux ratios).

For each image position , we obtain through the lens equation (5) a prediction of the source position , which we compare with the assumed free parameter values for the source position by evaluating

(7) |

where the weights in the source position are given by

(8) |

The conversion from the uncertainties in the observed image positions
in the image plane to uncertainties in
the source position in the source plane involves the magnification
tensor^{14}

(9) |

While the expression for the magnification is given in equation (6), those for the partial derivatives of the deflection potential follow as

(10) | |||||

Next, we turn the predicted magnifications through equation (6) per image into predictions for the flux ratios . We compare them to the observed flux ratios with corresponding errors by evaluating

(11) |

We can now find the best-fit parameters by minimizing

(12) |

with a constant. The last term guides the solution towards a smooth, realistic lens model by penalizing (strong) deviations from an elliptic shape.

For a range of slopes , the fits of the above lens model with Fourier terms are shown in Fig. 5. Whereas the predictions of the source position remain within the observed errors over the full range of slopes, the recovery of the observed radio fluxes places constraints on the slope. Dividing the corresponding values by the ten degrees of freedom (eleven constraints minus the slope parameter) provides a quality of fit, shown in the lower-right panel of Fig. 5. The minimum value shows that for a good overall best-fit model is obtained. The corresponding %-confidence interval (dotted line) is .

With Fourier terms the fits are significantly worse (), and with Fourier terms the fit does not improve (for the fitting problem becomes underdetermined). If we replace the radio fluxes by the optical fluxes from the CASTLES website, no acceptable lens model fit is possible (). On the other hand, the mid-infrared fluxes from Agol et al. (2000) yield similar good fits as the radio fluxes. The best-fit slope is lower because the fluxes of image A and C are respectively about 10 % and 30 % lower than in the radio. Image C is the faintest in both mid-infrared and radio, while image D is the faintest in the optical. Since the radio is expected to be (much) less affected by differential extinction and microlensing, we adopt the corresponding lens models.

In the investigation below, we use the lens models with slopes , shown in Fig. 6, and with best-fit parameters given in Table 1. The critical curves are obtained by solving the (quadratic) equation (6) for infinite magnification, i.e., for vanishing left-hand-side. The caustics then follow upon substitution of in the lens equation (5). We see in Fig. 6 that even the two models outside the %-confidence interval for are relatively smooth as expected for a realistic galaxy model.

The surface mass density of the lens models follows from Poisson’s equation as

(13) |

with the critical surface mass density defined as

(14) |

where is the speed of light and , and are the (angular diameter) distance to the lens galaxy, the quasar source and the distance from lens to source, respectively.

The (projected) mass within the critical curve

(15) |

is in general close to the projected mass within the Einstein radius. The latter is the radius for which the projected mass is equal to the Einstein mass defined as . Because is independent of , all higher order Fourier terms average out, and we are left with , so that the Einstein mass becomes

(16) |

Since we express all (projected) coordinates on the plane of the sky in arcseconds, including the radius , we multiply expressions (15) and (16) by to convert from arcseconds to pc for a given (angular diameter) distance to the lens galaxy in Mpc.

From the last two columns in Table 1, we see indeed that M, nearly independent of the slope . Taking into account an inverse scaling with the Hubble constant of km s Mpc, our values are within a few per cent of previous measurements (e.g. Rix et al., 1992; Wambsganss & Paczynski, 1994; Chae et al., 1998; Schmidt et al., 1998; Trott & Webster, 2002).

## 4. Results

### 4.1. Luminous versus total mass distribution

-band surface brightness | lens model | |||||
---|---|---|---|---|---|---|

1 | 4.329 | -1.564 | 0.700 | 5.114 | -1.398 | 0.645 |

2 | 3.935 | -0.942 | 0.700 | 4.556 | -1.037 | 0.660 |

3 | 3.606 | -0.627 | 0.700 | 4.238 | -0.779 | 0.662 |

4 | 3.293 | -0.239 | 0.700 | 3.928 | -0.544 | 0.678 |

5 | 3.005 | -0.043 | 0.700 | 3.628 | -0.350 | 0.673 |

6 | 2.845 | 0.230 | 0.700 | 3.459 | -0.186 | 0.675 |

7 | 2.261 | 0.587 | 0.700 | 3.354 | -0.001 | 0.675 |

8 | 2.160 | 0.917 | 0.414 | 3.208 | 0.220 | 0.675 |

9 | 1.334 | 1.135 | 0.700 | 3.069 | 0.496 | 0.675 |

10 | - | - | - | 2.991 | 1.000 | 0.675 |

Note. – The parameters of the Gaussians in the MGE fit to the HST/WFPC2/F814W -band image of the surface brightness (columns 2–4, “deconvolved” with the WFPC2 PSF), and to the surface mass density of the overall best-fit lens model with slope (columns 5–7) of the lens galaxy in the Einstein Cross (see also Fig. 1). Columns two and five give for each Gaussian component respectively the central surface mass density (in M pc) and the central surface brightness (in L pc), columns three and six the dispersion (in arcsec) along the major axis, and columns four and seven the observed flattening.

In the left and middle panel of Fig. 1 we show MGE fits (contours) to respectively the -band and the -band HST image, obtained with the software of Cappellari (2002), while masking the quasar images. For each fit we assume the same position angle (PA) for the Gaussians, consistent with the adopted axisymmetry in the dynamical models. Although the bar and spirals cannot be reproduced, it provides a good description of the disk in the outer parts (left panel) and reproduces well the bulge in the inner parts (middle panel). Similarly, we derive MGE fits to the surface mass density in equation (13) of the best-fit lens models with slopes from 0.8 to 1.2. In the right panel of Fig. 1, we show the MGE fit (solid contours) to the surface mass density (dashed contours) of the overall best-fit lens model with slope . The corresponding parameters are given in Table 2 (columns 5–7), together with the parameters of the MGE fit to the -band surface brightness (columns 2–4). The latter parameters are “deconvolved” using a MGE model fitted to the WFPC2 PSF as determined with Tiny Tim (Krist & Hook, 2004).

The Gaussians in the MGE fit to the -band surface brightness have a PA of , which is similar to the PA of in the MGE fit to the surface mass density of lens models. Both these values are consistent with the measurements by Yee (1988), who found a PA of for the axis through quasar images C and D, bracketed by a PA of for the outer disk and a PA of for the bar (see also Fig. 1 of Trott & Webster, 2002).

As shown in Fig. 7, not only the orientation, but also the shape of the luminous density distribution as inferred from the surface brightness, is very similar to that of the total mass density distribution as derived from the lens models. The left and middle panel in the top row show respectively the profile and the projected ellipticity of the surface mass density as function of the major axis on the sky-plane. The thin curves are for the lens models with fixed slope (as indicated in the bottom-right panel). The thick solid curve is for the -band surface brightness, multiplied with a constant mass-to-light ratio . The latter is derived from dividing the projected mass from the lens models (shown in the top-right panel of Fig. 7) within the Einstein radius, i.e., the Einstein mass , by the projected luminosity within the Einstein radius .

We see that around (indicated by the dashed vertical line),
where the lens models are best constrained by the observed quasar
image positions and relative flux ratios, both the slope and the
ellipticity of the *projected* luminous and total mass density
are nearly the same. The same holds true for the corresponding
*intrinsic* mass densities, of which the profiles and
ellipticities are shown in the left
and middle panel in the bottom row of Fig. 7.
These (analytic) deprojections (see Appendix A) are
under the assumption of oblate axisymmetry and for a given inclination
which we derive from the flattening of the disk in the outer parts. We
find from the MGE fit to the -band surface brightness (left panel
of Fig. 1), that Gaussian components with a measured
flattening as small as are required for an acceptable
fit. This sets a lower limit to the inclination of (significantly larger than found
by Irwin et al., 1989). Adopting a lower limit for the
intrinsic flattening of the disk of
(e.g. Lambas et al., 1992), we then obtain an inclination of
.

Around the luminosity density is close to isothermal ( intrinsic or in projection) like the mass density of the overall best-fit lens model. Towards the center the slope of the luminosity density becomes shallower. The (fixed) slope of the lens model is not anymore (well) constrained by the quasar images towards the center, neither for radii a few times . At these larger radii also the luminosity density is not only due to the bulge, but the bar and disk start contributing, as can also be seen from the increase in the ellipticity. Nevertheless, within a radius ″, the intrinsic luminous and total mass distribution are very similar, i.e., mass follows light. This means we can infer the gravitational potential either from the observed surface brightness adopting a (constant) total mass-to-light ratio , or directly from the best-fit lens model (see also Appendix A). The corresponding circular velocity curves, defined as in the equatorial plane, are shown in the bottom-right panel of Fig. 7, using as above.

### 4.2. Axisymmetric Jeans models

In Fig. 8, the second panel in the bottom row shows the square root of the combination of the observed (and symmetrized) mean line-of-sight velocity and dispersion maps of the lens galaxy (Fig. 4, two panels on the right). The other panels in Fig. 8 show the line-of-sight second-order velocity moment as predicted by the JAM models (Cappellari, 2008), which are based on the solution of the axisymmetric Jeans equations as summarized in Appendix A. With a central velocity dispersion of km s, the relation as given by Tremaine et al. (2002) predicts a central black hole of mass M for the lens galaxy. Since the corresponding sphere-of-influence ″ is not resolved by our kinematic data, we cannot fit for it, but do include the predicted BH in the Jeans models (as an additional Gaussian). Finally, we convolve the predicted second-order velocity moment with the Gaussian PSF of the kinematic data with FWHM″. The (linear) scale in km s, as indicated by the color bar and limits at the top of the panel with the observations, is the same in all panels of Fig. 8.

The axisymmetric Jeans models in the top row use the gravitational potential inferred from the -band surface brightness, assuming that mass follows light. The corresponding (constant) dynamical mass-to-light ratio that provides the best-fit to the observations is indicated at the top, together with the assumed value for the anisotropy parameter , increasing from left to right. The corresponding increase in the (constant) flattening of the velocity ellipsoid in the meridional plane (or decrease in ), has two main effects: a less pinched “butterfly” shape and a weaker gradient parallel to the (vertically aligned) short axis. Comparing with the observations, we see that at lower values the pinching is too strong, while at higher the gradient along the short axis is too shallow. Therefore, it is not unexpected that the “best-fit” is when and (in the -band), but with reduced significantly above unity.

Part of the latter mismatch might be due to (systematic) discrepancies in the data, in particular due to the challenging measurement of the velocity dispersion. The average uncertainty in the observed second-order velocity moment is about km s, which translates into an error of around per cent in , or about . However, the uncertainties in the observed second-order velocity moment vary quite a lot, from only km s in the center to km s a the edge of the GMOS field. And although we were very careful in disentangling the contribution from the quasar and excluded the regions at the quasar images, the measured velocity dispersion around these regions might still be systematically affected by (broadening due to) residual light from the quasar source. This possibly explains the observed excess in the observed second-order velocity moment around the quasar images with respect to the accurately measured central value. At the same time, such an excess might mimic a steeper gradient in the observed second-order velocity moment. Indeed, in all Jeans models in the top row of Fig. 8 the prediction at the center is significantly higher than the observed value. As a result, we expect that the dynamical mass-to-light ratio is overestimated.

If, alternatively, we infer the gravitational potential directly, i.e., without the need for a mass-to-light ratio, from the (deprojected) surface mass density of the best-fit lens models, the axisymmetric Jeans models predict maps of the second-order velocity moments that are similar in shape but with lower values near the center. This is shown in the middle row of Fig. 8 for the overall best-fit lens model with slope , and for the same range in anisotropy parameters . The formal best-fit is obtained for , but the observed central value is more closely matched in case of mild velocity anisotropy in the meridional plane with . For . i.e., density slopes steeper then isothermal, the observed central value is best matched when , with at the same time a steeper gradient along the short-axis, as shown in the first panel in the bottom row of Fig. 8. The last two panels show that for , velocity isotropy in the meridional plane () results in a central value much lower than the observed second-order velocity moment.

We conclude that for the lens models with slopes , as constrained by the lensing geometry, and for a mild anisotropic velocity distribution in the meridional plane with , the predictions of the second-order velocity moments match the accurately measured central value, and follow the gradient within the observational uncertainties outside the region that might be affected by (residual) quasar light. The latter results in broadening of the velocity dispersion, so that the best-fit dynamical mass-to-light ratios we derive when we infer the gravitational potential from the surface brightness, are likely overestimated. Fitting instead only the unaffected and accurately measured second-order velocity moment in the center, we find for a best-fit value of , with an error around . On the other hand, both the Einstein mass and luminosity within the Einstein radius are well constrained, resulting in at most a few per cent error in . Therefore, in what follows, we adopt the latter as the total mass-to-light ratio in the inner bulge-dominated region of the lens galaxy.

In the above axisymmetric Jeans models we fix the inclination to determined from the observed ellipticity of the outer disk (see § 4.1). When we vary the inclination over a reasonable range above the lower limit , there is no significant change in the above results, which is expected because the inclination only has a weak effect on the mass-to-light ratio (see also Fig. 4 of Cappellari et al., 2006).

### 4.3. Stellar mass-to-light ratio

We have shown that in the bulge-dominated inner region of the lens
galaxy, the total mass distribution closely follows the light
distribution. This does not mean there is no dark matter present,
since if dark matter follows (nearly) the same distribution it might
contribute to part of the above estimated total mass-to-light ratio
. To determine this possible
*constant* dark matter fraction, we need to measure the stellar
mass-to-light ratio . To this end we use the single stellar
population (SSP) models of Vazdekis et al. (1996), adopting as
initial reference the Kroupa (2001) initial mass function
(IMF) with lower mass cut-off M(faintest star
M) and upper mass cut-off of M.

We derive a high-S/N spectrum from our data cube by collapsing the
central spectra that are unaffected by the quasar images. Even so, the
Ca II triplet is equally well fitted by SSP models of nearly all ages
and metallicities of about around solar, which leaves
the stellar mass-to-light ratio nearly unconstrained. Therefore, we
also derived colors from archive HST images (PI: Kochanek).
After matching all images, using both the quasar positions and
isophotes of the lens galaxy, we compute the magnitudes in a circular
aperture of radius ″ (well within the quasar images at
″).
After correcting for a Galactic extinction of
(Schlegel et al., 1998), we derive in the F555W, F675W, F814W
(WFPC2), F160W and F205W (NICMOS2) filters magnitudes of 16.83, 16.03,
15.53, 13.80, and 13.51, respectively. We convert the latter to
*rest frame* magnitudes of 16.72, 16.13, 15.50, 13.70, 13.59 in
Johnson , , , and filters, respectively. Instead of
the standard K-correction, we use the method described in § 3 of
van de Ven et al. (2003) to take simultaneously into account the
filter responses and the redshift of the lens galaxy.

The resulting colors are best fitted by a SSP model with age Gyr and metallicity , corresponding to an -band stellar mass-to-light ratio . The confidence limits yield a range in from 7 to 14 Gyr, and in from to , corresponding to a range in from to . The lower limit implies at most per cent of dark matter within ″, but the results are fully consistent with no dark matter at all. The SSP models of Maraston (1998, 2005) with a Kroupa IMF yield similar results. Adopting a Salpeter (1955) IMF with the same lower and upper mass cut-offs, increases by about per cent. This not only implies no dark matter, but even for the lower limit , which is unphysical. This is consistent with evidences for both late-type galaxies (Bell & de Jong, 2001) and early-type galaxies (Cappellari et al., 2006) that, if the IMF is universal, a Salpeter IMF is excluded, whereas a Kroupa IMF matches the observations (see de Jong & Bell, 2007, for a review).

## 5. Discussion and conclusions

We used the GMOS-North integral-field spectrograph to obtain two-dimensional stellar kinematics of the lens galaxy in the Einstein Cross. In addition to the four bright quasar images and the distance of the lens galaxy ( Mpc), in particular the presence of sky lines in the observed Ca II triplet region made the extraction of the absorption line kinematics challenging. Even so, we were able to derive high-quality line-of-sight velocity and dispersion maps of the bulge-dominated inner region ″, reaching about two-thirds of the effective radius ″ of this early-type spiral galaxy. The map shows regular rotation up to km s around the minor axis of the bulge, consistent with axisymmetry. The map shows a weak gradient increasing towards a central (″) value of km s.

The only other direct (single) measurement of the velocity dispersion is km s by Foltz et al. (1992). Adopting their central aperture of 0.7″0.4″ at a position angle of 39 (along the bar), our extracted spectrum is shown in Fig. 2. The lower panel shows sky-subtracted galaxy spectrum (black curve) and our best-fit composite stellar population model (blue curve), yielding a velocity dispersion in this aperture of 172 km s. The other model (red curve) has the velocity dispersion fixed to 215 km s as measured by Foltz et al. (1992). Since the observing conditions were similar with seeing around ″, we can directly compare both measurements. We see from residuals at the bottom and the quoted reduced -values that 172 km sprovides a better fit than 215 km s, particularly in fitting the strongest of the three absorption lines.

In addition, for a singular isothermal sphere lens model, we can use the relation (e.g. Kochanek et al., 2000) with a separation ″ of the four quasar images, to obtain a simple estimate for the dispersion of km s. Adopting a spherical Jeans model with a range in velocity anisotropy, van de Ven et al. (2003) convert this to a central velocity dispersion of km s within a circular “Coma” aperture (with a diameter of 3.4″at the distance of the Coma cluster). At the (angular diameter) distance Mpc of the lens galaxy, this corresponds to a circular aperture with radius of 1″, so that it can be compared directly to our central value of km s. Finally, King and de Vaucouleurs models of Kent & Falco (1988) predict a similar value of km s, and also Barnes et al. (1999) find a value of km s based on their two H i rotation curve measurements. Hence, it is likely that the long-slit measurement by Foltz et al. (1992) is affected by the bright quasar images, whereas our spatially resolved measurements allow for a clean(er) separation of the quasar contribution.

A large variety of different lens models have been constructed for the Einstein Cross, most of which fit the positions of the quasar images but not their relative flux ratios. Adopting the scale-free lens model of Evans & Witt (2003), we found that fitting at the same time also the (radio) flux ratios constrained the slope of the total mass surface density to be . The total mass within the Einstein radius ″, i.e., the Einstein mass, is M, nearly independent of the slope , and consistent with previous measurements in the literature. Dividing by the projected luminosity within , as measured from the observed -band surface brightness, we obtained a mass-to-light ratio of , with an error of at most a few per cent.

We determined the mass-to-light ratio in an additional, independent way by fitting dynamical models to the (combined) observed and maps. We used the solution of the axisymmetric Jeans equations to predict this second-order velocity moment, with the gravitational potential inferred from the deprojected -band surface brightness, assuming that mass follows light. We expect the resulting best-fit constant mass-to-light ratio to be an overestimation due to possible residual contribution from the quasar light. Even so, restricting the fit to the unaffected and accurately measured central second-order velocity moment, we found . When we used instead the gravitational potential that follows directly, i.e., without mass-to-light ratio conversion, from the surface mass density of the best-fit lens models, we arrived at a similar prediction of the second-order velocity moment. We showed that the reason is that the luminous and total mass distribution, as inferred from respectively the surface brightness and lens model, are very similar. This implies that the inner region of the lens galaxy mass (closely) follows light, with a total mass-to-light ratio .

By fitting single stellar population models to measured colors of the
center of the lens galaxy, we estimated an -band stellar
mass-to-light ratio from 2.8 to 4.1 . Although a
*constant* dark matter fraction of 20 per cent is thus not
excluded, it is likely that dark matter does not play a significant
role in the inner region of this early-type spiral galaxy of
luminosity . This is consistent with indications that
less-luminous early-type galaxies (e.g. Gerhard et al., 2001; Rusin et al., 2003; Cappellari et al., 2006; Thomas et al., 2007; Bolton et al., 2008) and late-type galaxies
(e.g. Persic et al., 1996; Palunas & Williams, 2000; Bell & de Jong, 2001; Kassin et al., 2006; Williams et al., 2009) are
dominated by the stellar mass inside their central regions. This is
different from dwarf galaxies as well as in (at least the outer parts
of) giant elliptical and spiral galaxies where dark matter is expected
to be ubiquitous.

The constraint on the slope of the lens model implies that the intrinsic total mass density is close to isothermal, consistent with previous studies of lens galaxies (e.g. Koopmans et al., 2006, 2009). However, one has to be very careful not to over-interpret this result, not only due to assumptions in the modeling, but most of all because the slope is only (well) constrained around the Einstein radius . If indeed mass follows light in the inner region of this lens galaxy, the (deprojected) surface brightness indicates deviations from isothermal towards the center where the slope becomes shallower as well as a possible steepening at larger radii. Interestingly, the construction of realistic galaxy density profiles with a stellar and dark matter component that are both non-isothermal, shows that the combined slope and corresponding lensing properties are nevertheless consistent with isothermal around (van de Ven et al., 2009). With more extended images, e.g. in the case of galaxy-galaxy lensing, one can place a stronger constraint on the total mass distribution of the lens galaxy (e.g. Barnabè & Koopmans, 2007; Suyu et al., 2009). Since the extension in these cases is still mostly tangential, constraining a large radial range is in particular possible when lensing occurs in multiple (redshift) planes, resulting in images at different Einstein radii (e.g. Gavazzi et al., 2008).

Nevertheless, already the quasar images provide an accurate constraint on the total mass within , nearly independent of the details of the lens model (e.g. Kochanek, 1991; Evans & Witt, 2001). Therefore, at least around this radius no higher-order velocity moments are needed to break the mass-anisotropy degeneracy. Towards the center the surface brightness increases steeply, so that high enough S/N to measure velocity moments beyond and might be achievable even at higher redshift. In the outer parts, the degeneracy might be (partially) broken by combining the kinematics of stars and/or discrete tracers, such as globular clusters and planetary nebulae, with total mass estimates from hot X-ray gas and/or weak lensing measurement (e.g. Gavazzi et al., 2007). In the current and even more in the upcoming extensive and deep photometric surveys, numerous (strong) gravitational lensing systems will be discovered. They provide important and independent constraints on the total mass distribution in galaxies, especially in combination with (resolved) stellar kinematics, as we showed in this paper (see also Barnabè & Koopmans, 2007; Czoske et al., 2008; Barnabè et al., 2009).

## Appendix A Multi-Gaussian Expansion

We summarize the (numerically) convenient expressions of the axisymmetric intrinsic density and gravitational potential in the case of a Multi-Gaussian Expansion (MGE) of the surface density (Monnet et al., 1992; Emsellem et al., 1994). Next, we show that also the gravitational lensing properties can be readily computed, so that the MGE method can be used to efficiently construct general lens models. Finally, we give the line-of-sight second-order velocity moment as a solution of the axisymmetric Jeans equations which is discussed in detail in Cappellari (2008).

### a.1. Axisymmetric density and potential

We parameterize the surface brightness by a sum of Gaussian components

(A1) |

each with three parameters: the central surface brightness , the dispersion along the major -axis and the flattening . In case of an (oblate) axisymmetric system viewed at an inclination , the corresponding intrinsic luminosity density is

(A2) |

with intrinsic dispersion and intrinsic flattening given by .

The intrinsic mass density follows as , where the mass-to-light ratio per Gaussian component is a free parameter. Alternatively, if the surface mass density is available (e.g. from gravitational lensing), an MGE-parameterization as in (A1) directly yields as in (A2), but with in both expressions replaced by the central surface mass density .

The corresponding gravitational potential follows upon (numerical) evaluation of (Emsellem et al., 1994, eq. 39)

(A3) |

where we have introduced

(A4) |

and . The total mass per Gaussian component is given by , with when the mass density is inferred from the surface brightness.

In the latter case, we have implicitly assumed a *constant*
mass-to-light ratio per Gaussian component by taking it outside
the integral. Nevertheless, we can still mimic a (radially) varying
mass-to-light ratio by considering the of the Gaussian
components as free parameters (see e.g. van de Ven et al., 2006; van den Bosch et al., 2006). However, it is common in dynamical studies of
the inner parts of galaxies (e.g. Cappellari et al., 2006) to
assume the total mass-to-light ratio to be constant, i.e., for each Gaussian component . Since may be larger
than the stellar mass-to-light ratio , this still allows for
possible dark matter contribution, but with a constant fraction.

### a.2. Axisymmetric lens model

Under the thin-lens approximation, the gravitational lensing properties of a galaxy are characterized by the deflection potential and its (partial) derivatives. The deflection potential follows from projecting the potential along the line-of-sight or by solving the two-dimensional Poisson equation . Here, is the normalized surface mass density with the critical lensing value defined in equation (14).

Similar to the surface brightness in equation (A1), we parameterize this so-called “convergence” by a sum of Gaussian components

(A5) |

The corresponding deflection potential is then

(A6) |

where is the total mass per Gaussian component, and we have introduced

(A7) |

and .

The first-order partial derivatives of follow as

(A8) |

whereas the second-order partial derivatives are given by

(A9) |

The above single integrals can be readily evaluated numerically, so that the lensing properties follow in a straightforward way. The image positions are related to the source position by the lens equation