Twodimensional kinematics of SLACS lenses – II.
Combined lensing
and dynamics analysis of earlytype galaxies
at
Abstract
We present the first detailed analysis of the mass and dynamical structure of a sample of six earlytype lens galaxies, selected from the Sloan Lens ACS Survey, in the redshift range . Both Hubble Space Telescope (HST)/ACS highresolution imaging and VLT VIMOS integral field spectroscopy are available for these systems. The galaxies are modelled—under the assumptions of axial symmetry and twointegral stellar distribution function—by making use of the cauldron code, which selfconsistently combines gravitational lensing and stellar dynamics, and is fully embedded within the framework of Bayesian statistics. The principal results of this study are: (i) all galaxies in the sample are well described by a simple axisymmetric powerlaw profile for the total density, with a logarithmic slope very close to isothermal ( and an intrinsic spread close to per cent) showing no evidence of evolution over the probed range of redshift; (ii) the axial ratio of the total density distribution is rounder than and in all cases, except for a fast rotator, does not deviate significantly from the flattening of the intrinsic stellar distribution; (iii) the dark matter fraction within the effective radius has a lower limit of about to per cent; (iv) the sample galaxies are only mildly anisotropic, with ; (v) the physical distinction among slow and fast rotators, quantified by the ratio and the intrinsic angular momentum, is already present at . Altogether, earlytype galaxies at are found to be markedly smooth and almost isothermal systems, structurally and dynamically very similar to their nearby counterparts. This work confirms the effectiveness of the combined lensing and dynamics analysis as a powerful technique for the study of earlytype galaxies beyond the local Universe.
keywords:
gravitational lensing — galaxies: elliptical and lenticular, cD — galaxies: kinematics and dynamics — galaxies: structure.1 Introduction
The currently favoured cosmological scenario, the socalled CDM (cold dark matter) paradigm, has been remarkably successful at explaining the large scale structure of the Universe. In the nonlinear regime, below several Mpc, however, the situation is less certain, and a full understanding of the galaxy formation and evolution processes remains a work in progress.
Within the standard paradigm, massive earlytype galaxies are thought to be the endproduct of hierarchical merging of lower mass galaxies, and to be embedded in extended dark matter haloes (e.g. Toomre, 1977; White & Frenk, 1991; Barnes, 1992; Cole et al., 2000). Numerical studies of merging galaxies (e.g. Naab et al., 2006; Jesseit et al., 2007) have managed to reproduce a number of observational characteristics of massive ellipticals, and have made clear that stringent tests of galaxy formation models require a detailed and reliable description of the intrinsic physical properties of real earlytype galaxies, such as their mass density distribution and orbital structure. Furthermore, knowledge of how these galaxy properties evolve through time would provide much needed information and even stronger constraints on the theoretical predictions.
Galaxy name  

(km s)  (kpc)  (mag)  (kpc)  (mag)  (deg)  (kpc)  ()  
SDSS J00370942  0.1955  0.6322  23.11  22.25  0.76  9.5  4.77  27.3  
SDSS J02160813  0.3317  0.5235  23.94  23.06  0.85  79.2  5.49  48.2  
SDSS J09120029  0.1642  0.3240  23.41  22.55  0.67  13.2  4.55  39.6  
SDSS J09590410  0.1260  0.5349  21.45  20.58  0.68  57.4  2.25  7.7  
SDSS J16270053  0.2076  0.5241  22.68  21.71  0.85  5.6  4.11  22.2  
SDSS J23210939  0.0819  0.5324  22.59  21.72  0.77  126.5  2.43  11.7 
In the last decades, local earlytype galaxies have been the object of substantial observational and modelling efforts. These studies have employed a variety of tracers ranging from stellar kinematics (e.g. Saglia, Bertin & Stiavelli 1992, Bertin et al. 1994, Franx et al. 1994, Carollo et al. 1995, Rix et al. 1997, Loewenstein & White 1999, Gerhard et al. 2001, Borriello et al. 2003, Thomas et al. 2007 and the SAURON collaboration: see e.g. de Zeeuw et al. 2002, Emsellem et al. 2004, Cappellari et al. 2006) and kinematics of discrete tracers such as globular clusters (e.g. Mould et al., 1990; Côté et al., 2003) or planetary nebulae (e.g. Arnaboldi et al. 1996, Romanowsky et al. 2003; also in combination with the kinematics of stars: de Lorenzi et al. 2008) to hot Xray gas (e.g. Fabbiano, 1989; Matsushita et al., 1998; Fukazawa et al., 2006; Humphrey et al., 2006), usually finding evidence for a dark matter halo component and for a total mass density profile close to isothermal (i.e. ) in the inner regions.
On the other hand, because of the severe observational limitations, thorough studies of distant earlytype galaxies (at redshift ) are still in their infancy. Traditional analyses based on dynamics alone are hindered by the lack of tracers at large radii and by the mass–anisotropy degeneracy, i.e. a change in the mass profile of the galaxy or in the anisotropy of the velocity dispersion tensor can both determine the same effect in the measured velocity dispersion map. Higherorder velocity moments, which potentially allow one to disentangle this degeneracy by providing additional constraints (Gerhard, 1993; van der Marel & Franx, 1993), can only be measured with sufficient accuracy in the inner parts of nearby galaxies with the current instruments. Fortunately, valuable additional information on distant earlytype galaxies can be provided by gravitational lensing (see e.g. Schneider, Ehlers & Falco, 1992), when the galaxy happens to act as a gravitational lens with respect to a luminous background source at higher redshift. Strong gravitational lensing allows one to determine the total mass within the Einstein radius in an accurate and almost modelindependent way (Kochanek, 1991) although, due to the masssheet or massslope degeneracies (Falco, Gorenstein & Shapiro, 1985; Wucknitz, 2002), it does not permit the univocal recovery of the mass density profile of the lens.
Gravitational lensing and stellar dynamics are particularly effective when they are applied in combination to the analysis of distant earlytype galaxies. The complementarity of the two approaches is such that the masssheet and mass–anisotropy degeneracies can to a large extent be disentangled and the mass profile of the lens galaxy can be robustly determined (see e.g. Koopmans & Treu 2002, 2003, Treu & Koopmans 2002, 2003, 2004, Barnabè & Koopmans 2007, Czoske et al. 2008, Czoske, Barnabè & Koopmans 2008, van de Ven et al. 2008, Trott et al. 2008). Recently, the dedicated Sloan Lens ACS Survey (SLACS; Bolton et al., 2006; Treu et al., 2006; Koopmans et al., 2006; Gavazzi et al., 2007; Bolton et al., 2008; Gavazzi et al., 2008; Bolton et al., 2008; Treu et al., 2008) has discovered a large and homogeneous sample of 70 strong gravitational lenses^{1}^{1}1A further 19 systems are possible gravitational lenses, but the multiple imaging is not secure., for the most part earlytype galaxies, at redshift between and . Koopmans et al. (2006) have applied a joint analysis to the SLACS galaxies, finding an average total mass profile very close to isothermal with no sign of evolution to redshifts approaching unity (when including also systems from the Lenses Structure and Dynamics Survey, e.g. Treu & Koopmans 2004). A limitation of the method used in those works is that lensing and dynamics are treated as independent problems, and all the kinematic constraints come from a single apertureaveraged value of the stellar velocity dispersion, which could potentially lead to biased results. Barnabè & Koopmans (2007, hereafter BK07) have expanded the technique for the combined lensing and dynamics analysis into a more general and selfconsistent method, embedded within the framework of Bayesian statistics. Implemented as the cauldron algorithm, under the only assumptions of axial symmetry and twointegral stellar distribution function (DF) for the lens galaxy, this method makes full use of all the available data sets (i.e. the surface brightness distribution of the lensed source, and the surface brightness and kinematic maps of the lens galaxy) in order to recover the lens structure and properties in the most complete and reliable way allowed by the data. The ideal testingground for an indepth analysis is represented by a subsample of 17 SLACS lenses which have been observed with the VIMOS integralfield unit (IFU) mounted on the VLT as part of a pilot and a large programme (ESO programmes 075.B0226 and 177.B0682, respectively; PI: Koopmans), obtaining detailed twodimensional kinematic maps (first and second velocity moments) in addition to the HST imaging data. The first joint study conducted with the cauldron code of one of these systems, SDSS J2321 at , has been presented in Czoske et al. (2008, hereafter C08).
In this paper, we extend the study of C08 to a total of six SLACS lenses for which kinematic data sets are now available: SDSS J0037, SDSS J0216, SDSS J0912, SDSS J0959 and SDSS J1627, in addition to the already mentioned SDSS J2321. This sample was chosen to cover a range in redshift, mass and importance of rotation which is representative of the SLACS sample. Since the SLACS galaxies have been shown to be statistically indistinguishable from control samples in terms of any of their known observables, such as size, luminosity, surface brightness (Bolton et al., 2008), location on the Fundamental Plane (Treu et al., 2006) and environment (Treu et al., 2008), we expect that the results of the combined lensing and dynamics analysis described in this work can be generalized to the massive earlytype population, nicely complementing the work done by, e.g., the SAURON collaboration on lower redshift and lower mass earlytype galaxies. Basic information on the six systems under study are listed in Table 1. The VIMOS and HST observations of these systems, together with a description of the data reduction, will be detailed in a forthcoming paper (Czoske et al., in preparation).
This paper is organized as follows: in Section 2 we give a brief overview of the available data sets. In Section 3 we recall the basic features of the cauldron algorithm and the adopted mass model. The results of the combined lensing and dynamics analysis of the SLACS subsample are presented in Sections 4 and 5, with the latter focusing on the recovered dynamical structure of the lenses. In Section 6 we summarize our findings and draw conclusions. Throughout this paper we adopt a concordance CDM model described by , and with , unless stated otherwise.
2 Overview of the data sets
Galaxy  (s)  Grism  [Å]  

SDSS J0037  33  18 315  HR_Blue  
SDSS J0216  14  28 840  HR_Orange  
SDSS J0912  12  6 660  HR_Blue  
SDSS J0959  5  10 300  HR_Orange  
SDSS J1627  12  24 720  HR_Orange  
SDSS J2321  15  8 325  HS_Blue 
2.1 Spectroscopy
Integralfield spectroscopy for seventeen lens systems was obtained using the integralfield unit of VIMOS on the VLT, UT3. All observations, split in Observing Blocks (OB) of roughly one hour, including calibration, were done in service mode.
Of the six systems analyzed in this paper, three were observed in the course of a normal ESO programme, a pilot, 075.B0226 (PI: Koopmans). For this programme we used the HRBlue grism with a resolution of Å ( Å full width at half maximum, FWHM), covering an observed wavelength range of 4000 to 6200 Å. Each OB was split into three dithered exposures of seconds each. For the large programme, we switched to the HROrange grism with mean resolution Å, covering the range 5050 to 7460 Å. Only one long exposure of seconds was obtained for each observing block; the number of observing blocks was sufficient to fill in on gaps in the data due to bad instrument fibres through pointingoffset between subsequent OBs.
The data were reduced using the VIPGI package (Scodeggio et al., 2005; Zanichelli et al., 2005). For more details on the procedure and tests of the quality of the reduced data we refer to Czoske et al. (2008) and Czoske et al. (in preparation).
The kinematic parameters and were determined from the individual spectra using a direct pixelfitting routine. Compared to Czoske et al. (2008), we have made a number of modifications in the algorithm. In particular, we now use almost the entire wavelength range that is available from the spectra; noisy parts at the blue and red ends of the spectra were cut off. Due to the varying redshifts of the lenses, the restframe wavelength ranges and hence the spectral features that were used in the kinematic analysis varied from lens to lens. The template used was a spectrum of the K2 giant HR 19, taken from the IndoUS survey (Valdes et al., 2004). The native resolution of the template spectrum is Å FWHM ( Å). The template is first smoothed to the instrumental resolution of the VIMOS spectra, corrected to the rest frame of the lens. Due to the low signaltonoise ratio of the spectra from individual spaxels, we assume the lineofsight velocity distribution (LOSVD) to be described by a Gaussian. Tests show that including GaussHermite terms and (van der Marel & Franx, 1993) does not improve the fit in terms of per degree of freedom and does not yield robust results for and and consequently for the velocity moments and (see also Cappellari & Emsellem, 2004). The larger wavelength range requires us to modify the linear correction function used in Czoske et al. (2008) by multiplicative and additive polynomial corrections (following Kelson et al., 2000). Extensive testing shows that choosing polynomial orders of five ensures good fits for the continuum shape without affecting the structure of the small scale absorption features.
A number of spectral features that are not well reproduced by the stellar template are masked. This includes in particular the Mg b line which is enhanced in the lens galaxy, possibly the result of an enhancement, as compared to the Galactic star HR 19 (Barth et al., 2002) and the Balmer lines which may be partially filled in by emission.
2.2 Imaging
We use Hubble Space Telescope (HST) imaging data from ACS and NICMOS to obtain information on the surface brightness distributions of the lens galaxies and the gravitationally lensed background galaxies.
ACS images taken through the F814W filter form the basis of the lens modelling. For four of the systems, we have deep (fullorbit) imaging (program 10494, PI: Koopmans); for the remaining two (SDSS J2321, SDSS J0037) we use singleexposure images from our snapshot program (10174, PI: Koopmans). Nonparametric elliptical Bspline models of the lens galaxies were subtracted off the images in order to obtain a clean representation of the source structure without contamination from the lens galaxy (Bolton et al., 2008).
For the dynamical analysis the kinematic maps are weighted by the surface brightness of the lens galaxy. We use the reddest band possible since this gives the most reliable representation of the stellar light. For four of the six systems described here, NICMOS images taken through the F160W filter are available. For SDSS J2321 and SDSS J1627, we had instead to resort to the F814W ACS images. Since the lensed source is typically stronger in F814W than in F160W, we start from the Bspline model of the lens galaxy to which we add random Gaussian noise according to the variance map of the images. The images are convolved to the spatial resolution of the VIMOS data (the seeing limit of the service mode observations, ) and resampled to the grid of the kinematic data using swarp^{2}^{2}2http://terapix.ia.fr/soft/swarp (Bertin, 2008).
3 Combining lensing and dynamics
In this Section we recall the main features of the cauldron algorithm and describe the adopted family of galaxy models. The reader is referred to BK07 for a detailed description of the method.
3.1 Overview of the cauldron algorithm
The central premise of a selfconsistent joint analysis is to adopt a total gravitational potential (or, equivalently, the total density profile , from which is calculated via the Poisson equation), parametrized by a set of nonlinear parameters, and use it simultaneously for both the gravitational lensing and the stellar dynamics modelling of the data. While different from a physical point of view, these two modelling problems can be formally expressed in an analogous way as a single set of coupled (regularized) linear equations. For any given choice of the nonlinear parameters, the equations can be solved (in a direct, noniterative manner) to obtain as the best solution for the chosen potential model: (i) the surface brightness distribution of the unlensed source, and (ii) the weights of the elementary stellar dynamics building blocks (e.g. orbits or twointegral components, TICs, Schwarzschild, 1979; Verolme & de Zeeuw, 2002). This linear optimization scheme is consistently embedded within the framework of Bayesian statistics. As a consequence, it is possible to objectively assess the probability of each model by means of the evidence merit function (also called the marginalized likelihood) and, therefore, to compare and rank different models (see MacKay, 1992, 1999, 2003). In this way, by maximizing the evidence, one recovers the set of nonlinear parameters corresponding to the “best” potential (or density) model, i.e. that model which maximizes the joint posterior probability density function (PDF), hence called maximum a posteriori (MAP) model. In the context of Bayesian statistics, the MAP model is the most plausible model in an Occam’s razor sense, given the data and the adopted form of the regularization (the optimal level of the regularization is also set by the evidence). However, in a Bayesian context, the MAP parameters individually do not necessarily have to be probable: they have a high joint probability density, but might only occupy a small volume in parameter space. This situation can arise if the MAP solution does not lie in the bulk of the posterior probability distribution (roughly parameterspace volume times likelihood density). This can be particularly severe if the PDF is nonsymmetric in a highdimensional space where volume can dramatically increase with distance from the MAP solution. In that case, the MAP solution for each parameter could easily lie fully outside the PDFs of each individual parameter, when marginalized over all other parameters. We will see this happening in some cases, and, although somewhat counterintuitive, it is fully consistent with Bayesian statistics and a peculiar aspect of statistics in highdimensional spaces.
As discussed in BK07, the method is extremely general and can in principle be applied to any potential shape. However, its current practical implementation, the cauldron algorithm, is more restricted in order to make it computationally efficient and applies specifically to axially symmetric potentials, , and twointegral DFs , where and are the two classical integrals of motion, i.e., respectively, energy and angular momentum along the rotation axis. Under these assumptions, the dynamical model can be constructed by making use of the fast Monte Carlo numerical implementation by BK07 of the twointegral Schwarzschild method described by Cretton et al. (1999) and Verolme & de Zeeuw (2002), whose building blocks are not stellar orbits (as in the classical Schwarzschild method) but TICs.^{3}^{3}3A TIC can be visualized as an elementary toroidal system, completely specified by a particular choice of energy and axial component of the angular momentum . TICs have simple radial density distributions and analytic unprojected velocity moments, and by superposing them one can build models for arbitrary spheroidal potentials (cf. Cretton et al., 1999): all these characteristics contribute to make TICs particularly valuable and inexpensive building blocks when compared to orbits. The weights map of the optimal TIC superposition that best reproduces the observables, in a Bayesian sense, is yielded as an outcome of the joint analysis.
The code is remarkably robust and its applicability is not drastically limited by these assumptions. Barnabè et al. (2008) have tested cauldron on a galaxy model (comprising a stellar component and a dark matter halo) resulting from a numerical Nbody simulation of a merger process, i.e. a complex nonsymmetrical system which departs significantly from the algorithm’s assumptions. Despite this, it is found that several important global properties of the system, including the total density slope and the dark matter fraction, are reliably recovered.
3.2 The galaxy model
Koopmans et al. (2006) have shown that a powerlaw model, despite its simplicity, seems to provide a satisfactory description for the total mass density profile of the inner regions of SLACS lens galaxies, to the level allowed by the data. This has been further confirmed by the C08 study of SDSS J2321, the first case where a fully selfconsistent analysis was performed. Therefore, in the present work, we still adopt a powerlaw model. If this description is oversimplified for the galaxy under analysis, this will usually have a clearly disruptive effect on the quality of the reconstruction, such as very large residuals (compared to the noise level) for the best model lensed image, and a patchy or strongly pixelized reconstructed source (Barnabè et al., 2008).
Galaxy name  

(deg)  (deg)  (  
SDSS J00370942  0.434  1.968  0.693  0.10  0.23  3.35  5.40  
SDSS J02160813  0.344  1.973  0.816  0.19  0.17  12.20  9.35  
SDSS J09120029  0.412  1.877  0.672  0.16  0.30  7.40  9.08  
SDSS J09590410  0.323  1.873  0.930  0.25  0.30  0.95  7.17  
SDSS J16270053  0.369  2.122  0.851  0.11  0.21  2.23  5.93  
SDSS J23210939  0.468  2.061  0.739  0.13  0.29  1.98  5.22 
The total mass density distribution of the galaxy is taken to be a powerlaw stratified on axisymmetric homoeoids:
(1) 
where is a density scale, will be referred to as the (logarithmic) slope of the density profile, and
(2) 
where and are lengthscales and .
The (inner) gravitational potential associated with a homoeoidal density distribution is given by Chandrasekhar (1969) formula. In our case, for , one has
(3) 
while for
(4) 
where and
(5) 
There are three nonlinear parameters in the potential to be determined via the evidence maximization: (or equivalently, through equation [B4] of BK07, the adimensional lens strength ), the logarithmic slope and the axial ratio . A core radius can be straightforwardly included in the density distribution, if necessary. Beyond these parameters, there are four additional parameters which determine the geometry of the observed system: the position angle , the inclination and the coordinates of the lens galaxy centre with respect to the sky grid. Usually, the position angle and the lens centre (which are very well constrained by the lensed image brightness distribution) can be accurately determined by means of fast preliminary explorations and kept fixed afterwards in order to reduce the number of free parameters during the more computationally expensive optimization and error analysis runs. Finally, a proper modelling of the lensed image can occasionally require the introduction of two additional parameters, namely the shear strength and the shear angle, in order to account for external shear.
We employ a curvature regularization (defined as in Suyu et al. 2006 and Appendix A of BK07) for both the gravitational lensing and the stellar dynamics reconstructions. The level of the regularization is controlled by three socalled “hyperparameters” (one for lensing and two for dynamics, see discussion in BK07), whose optimal values are also set via maximization of the Bayesian evidence. The starting values of the hyperparameters are chosen to be quite large, since the convergence to the maximum is found to be faster when starting from an overregularized system.
4 Analysis and results
4.1 Best model reconstruction
We have applied the combined cauldron code to the analysis of the six SLACS lens galaxies in our current sample with available kinematic maps. The recovered nonlinear parameters for each galaxy best model are presented in Table 3. The listed parameters are the inclination (expressed in degrees), the lens strength , the logarithmic density slope and the axial ratio of the total density distribution. The Bayesian statistical errors on the parameters (i.e. the posterior probability distributions) are presented in Section 4.3.
Our analysis shows that, given the current data, there is no need to include external shear or core radius in the modelling of any of the six galaxies: there is no significant improvement in the evidence when these parameters are allowed to vary, and their final values are found to be very close to zero. As mentioned in Section 3, for each system the lens centre and position angle are evaluated in a preliminary run and then kept constant. The best model position angles, relative to the total mass distribution, are found to depart less than (and, with the exception of SDSS J0959 and SDSS J2321, less than ) from the observed values obtained from the light distribution. This suggests that there is at most a small misalignment between the dark and luminous components. A similar conclusion was also drawn in Koopmans et al. (2006) and used to set an upper limit on the level of external shear.
The reconstructed observables corresponding to the best model for each galaxy (with the exception of SDSS J2321, for which we refer to the plots presented in C08) are shown and compared to the data in Figs 2 to 10. The reconstruction of both lensing and kinematic quantities appears in general to be very accurate. The residuals in the reconstructed lensed image are typically consistent with the noise level, and there is no sign of substantial discrepancies: this indicates that the underlying total density distribution is not significantly more complex (e.g. strongly triaxial) than the adopted axisymmetric powerlaw model, as discussed in Barnabè et al. (2008). The surface brightness distribution is also well reconstructed. The lowlevel ripples which are sometimes visible in the residuals map are due to the discrete nature of the TICs, and can be easily remedied by increasing the number of the employed TIC components, at the cost of a considerable slowdown of the optimization process, and without changing significantly the results. In the case of SDSS J0037, SDSS J0216 and SDSS J0959, the residuals clearly reveal the presence of the lensed image, which is not prominent in the surface brightness data map. The situation is more complicated for the kinematic maps, where the noise level is higher and the data sets composed by a smaller number of pixels (with as few as ten pixels for SDSS J0912). The models are fairly successful in reproducing the observables, in particular the velocity maps, with the exception of the velocity dispersion maps of SDSS J1627 (the reconstructed profile declines too rapidly) and possibly SDSS J0959 (there appears to be a stronger gradient in the data than in the model). Such difficulties in reproducing the velocity dispersion maps of these two systems might well reflect the shortcomings of the twointegral DF assumption (which implies isotropic velocity dispersion tensor in the meridional plane, i.e. ), whereas a more sophisticated threeintegral dynamical modelling could be required. We discuss this issue further in Section 6.
4.2 Mass distribution and dark matter fraction
For each galaxy, we have calculated from the best model the spherically averaged total mass profile, shown in Fig. 11. In order to assess the dark matter content of the analyzed systems, one also requires the corresponding profile for the luminous component, which can be obtained from the bestmodel reconstructed stellar DF. However, within the framework of the method, stars are just tracers of the total potential, and therefore an additional assumption is needed in order to constrain normalization for the luminous mass distribution (which is arbitrary unless the stellar masstolight ratio can be determined independently). Following C08, in this work we adopt the socalled maximum bulge approach, which consists of maximizing the contribution of the luminous component. In other words, the stellar density profile is maximally rescaled without exceeding the total density profile^{4}^{4}4This approach is effectively the earlytype galaxies equivalent of the classical maximum disk hypothesis (van Albada & Sancisi, 1986) frequently used in modelling the rotation curves of latetype galaxies., positing a positionindependent stellar masstolight ratio. In real galaxies the stellar masstolight ratio might not be uniform, although, based on observed colour gradients, the effect is not expected to be strong (e.g. Kronawitter et al., 2000). This procedure provides a uniform way to determine a lower limit for the dark matter fraction in the sample galaxies. Moreover, as shown in Barnabè et al. (2008), the maximum bulge approach is quite robust and allows a reliable determination of the dark matter fraction (within approximately 10 per cent of the total mass) even when the model assumptions are violated. The spherically averaged stellar mass profiles obtained under the maximum bulge hypothesis are also presented in Fig. 11, where we also indicate the threedimensional radii corresponding to the effective radius , the Einstein radius and the outermost boundary of the kinematic data . The boundary of the surface brightness map is at least comparable to , and often larger^{5}^{5}5The surface brightness maps considered for the combined analysis and presented in Figs 2, 4, 6, 8 and 10 are actually cutouts of HST images extending over 10 arcsec.. Although the spatial coverage of the lensing and kinematic data, in most cases, does not extend up to (the only exception being SDSS J0959, for which ), it should be noted that the more distant regions of the galaxy which are situated along the line of sight—and therefore observed in projection—also contribute fairly significantly in constraining the mass model, as extensively discussed in C08.
We find that the total mass profile closely follows the light in the very inner regions, which are presumably dominated by the stellar component, while dark matter typically starts playing a role in the vicinity of the (threedimensional) radius , where its contribution in total mass is of order to per cent, and becomes progressively more important when moving outwards (the system SDSS J0216, however, constitutes an exception, with its dark matter fraction remaining roughly constant over the probed region for kpc). Within a sphere of radius , approximately to per cent of the mass is dark. This result is in general agreement with the conclusions of previous dynamical studies of earlytype galaxies in the local Universe, in particular the analysis of nearly round and slowlyrotating ellipticals by Gerhard et al. (2001), the modelling of SAURON systems (under the assumption that mass follows light, Cappellari et al. 2006), and the study by Thomas et al. (2007) of earlytype galaxies in the Coma cluster.
From the value of the luminous mass inside the effective radius, obtained under the maximum bulge hypothesis, we also calculate for each system the corresponding upper limit for the stellar masstolight ratio (see Table 3), finding . This is in agreement with stellar population studies, e.g. Trujillo et al. (2004).
4.3 Error analysis
In this Section we present, for each galaxy in the sample, the corresponding model uncertainties, i.e. the errors on the recovered nonlinear parameters , , and . The uncertainties are calculated within the framework of Bayesian statistics by making use of the recently developed nested sampling technique (Skilling 2004, Sivia & Skilling 2006; see also Vegetti & Koopmans 2008 for the first astrophysical application in the context of gravitational lensing). Nested sampling is a Monte Carlo method aimed at calculating the Bayesian evidence, i.e. the fundamental quantity for model comparison, in a computationally efficient way. The marginalized posterior probability distribution functions (PDFs) of the model parameters, which are used to estimate the uncertainties, are obtained as very valuable byproducts of the method.
Within the context of Bayesian statistics, a priori assumptions or knowledge on each parameter are made explicit and formalized by defining the prior function . We assign a uniform prior within the interval , symmetrical around the recovered best model value and wide enough to include the bulk of the likelihood (very conservative estimates of are obtained by means of fast preliminary runs), that is:
(6) 
This choice of an uniform prior is aimed at formalizing the absence of any a priori information within the interval (see e.g. Cousins, 1995). We find, however, that the errors on the parameters are very small in comparison with , so that the prior, largely independently of the adopted functional form, is nearly constant over the likelihood. Therefore, the specific choice for is not critical in our case.
For each analyzed galaxy, the histograms in Fig. 12 show the marginalized posterior PDF of the powerlaw model nonlinear parameters. Because of the marginalization involved in their evaluation, these distributions constitute the most conservative estimate of statistical errors on the parameters, given the data and all the assumptions (i.e. positing that the adopted model is the “true” description underlying the data; cf. MacKay 1992). These errors are relatively small, as a consequence of the numerous constraints provided by the data: the typical data set for each of the sample galaxies consists of data points or more, most of them (in the lensed image and surface brightness maps) characterized by fairly good signal to noise level. Furthermore, the maximum of the posterior PDF (and, more generally, the bulk of the posterior probability) for the th parameter is often found to be somewhat skewed with respect to the corresponding best model value . This is a wellknown projection effect arising from the marginalization over a single parameter of a complicated highdimensional multivariate function such as the total posterior PDF, as discussed above (§ 3.1).
The analysis conducted in this Section does not take into account systematic uncertainties, which are frequently larger than the statistical errors, and more difficult to quantify. They can arise from a variety of sources, including incorrect modelling assumptions and problems associated with the generation of the data sets (see e.g. Marshall et al. 2007 for an indepth treatment of the systematic uncertainties connected with the lens galaxy subtraction process and the incomplete knowledge of the PSF). The study of Barnabè et al. (2008) provides a more quantitative feel for the systematic errors introduced by the adoption of an oversimplified galaxy model, showing that even in a quite extreme case (where the reconstruction of lensing observables is clearly unsatisfactory) the systematic error on remains within about per cent. Other parameters such as inclination and oblateness, however, are less robust and actually become illdefined if the assumption of axial symmetry does not hold. We note, however, that in none of the six systems under study are the model residuals as severe as in the simulations in Barnabè et al. (2008). Therefore, we expect systematic uncertainties to remain within a few per cent level.
4.4 The density profile of the ensemble
From the combined analysis, we have found that all the galaxies in the ensemble have a total density profile very close to isothermal, with an average logarithmic slope , in agreement with the results of Koopmans et al. (2006). There is also no evidence of evolution of the density slope with redshift (see Fig. 13).
No correlation is found between the logarithmic slope and the axial ratio , the effective radius, the normalized Einstein radius (i.e., the ratio ) and the aperture averaged velocity dispersion .
We now want to calculate, on the basis of the analyzed systems, the intrinsic spread around the average slope. If we assume that the slope of each galaxy is extracted from a parent Gaussian distribution of centre and variance , then the joint posterior probability for the sampling is given by
(7) 
where is the prior on and (for which we adopt a uniform distribution) and are the errors on , calculated by considering the region around the peak which contains 68 per cent of the posterior probability. From Eq. (7), the maximum likelihood solution for is simply the average slope , while for is obtained from the equation
(8) 
Inserting the values for and determined in our analysis, we solve Eq. (8) to find, for our sample of six distant ellipticals, an intrinsic spread around the average logarithmic slope, corresponding to per cent. A measure of the joint posterior of these quantities is provided by Fig. 14, where we plot and we draw contours corresponding to posterior (or likelihood in the case of flat prior) ratios , with . We note that these contours are only for indication, and formally have a proper meaning only in the case of a multivariate Gaussian, in which case is the usual chisquare.
If we consider a different prior in Eq. (7), the outcome is only slightly modified. For instance, if we adopt , which formalizes the absence of a priori information on the scale (i.e. the order of magnitude) of , we find a maximum likelihood value of . This points out that, despite the fact that we only have a handful of systems, the results are essentially driven by the data, with the choice of the prior playing only a minimal role.
4.5 Axial ratio of the density distribution
The axial ratio of the total density distribution is found to be always rounder than . It is interesting to compare this quantity with the intrinsic axial ratio of the luminous distribution, obtained by deprojecting the observed isophotal axial ratio , i.e.
(9) 
where we use the best model value for the inclination . The results are illustrated in Fig. 15. For four of the galaxies in the sample, the flattenings coincide closely, while the total distribution is rounder in the case of SDSS J1627 and SDSS J0959. The discrepancy is particularly conspicuous for SDSS J0959, where , while it is only for SDSS J1627. Intriguingly, SDSS J0959 is the only clearly fastrotating galaxy in the sample (see the velocity map in Fig. 8 and the strongly asymmetric DF in Fig. 16(d)), and has also peculiar dynamical properties when compared with the rest of the ensemble, as discussed in Section 5.
5 Recovered dynamical structure
As reviewed in Section 3, for any given total gravitational potential the cauldron algorithm determines the bestfitting dynamical model by means of TICs superposition. Therefore, the best model of each galaxy has an associated best reconstructed map of the relative TICs weights, which is a representation in the integral space of the corresponding (weighted) twointegral DF. The maps are shown in Fig. 16. In this work, consistently with C08, we employ a library of TICs for the dynamical modelling. The TIC grid is constructed by considering elements logarithmically sampled in the circular radius , and, for each of them, elements linearly sampled in angular momentum between and , mirrored in the negative values^{6}^{6}6The grid in the radial coordinate corresponds to a grid in the energy , where (10) and the circular speed is given by: (11) The choice for the circular radius also sets the corresponding maximum angular momentum along the axis: . . Even though this grid may appear coarse due to the limited number of TICs employed, it represents—for the applications described in this paper—an excellent compromise between computational efficiency and quality of the observables reconstruction. Increasing the number of grid elements has the effect of improving the surface brightness reconstruction, but does not change significantly the values of the recovered parameters. It also makes the reconstructed TIC weights map smoother, while preserving the main features already visible in the corresponding coarse map.
While the phasespace DF completely determines the structure and dynamics of a (collisionless) stellar system, it is not immediately intuitive to interpret. Therefore, from the best reconstructed DFs we derive dynamical characterizations of the galaxies, such as the anisotropy parameters (§ 5.1) the diagram (§ 5.2), and the angular momentum (§ 5.3) which can be more easily related with the observations.
5.1 Anisotropy parameters
The global anisotropy of an elliptical galaxy is related to the distribution of its stellar orbits and is often considered an important indicator of the assembly mechanism of the system (see e.g. Burkert & Naab, 2005).
For an axisymmetric system, the global shape of the velocity dispersion tensor can be quantified by using the three anisotropy parameters (Cappellari et al., 2007; Binney & Tremaine, 2008)
(12) 
(13) 
and
(14) 
where
(15) 
denotes the total unordered kinetic energy in the coordinate direction , and is the velocity dispersion along the direction at any given location in the galaxy. For an isotropic system (e.g. the classic isotropic rotator) the three parameters are all zero. Stellar systems supported by a twointegral DF, as assumed in our case, are isotropic in the meridional plane, i.e. everywhere, which implies .
For each object in the sample, we have calculated the anisotropy parameters by integrating up to half the effective radius, which is the typical region inside which the galaxy models are more strongly constrained. The results are reported in Table 4. All the systems, with the exception of SDSS J0959, have slightly positive , i.e. are mildly anisotropic in the sense of having larger pressure parallel to the equatorial plane than perpendicular to it. This is quite similar to what Cappellari et al. (2007) and Thomas et al. (2008), using threeintegral axisymmetric orbitsuperposition codes, find for local ellipticals; however, their samples also display a few galaxies with clearly higher anisotropy (). The fastrotating galaxy SDSS J0959, instead, is anisotropic in the opposite sense, due to the fact that for this system over most of the densityweighted volume, which translates into a negative parameter. This property is uncommon, although not unprecedented, for nearby earlytype galaxies: two systems out of the analyzed by Thomas et al. (2008) have , while no case is reported from the SAURON sample.
Since for our models by construction, the two remaining anisotropy parameters are univocally related by Eq. 14 so that, for , is necessarily negative. The twointegral DF assumption in general does not hold for nearby ellipticals (e.g. Merrifield, 1991; Gerssen et al., 1997; Thomas et al., 2008); if this is the case also for distant ellipticals, then the recovered could be significantly in error. On the other hand, as shown in Barnabè et al. (2008), the global parameter is more robust, and can be reliably recovered by cauldron (typically within per cent) even when the assumptions of twointegral DF and axial symmetry are both violated.
5.2 The global and local
The diagram provides a classic indicator to estimate the importance of rotation with respect to random motions in earlytype galaxies (see Binney, 1978). For each system, we calculate the “intrinsic” (i.e. inclination corrected) global quantity from the best reconstructed DF by integrating up to (the results are listed in Table 4), and we plot it against the intrinsic ellipticity of the luminous distribution . The diagram is presented in Fig. 17 and compared with the findings for SAURON galaxies, corrected for inclination (Cappellari et al., 2007). There is a sharp dichotomy in the SLACS subsample between the obviously fastrotating system SDSS J0959 and the remaining galaxies, four of which have (two of these systems have clear characteristics of slow rotators, as discussed in Section 5.3).
Whereas is a global quantity which provides information on the general properties of the galaxy, important insights on the characteristics of the system at different locations of the meridional plane are offered by its local analogue, i.e. the ratio between the mean rotation velocity around the axis and the mean velocity dispersion . We illustrate this quantity in Fig. 18 for the galaxies in our sample. For visualization purposes, the plot was produced from a weighted DF map of elements ( and , mirrored). In order to do this, we reoptimized the best model for the dynamics hyperparameters only, determining the new optimal level of the regularization, while all the other parameters were kept fixed. This procedure amounts to some extent to an interpolation over the weighted DF map of Fig. 16 with the aim of obtaining a smoother distribution.
5.3 Angular momentum
Another robust way to characterize the global velocity structure of a galaxy is provided by its angular momentum content. For each galaxy in the ensemble we calculate the (massnormalized) component of the angular momentum parallel to the axis of symmetry as
(16) 
where is the radial coordinate, denotes the mean azimuthal stellar velocity at position and is the spatial density of stars as obtained by the best reconstructed DF, i.e. . The results obtained by integrating inside (the region most strongly constrained by the observations for all the systems) are reported in Table 4 in units of kpc km s.
Whereas the dimensional parameter has a direct physical interpretation as angular momentum, it is not the most practical way to quantify and to compare the level of ordered rotation in elliptical galaxies. To this purpose, we define a more convenient adimensional parameter as:
(17) 
This quantity is effectively the intrinsic equivalent of the observationallydefined parameter introduced by Emsellem et al. (2007) as an objective criterion for the kinematic classification of earlytype galaxies. Analogously to , tends to unity for systems which display largescale ordered rotation, and conversely it goes to zero for galaxies globally dominated by random motions, whereas the same galaxies might have a moderate ratio due to the presence of smallscale rotation patterns in the highdensity central regions.
We have computed for the SLACS subsample by integrating inside half the effective radius, listing the results in Table 4. The small number of galaxies in the sample and the limited spatial coverage and quality of the kinematic data sets do not allow us to trace a sharp demarcation line between slow and fast rotators; nevertheless, there is a clear indication that SDSS J0959 belongs to the latter, while SDSS J0216 and SDSS J2321 are part of the first group, although they could not be straightforwardly identified as slow rotators on the basis of the diagram only. The remaining three galaxies lie somewhat in between.
Galaxy  

J0037  0.16  0.37  112  0.248  
J0216  0.08  0.22  11  0.116  
J0912  0.07  0.24  0.229  
J0959  0.27  0.76  0.645  
J1627  0.16  0.23  0.181  
J2321  0.14  0.20  0.075 
6 Summary and conclusions
We have conducted, for the first time, an indepth analysis of the mass distribution and dynamical structure of a sample of massive earlytype galaxies beyond the very local Universe, with a redshift range of .
The examined systems are six earlytype lens galaxies from the SLACS survey for which both HST/ACS highresolution imaging and VLT VIMOS integral field spectroscopy are available. These unique, highquality data sets of earlytype galaxies beyond the local Universe have enabled us to carry out a joint analysis of these systems, by combining gravitational lensing and stellar dynamics in a fully selfconsistent way (using the specifically designed code cauldron). The method is completely embedded within the framework of Bayesian statistics, permitting an objective datadriven determination of the “best model”, given the observations and our priors (formalized by the choice of the regularization). This technique makes it possible—under the assumptions of axial symmetry and twointegral stellar DF—to disentangle to a large extent several classical degeneracies and to effectively “dissect” the investigated galaxies, recovering their intrinsic structure. We summarize and discuss as follows the main results of this study:

The global density distribution of massive earlytype galaxies within approximately is remarkably well described by a simple axisymmetric powerlaw profile. Despite being very sensitive to the features of the underlying mass distribution (as shown by simulations of non axially symmetric systems, e.g. Barnabè et al. 2008; see also Koopmans 2005 and Vegetti & Koopmans 2008), the lensed images can be reconstructed almost to the noise level by adopting a model (not even a weak external shear is required), indicating a surprising degree of smoothness in the mass structure of ellipticals. While this conclusion could already be envisioned from the results of the Koopmans et al. (2006) study of SLACS lenses, here we have shown that such smooth models are also consistent with the observed surface brightness and kinematics maps. We suggest that this significantly smooth structure might be related to the formation mechanisms of earlytype galaxies.

The average logarithmic slope of the total mass density distribution is , with an intrinsic spread of per cent. The galaxies in the sample have therefore a density profile consistent with isothermal, corresponding to flat rotation curves, inside a range of Einstein radii of . This is in agreement with the findings of previous studies of ellipticals both in the local Universe and up to redshift of (e.g. Gerhard et al. 2001, Koopmans et al. 2006, Thomas et al. 2007).

There is no evidence for evolution of the logarithmic total density slope within the probed range of redshifts, although this is not surprising given the findings of the (non selfconsistent) combined lensing and dynamics study of Koopmans et al. (2006). However, it does show that we have systematics well under control.

The shape of the total density distribution is fairly round, with an axial ratio , and does not differ much from the intrinsic axial ratio of the luminous distribution (obtained via deprojection of the observed isophotal ratio, by making use of the recovered best model value for the inclination). The only exception is represented by the case of the lenticular galaxy SDSS J0959, which has a total density profile much rounder that the luminous one.

The lower limit for the dark matter fraction, calculated within spherical shells under the hypotheses of “maximum bulge” and positionindependent stellar masstolight ratio, falls in the range per cent (of the total mass) at half the effective radius, and rises to per cent at . This is fully consistent with the results of several studies of the dark and luminous mass distribution in nearby ellipticals (in particular Gerhard et al., 2001; Cappellari et al., 2006; Thomas et al., 2007).

The SLACS galaxies in our subsample are only mildly anisotropic, with the global parameter . Five out of six systems are slightly flattened (beyond the contribution of rotation) by having a larger pressure support parallel to the equatorial plane rather than perpendicular to it; the situation is reversed in the case of SDSS J0959, which shows negative .

From the inspection of the stellar velocity maps, the global and local ratio and, more decisively, the intrinsic rotation parameter (directly related to the angular momentum of the stellar component of the galaxy), it is possible to objectively quantify the level of ordered rotation with respect to the random motions. Two of the systems, namely SDSS J0216 and SDSS J2321, are identified as slow rotators, while SDSS J0959 unambiguously shows the characteristics of fast rotators, including the relatively low mass and luminosity: with it is the only system in our sample to lie within the typical range of absolute magnitudes for fast rotators determined by Emsellem et al. (2007). The remaining three galaxies exhibit a moderate amount of rotation.
Overall, the earlytype lens galaxies analyzed in this paper, in the redshift range , are found to be very similar to the local massive ellipticals, in terms of structural and dynamical properties of their inner regions. Our study shows, for the first time, that also the physical distinction between slow and fast rotators (originally revealed by the SAURON survey: Emsellem et al. 2007; Cappellari et al. 2007) is already in place at redshift , although a larger sample is necessary in order to quantify this more precisely to higher redshifts. Since a series of studies (Treu et al., 2006; Bolton et al., 2008; Treu et al., 2008) has shown that the SLACS systems are statistically identical—in terms of observational properties and environment—to nonlens galaxies of comparable size and luminosity, we can generalize the results of this work and conclude that at least the most massive elliptical galaxies did not experience any major evolutionary process in their global structural properties between redshift and . Since the lookback time is only Gyr, this might not be surprising. However, pushing lensing and dynamics techniques back in redshift and cosmic time is crucial if we ever wish to fully understand the structural evolution of earlytype galaxies. In this paper a first step has been taken, using more sophisticated selfconsistent techniques. This goes beyond what is possible with the use of lensing or dynamics alone.
On the other hand, we also find differences between the analyzed systems and nearby galaxies when we compare the respective global anisotropy parameters. Although the values of recovered for the SLACS subsample fall within the typical range for local galaxies, we note that we do not find any systems with , which are instead quite common in the SAURON and Coma samples. This might be due, however, to the very modest size of our sample, particularly in terms of fastrotating objects. A much more drastic discrepancy is evident in the distributions of the anisotropy parameters and , plausibly due to the limitations of our assumption of twointegral DF, which imposes at every location. Isotropy in the meridional plane is not observed, in general, for local ellipticals, where usually (e.g. Gerssen et al., 1997; Thomas et al., 2008). If this applies also to more distant systems, then the anisotropy will not be correctly estimated by our method, and more sophisticated axisymmetric threeintegral models might provide a better dynamical description of the galaxy. We foresee future developments of the cauldron algorithm in this direction, concurrently with the availability of improved kinematic data sets. As for the present analysis, however, the assumption generally appears to work satisfactorily, providing generally correct reconstructions of the observables. Furthermore, as tested in Barnabè et al. (2008), the current method is robust enough to recover in a reliable way several important global properties of the analyzed galaxy (such as the density slope, dark matter fraction, angular momentum and parameter, although not the flattening and the and parameters) even if it is applied to a complex, nonsymmetric system which departs significantly from the idealized assumptions of axisymmetry and two or threeintegral DF, and which is likely more extreme than the typical ellipticals under study.
In future papers in this series, we plan to extend the current study by applying the combined analysis to all the systems for which twodimensional kinematic maps are or will become available. This includes the entire sample of SLACS earlytype lens galaxies for which integral field spectroscopy has been obtained. Twodimensional kinematic maps can also be obtained for a further lens galaxies for which longslit spectroscopic observations have been conducted with the instrument LRIS mounted on KeckI, with the slit positions aligned with the major axis of the system and offset along the minor axis in order to mimic integral field spectroscopy.
Acknowledgments
M. B. is grateful to Mattia Righi for many lively and stimulating discussions. M. B. acknowledges the support from an NWO program subsidy (project number 614.000.417). OC and LVEK are supported (in part) through an NWOVIDI program subsidy (project number 639.042.505). We also acknowledge the continuing support by the European Community’s Sixth Framework Marie Curie Research Training Network Programme, Contract No. MRTNCT2004505183 ‘ANGLES’. TT acknowledges support from the NSF through CAREER award NSF0642621, by the Sloan Foundation through a Sloan Research Fellowship, and by the Packard Foundation through a Packard Fellowship. Support for programs #10174 and #10494 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 526555.
References
 Arnaboldi et al. (1996) Arnaboldi M., Freeman K. C., Mendez R. H., Cappaccioli M., Ciardullo R., Ford H., Gerhard O., Hui X., Jacoby G. H., Kudritzki R. P., Quinn P. J., 1996, ApJ, 472, 145
 Barnabè & Koopmans (2007) Barnabè M., Koopmans L. V. E., 2007, ApJ, 666, 726
 Barnabè et al. (2008) Barnabè M., Nipoti C., Koopmans L. V. E., Vegetti S., Ciotti L., 2008, MNRAS, in press
 Barnes (1992) Barnes J. E., 1992, ApJ, 393, 484
 Barth et al. (2002) Barth A. J., Ho L. C., Sargent W. L. W., 2002, AJ, 124, 2607
 Bertin (2008) Bertin E., 2008, SWARP User’s Guide. v2.17.0 edn
 Bertin et al. (1994) Bertin G., Bertola F., Buson L. M., Danzinger I. J., Dejonghe H., Sadler E. M., Saglia R. P., de Zeeuw P. T., Zeilinger W. W., 1994, A&A, 292, 381
 Binney (1978) Binney J., 1978, MNRAS, 183, 501
 Binney (2005) Binney J., 2005, MNRAS, 363, 937
 Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics. Princeton University Press
 Bolton et al. (2008) Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Gavazzi R., Moustakas L. A., Wayth R., Schlegel D. J., 2008, ApJ, 682, 964
 Bolton et al. (2006) Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Moustakas L. A., 2006, ApJ, 638, 703
 Bolton et al. (2008) Bolton A. S., Treu T., Koopmans L. V. E., Gavazzi R., Moustakas L. A., Burles S., Schlegel D. J., Wayth R., 2008, ApJ, 684, 248
 Borriello et al. (2003) Borriello A., Salucci P., Danese L., 2003, MNRAS, 341, 1109
 Burkert & Naab (2005) Burkert A., Naab T., 2005, MNRAS, 363, 597
 Cappellari et al. (2006) Cappellari M., Bacon R., Bureau M., Damen M. C., Davies R. L., de Zeeuw P. T., Emsellem E., FalcónBarroso J., et al., 2006, MNRAS, 366, 1126
 Cappellari & Emsellem (2004) Cappellari M., Emsellem E., 2004, PASP, 116, 138
 Cappellari et al. (2007) Cappellari M., Emsellem E., Bacon R., Bureau M., Davies R. L., de Zeeuw P. T., FalcónBarroso J., Krajnović D., Kuntschner H., McDermid R. M., Peletier R. F., Sarzi M., van den Bosch R. C. E., van de Ven G., 2007, MNRAS, 379, 418
 Carollo et al. (1995) Carollo C. M., de Zeeuw P. T., van der Marel R. P., Danziger I. J., Qian E. E., 1995, ApJ, 441, L25
 Chandrasekhar (1969) Chandrasekhar S., 1969, Ellipsoidal Figures of Equilibrium. Yale University Press
 Cole et al. (2000) Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168
 Côté et al. (2003) Côté P., McLaughlin D. E., Cohen J. G., Blakeslee J. P., 2003, ApJ, 591, 850
 Cousins (1995) Cousins R. D., 1995, American Journal of Physics, 63, 398
 Cretton et al. (1999) Cretton N., de Zeeuw P. T., van der Marel R. P., Rix H.W., 1999, ApJS, 124, 383
 Czoske et al. (2008) Czoske O., Barnabè M., Koopmans L. V. E., 2008, in Probing Stellar Populations out to the Distant Universe Conference, (arXiv:0811.2391) Integralfield spectroscopy of SLACS lenses
 Czoske et al. (2008) Czoske O., Barnabè M., Koopmans L. V. E., Treu T., Bolton A. S., 2008, MNRAS, 384, 987
 de Lorenzi et al. (2008) de Lorenzi F., Gerhard O., Saglia R. P., Sambhus N., Debattista V. P., Pannella M., Méndez R. H., 2008, MNRAS, 385, 1729
 de Zeeuw et al. (2002) de Zeeuw P. T., Bureau M., Emsellem E., Bacon R., Carollo C. M., Copin Y., Davies R. L., Kuntschner H., et al., 2002, MNRAS, 329, 513
 Emsellem et al. (2007) Emsellem E., Cappellari M., Krajnović D., van de Ven G., Bacon R., Bureau M., Davies R. L., de Zeeuw P. T., et al., 2007, MNRAS, 379, 401
 Emsellem et al. (2004) Emsellem E., Cappellari M., Peletier R. F., McDermid R. M., Bacon R., Bureau M., Copin Y., Davies R. L., Krajnović D., Kuntschner H., Miller B. W., de Zeeuw P. T., 2004, MNRAS, 352, 721
 Fabbiano (1989) Fabbiano G., 1989, ARAA, 27, 87
 Falco et al. (1985) Falco E. E., Gorenstein M. V., Shapiro I. I., 1985, ApJ, 289, L1
 Franx et al. (1994) Franx M., van Gorkom J. H., de Zeeuw T., 1994, ApJ, 436, 642
 Fukazawa et al. (2006) Fukazawa Y., BotoyaNonesa J. G., Pu J., Ohto A., Kawano N., 2006, ApJ, 636, 698
 Gavazzi et al. (2008) Gavazzi R., Treu T., Koopmans L. V. E., Bolton A. S., Moustakas L. A., Burles S., Marshall P. J., 2008, ApJ, 677, 1046
 Gavazzi et al. (2007) Gavazzi R., Treu T., Rhodes J. D., Koopmans L. V. E., Bolton A. S., Burles S., Massey R. J., Moustakas L. A., 2007, ApJ, 667, 176
 Gerhard et al. (2001) Gerhard O., Kronawitter A., Saglia R. P., Bender R., 2001, AJ, 121, 1936
 Gerhard (1993) Gerhard O. E., 1993, MNRAS, 265, 213
 Gerssen et al. (1997) Gerssen J., Kuijken K., Merrifield M. R., 1997, MNRAS, 288, 618
 Humphrey et al. (2006) Humphrey P. J., Buote D. A., Gastaldello F., Zappacosta L., Bullock J. S., Brighenti F., Mathews W. G., 2006, ApJ, 646, 899
 Jesseit et al. (2007) Jesseit R., Naab T., Peletier R. F., Burkert A., 2007, MNRAS, 376, 997
 Kelson et al. (2000) Kelson D. D., Illingworth G. D., van Dokkum P. G., Franx M., 2000, ApJ, 531, 159
 Kochanek (1991) Kochanek C. S., 1991, ApJ, 373, 354
 Koopmans et al. (2006) Koopmans L., Treu T., Bolton A. S., Burles S., Moustakas L. A., 2006, ApJ, 649, 599
 Koopmans & Treu (2002) Koopmans L. V. E., Treu T., 2002, ApJ, 568, L5
 Koopmans & Treu (2003) Koopmans L. V. E., Treu T., 2003, ApJ, 583, 606
 Koopmans (2005) Koopmans Léon V. E., 2005, MNRAS, 363, 1136
 Kronawitter et al. (2000) Kronawitter A., Saglia R. P., Gerhard O., Bender R., 2000, A&AS, 144, 53
 Loewenstein & White (1999) Loewenstein M., White III R. E., 1999, ApJ, 518, 50
 MacKay (1992) MacKay D. J. C., 1992, PhD Thesis
 MacKay (1999) MacKay D. J. C., 1999, Neural Comp, 11, 1035
 MacKay (2003) MacKay D. J. C., 2003, Information Theory, Inference and Learning Algorithms. Cambridge University Press
 Marshall et al. (2007) Marshall P. J., Treu T., Melbourne J., Gavazzi R., Bundy K., Ammons S. M., Bolton A. S., Burles S., Larkin J. E., Le Mignant D., Koo D. C., Koopmans L. V. E., Max C. E., Moustakas L. A., Steinbring E., Wright S. A., 2007, ApJ, 671, 1196
 Matsushita et al. (1998) Matsushita K., Makishima K., Ikebe Y., Rokutanda E., Yamasaki N., Ohashi T., 1998, ApJ, 499, L13
 Merrifield (1991) Merrifield M. R., 1991, AJ, 102, 1335
 Mould et al. (1990) Mould J. R., Oke J. B., de Zeeuw P. T., Nemec J. M., 1990, AJ, 99, 1823
 Naab et al. (2006) Naab T., Khochfar S., Burkert A., 2006, ApJ, 636, L81
 Rix et al. (1997) Rix H.W., de Zeeuw P. T., Cretton N., van der Marel R. P., Carollo C. M., 1997, ApJ, 488, 702
 Romanowsky et al. (2003) Romanowsky A. J., Douglas N. G., Arnaboldi M., Kuijken K., Merrifield M. R., Napolitano N. R., Capaccioli M., Freeman K. C., 2003, Science, 301, 1696
 Saglia et al. (1992) Saglia R. P., Bertin G., Stiavelli M., 1992, ApJ, 384, 433
 Schneider et al. (1992) Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses. Berlin, SpringerVerlag
 Schwarzschild (1979) Schwarzschild M., 1979, ApJ, 232, 236
 Scodeggio et al. (2005) Scodeggio M., Franzetti P., Garilli B., Zanichelli A., Paltani S., Maccagni D., Bottini D., Le Brun V., et al., 2005, PASP, 117, 1284
 Sivia & Skilling (2006) Sivia D. S., Skilling J., 2006, Data Analysis: A Bayesian Tutorial. Oxford University Press
 Skilling (2004) Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, American Institute of Physics Conference Series Vol. 735 of American Institute of Physics Conference Series, Nested Sampling. pp 395–405
 Suyu et al. (2006) Suyu S. H., Marshall P. J., Hobson M. P., Blandford R. D., 2006, MNRAS, 371, 983
 Thomas et al. (2007) Thomas J., Jesseit R., Naab T., Saglia R. P., Burkert A., Bender R., 2007, MNRAS, 381, 1672
 Thomas et al. (2008) Thomas J., Jesseit R., Saglia R. P., Bender R., Burkert A., Corsini E. M., Gebhardt K., Magorrian J., Naab T., Thomas D., Wegner G., 2008, MNRAS, in press
 Thomas et al. (2007) Thomas J., Saglia R. P., Bender R., Thomas D., Gebhardt K., Magorrian J., Corsini E. M., Wegner G., 2007, MNRAS, 382, 657
 Toomre (1977) Toomre A., 1977, in Tinsley B. M., Larson R. B., eds, Evolution of Galaxies and Stellar Populations Mergers and Some Consequences. p. 401
 Treu et al. (2008) Treu T., Gavazzi R., Gorecki A., Marshall P. J., Koopmans L. V. E., Bolton A. S., Moustakas L. A., Burles S., 2008, ApJ, in press
 Treu & Koopmans (2002) Treu T., Koopmans L. V. E., 2002, MNRAS, 337, L6
 Treu & Koopmans (2003) Treu T., Koopmans L. V. E., 2003, MNRAS, 343, L29
 Treu & Koopmans (2004) Treu T., Koopmans L. V. E., 2004, ApJ, 611, 739
 Treu et al. (2006) Treu T., Koopmans L. V. E., Bolton A. S., Burles S., Moustakas L., 2006, ApJ, 640, 662
 Trott et al. (2008) Trott C. M., Treu T., Koopmans L. V. E., Webster R. L., 2008, preprint (arXiv:0812.0748)
 Trujillo et al. (2004) Trujillo I., Burkert A., Bell E. F., 2004, ApJ, 600, L39
 Valdes et al. (2004) Valdes F., Gupta R., Rose J. A., Singh H. P., Bell D. J., 2004, ApJS, 152, 251
 van Albada & Sancisi (1986) van Albada T. S., Sancisi R., 1986, Royal Society of London Philosophical Transactions Series A, 320, 447
 van de Ven et al. (2008) van de Ven G., FalconBarroso J., McDermid R. M., Cappellari M., Miller B. W., de Zeeuw P. T., 2008, preprint (arXiv:0807.4175)
 van der Marel & Franx (1993) van der Marel R. P., Franx M., 1993, ApJ, 407, 525
 Vegetti & Koopmans (2008) Vegetti S., Koopmans L. V. E., 2008, MNRAS, in press
 Verolme & de Zeeuw (2002) Verolme E. K., de Zeeuw P. T., 2002, MNRAS, 331, 959
 White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
 Wucknitz (2002) Wucknitz O., 2002, MNRAS, 332, 951
 Zanichelli et al. (2005) Zanichelli A., Garilli B., Scodeggio M., Franzetti P., Rizzo D., Maccagni D., Merighi R., Picat J.P., et al., 2005, PASP, 117, 1271