EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 643 (November 2020) A&A, 643 (2020) A135 Full HTML
Open Access
Issue
A&A
Volume 643, November 2020
Article Number A135
Number of page(s) 16
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202037929
Published online 16 November 2020

© G. Chirivì et al. 2020

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1. Introduction

The progenitors of today’s massive galaxy population are thought to be small and dense (Daddi et al. 2005; Trujillo et al. 2006; Zirm et al. 2007; van der Wel et al. 2008, 2014; van Dokkum et al. 2008; Szomoru et al. 2010, 2012), discy (Toft et al. 2005; Trujillo et al. 2006; van der Wel et al. 2011; Chang et al. 2013), with quenched star formation and old stellar populations (Kriek et al. 2006, 2008, 2009; Toft et al. 2007; Cimatti et al. 2008; van Dokkum et al. 2010), small Sérsic indices, and high stellar velocity dispersion (van Dokkum et al. 2009; Bezanson et al. 2011; Toft et al. 2012; van de Sande et al. 2013) (see Cappellari 2016, for a review). How they came to be nowadays’ most massive galaxies, enclosing most of the stellar mass in the Universe (Fukugita et al. 1998; Hogg et al. 2002; Bell et al. 2003; Baldry et al. 2004), is still a debated topic. Significant progress on the numerical side now suggests a two-phase scenario for the formation and evolution of the massive galaxy population (Oser et al. 2010; Wellons et al. 2016), while our comprehension on the observational side has developed immensely by means of integral field unit (IFU) observations. These have proven to be a groundbreaking tool to unveil the structural and dynamical properties of galaxies (see Cappellari 2016, for a review). In fact, integral field spectroscopy (IFS) has become an essential tool in astrophysics, allowing us to obtain a spectrum in every spaxel on the sky covering the entire galaxy field and, therefore, obtaining a 3D data cube as outcoming data (Bacon et al. 1995; de Zeeuw et al. 2002). From these data cubes, we can reconstruct the stellar and gas kinematic maps, which are important tracers of the underlying gravitational potential due to visible (stars and gas) and non-visible (dark) matter, allowing us to explore the matter content, matter distribution, and internal dynamics of galaxies, which hold clues regarding their assembly history. A way in which this information can be quantified is via dynamical modelling, which allows for a detailed description of the galaxy’s dynamics. Key results were obtained using dynamical models fitted to stellar and gas kinematics, including, for example, mass determinations of supermassive black holes in galaxies (e.g. Gebhardt et al. 2000; Barth et al. 2001; Cappellari et al. 2002; Sarzi et al. 2001; Gültekin et al. 2009a,b, 2012; van den Bosch & de Zeeuw 2010; McConnell et al. 2012; van den Bosch et al. 2012; Walsh et al. 2013; Thomas et al. 2016), the determination of their stellar mass-to-light ratios, dark matter fractions, total mass profiles, and slopes (e.g. Kronawitter et al. 2000; Weijmans et al. 2009; Murphy et al. 2011; Cappellari et al. 2013, 2015; Yıldırım et al. 2017). The revelation of fast and slow rotators (Emsellem et al. 2007) and how these are potentially linked to the formation and evolution histories of today’s most massive galaxies are also important findings (Naab et al. 2014). However, obtaining resolved gas and stellar kinematic maps is mainly possible for galaxies in the local Universe.

To unveil the properties of galaxies in the high-redshift Universe, we can exploit natural cosmic telescopes such as gravitational lenses. Gravitational lensing is a relativistic effect for which the light travelling from a source towards the observer is bent by the presence of matter (baryonic and dark) in between them. Consequently, the source will be observed at a different position than it actually is, distorted in shape, and in some cases also multiply imaged (in the so-called “strong lensing” regime). It will also appear magnified by a factor μ (magnification). By modelling the lens mass distribution, one can reconstruct the image positions and magnification μ of the source and, thanks to surface brightness conservation, the surface brightness distribution of the source galaxy can be reconstructed. This will allow one to study high-redshift galaxies, providing crucial probes of structure formation and galaxy evolution (e.g. Oldham et al. 2017, among others). Combined lensing and stellar dynamic techniques were employed in previous works (e.g. Treu & Koopmans 2002, 2004; Barnabè & Koopmans 2007; Koopmans et al. 2009; van de Ven et al. 2010; Barnabè et al. 2011, 2012, among others), mainly to study properties of lens galaxies.

In this paper we present GLaD (Gravitational Lensing and Dynamics), a software that is able to unify dynamical modelling of galaxies with the magnification power of strong gravitational lensing, to reconstruct the dynamical properties of high-redshift source galaxies. GLaD is able to model the source stellar kinematics using axisymmetric Jeans modelling, assuming a source surface brightness and surface mass density profile, and then predict how these maps will look like when lensed and distorted by a strong gravitational lens. This allows us to compare the predicted arcs and images directly to the observed strong gravitational lensing systems, in order to provide improved constraints for both the source and the deflector.

Previous works have applied similar methods exploiting the gravitational lensing magnification to study spatially resolved kinematics of background sources (e.g. Jones et al. 2010; Dye et al. 2015; Rybak et al. 2015; Swinbank et al. 2015; Newman et al. 2017, 2018; Di Teodoro et al. 2018; Girard et al. 2018; Patrício et al. 2018; Rizzo et al. 2018). However, these methods mostly rely on a pre-modelled lens model, which is kept fixed during the dynamical analysis. This technique is suboptimal for it does not quantify degeneracies between the lens mass and source kinematics. Another trait which is common to these works is performing the dynamical analysis on the source plane, by “de-lensing” the kinematic map from the lens plane to the source plane. A disadvantage of this approach is the induced pixel correlation derived from the “de-lensing”. Moreover, it is not clear how to properly characterise noise properties on the source plane. Finally, the resolution on the source plane is dependent on differential magnification, which is an effect that must be taken into account. Some recent works overcame these issues partially or totally (Patrício et al. 2018; Rizzo et al. 2018), but differ from our method for the use of different dynamical modelling techniques and scientific goal, since they focus on the star-forming population by tracing the gas kinematics, which are not necessarily a pure tracer of the gravitational potential. Moreover, and unlike previous methods, we do not rely on a pixellated source reconstruction but instead assume parametrised profiles for the source, from which we can easily recover source properties of interest such as its total mass, ellipticity, Sérsic index and effective size.

The paper is organised as follows: we describe our method for both the lensing and the dynamical analysis in Sect. 2. We present a test case based on mock data of RX J1131−1231 and show how well we are able to recover the lensing and dynamical parameters using GLaD in Sect. 3. We discuss and conclude in Sect. 4. Throughout the paper, parameter constraints are given by the median values with the uncertainties given by the 16th and 84th percentiles (corresponding to 68% credible intervals (CI)) of the marginalised probability density distributions. We assume a flat ΛCDM cosmology with H0 = 70 km s−1 Mpc−1 and ΩΛ = 1 − ΩM = 0.73. From the redshifts of the lens and the source galaxies in Sect. 3, one arcsecond at the lens (source) plane in RX J1131−1231 corresponds to 4.43 (7.04) kpc.

2. GLaD methodology

GLaD is a software developed to combine the stellar dynamical and strong lensing analyses with the aim of reconstructing the light and mass properties of the source galaxy. We use axisymmetric Jeans modelling to create a dynamical model of the source galaxy assuming a mass and a surface brightness distribution, and then predict how the source’s surface brightness and kinematics would look like when lensed into arcs. We then compare our predictions of the source’s distorted surface brightness and projected second order velocity moment v LOS 2 ¯ $ \overline{{v}_{\mathrm{LOS}}^{2}} $ with that of the observed arcs. In this method, we simultaneously fit the source and deflector properties, consisting of the source and deflector’s light and mass distributions, the source orbital anisotropy parameter β and inclination i. We use Bayesian analysis to infer the best-fit parameter values together with their uncertainties and degeneracies.

The GLaD method can be schematised as follows:

1. Assumption of a light profile, mass profile, anisotropy of stellar orbits and inclination of the mass profile for the source galaxy.

2. Prediction of the kinematic map and surface brightness map of the source on the source plane.

3. Assumption of a lens mass profile (and surface brightness, if desired).

4. Computation of the deflection angle for the particular the lens mass profile.

5. Mapping of the kinematic map and surface brightness map of the source onto the lens plane.

6. PSF convolution of the data on the image plane.

7. Comparison of predicted kinematic map and surface brightness map with observational data.

8. Iterative repetition of this process in case of optimisation and sampling of model parameters.

In this section we explore the different functionalities of GLaD: we introduce the lens profiles we employ in our analysis (Sect. 2.1), we describe the dynamics and lensing analysis (Sects. 2.2 and 2.3) and explain how these are combined to obtain a consistent joint analysis (Sect. 2.4).

2.1. Mass and light profiles

To describe the source and lens galaxy’s mass distribution, we use simply parametrised profiles. We describe the source’s total mass distribution with a pseudoisothermal elliptical mass distribution (PIEMD; Kassiola & Kovner 1993) with dimensionless surface mass density

κ piemd ( x , y ) | z s = ∞ = θ E , ∞ 2 R em 2 + r c 2 , $$ \begin{aligned} \kappa _{\rm piemd}(x,{ y})\Big |_{z_{\rm s}=\infty }= \frac{\theta _{\rm E,\infty }}{2 \sqrt{R_{\rm em}^2 +r_{\rm c}^2}}, \end{aligned} $$(1)

where (x, y) are the coordinates in the galaxy’s plane (along the semi-major and semi-minor axes), Rem is the elliptical mass radius,

R em = x 2 ( 1 + e ) 2 + y 2 ( 1 − e ) 2 , $$ \begin{aligned} R_{\rm em}= \root \of {\frac{x^2}{(1+e)^2} + \frac{{ y}^2}{(1-e)^2}}, \end{aligned} $$(2)

e is the ellipticity e = 1 − q 1 + q $ e= \frac{1-q}{1+q} $ with q the axis ratio, θE, ∞ is the lens Einstein radius for source at redshift infinity and rc is the core radius. The mass distribution is then suitably rotated by its position angle θ and shifted by the centroid position of the coordinate system used. The parameters that identify this profile are its centroid position (xc, yc), its axis ratio q, its position angle θ, its Einstein radius θE, ∞, and its core radius rc.

To represent the mass distribution of the lens we use a softened power-law elliptical mass distribution (SPEMD; Barkana 1998) which contains an additional parameter as compared to the PIEMD profile, namely the slope γ. Its convergence is

κ spemd ( x , y ) | z s = ∞ = θ E , ∞ ( x 2 + y 2 q 2 + r c 2 ) − γ , $$ \begin{aligned} \kappa _{\rm spemd}(x,{ y})\Big |_{z_{\rm s}=\infty }= \theta _{\rm E,\infty } \left(x^2+\frac{{ y}^2}{q^2}+r_{\rm c}^2 \right) ^{-\gamma }, \end{aligned} $$(3)

where q is the axis ratio, rc is the core radius, γ is the power law index, which is 0.5 for an isothermal profile (Barkana 1998). To the lens mass we also add a constant external shear, described by the lens potential parametrised by

ψ ext ( ϑ , φ ) = 1 2 γ ext , ∞ ϑ 2 cos ( 2 ( φ − ϕ ext ) ) , $$ \begin{aligned} \psi _{\rm ext}(\vartheta ,\varphi ) = \frac{1}{2} \gamma _{\rm ext, \infty } \vartheta ^2 \cos (2(\varphi -\phi _{\rm ext})), \end{aligned} $$(4)

where ϑ and φ are polar coordinates such that x = ϑcos(φ) and y = ϑsin(φ), γext, ∞ is the shear strength for source at redshift infinity and ϕext is the shear angle (ϕext = 0° corresponds to a shearing along the x-direction while ϕext = 90° corresponds to a shearing along the y-direction). The shear centre is arbitrary since it is not observable. Finally, to disentangle the baryonic mass from the dark component, we alternatively model the lens mass with a composite mass model. The composite model consists of multiple aforementioned PIEMD profiles (see “chameleon” profile below), to account for the contribution of the stellar density, and a Navarro et al. (1997, NFW) profile for the dark matter density distribution

ρ ( r ) = ρ s r r s ( 1 + r r s ) 2 , $$ \begin{aligned} \rho (r) = \frac{\rho _{\rm s}}{\frac{r}{r_{\rm s}} \left( 1 + \frac{r}{r_{\rm s}} \right)^{2}}, \end{aligned} $$(5)

where ρs is the characteristic overdensity and rs is the scale radius.

To represent both the lens and source light profiles, we use a Sérsic profile, whose intensity is given by

I ( R q ) = I e exp { − k [ ( R q R eff ) 1 n − 1 ] } , $$ \begin{aligned} I(R_{\rm q}) = I_{\rm e} \exp \Bigg \{ -k\ \bigg [ \bigg (\frac{R_{\rm q}}{R_{\rm eff}}\bigg )^{\frac{1}{n}}-1 \bigg ] \Bigg \}, \end{aligned} $$(6)

where k is approximately 2 n − 1 3 $ 2n-\frac{1}{3} $, R q = x 2 + y 2 / q 2 $ R_{\mathrm{q}}=\sqrt{x^2+\mathit{y}^2/q^2} $ and n is the Sérsic index which, for most galaxies, spans values between 1 2 < n < 10 $ \frac{1}{2} < n < 10 $ and whose value is correlated to the size and the magnitude of the galaxy. The Sérsic profile is defined by the centroid position xc, s and yc, s, the axis ratio q, the position angle θ, the Sérsic amplitude Ie, the effective radius Reff and the Sérsic index n. In addition to the Sérsic profile, we also employ the so-called “chameleon” profile, which is a profile that mimics the Sérsic and allows analytic computations of lensing quantities (e.g. Maller et al. 2000; Dutton et al. 2011; Suyu et al. 2014). This profile is composed of a difference of two isothermal profiles, namely

L ( x , y ) = L 0 1 + q l ( 1 x 2 + y 2 / q l 2 + 4 w c 2 / ( 1 + q l ) 2 − 1 x 2 + y 2 / q l 2 + 4 w t 2 / ( 1 + q l ) 2 ) , $$ \begin{aligned} L(x,{ y}) =&\frac{L_0}{1+q_{\rm l}} \left( \frac{1}{\sqrt{x^2 + { y}^2 / q_{\rm l}^2 + 4 { w}_{\rm c}^2 / (1+q_{\rm l})^2}} \right. \nonumber \\& - \left. \frac{1}{\sqrt{x^2 + { y}^2 / q_{\rm l}^2 + 4 { w}_{\rm t}^2 / (1+q_{\rm l})^2}} \right), \end{aligned} $$(7)

where ql is the axis ratio, wt and wc are parameters of the profile (with wt >  wc to keep L >  0). To represent the baryonic mass we scale the chameleon light profile by a constant mass-to-light ratio.

2.2. Dynamical modelling: axisymmetric Jeans modelling

The distribution function (DF) f(x, v) describing the positions x = (x, y, z) and velocities v = (vx, vy, vz) of a large system of stars must satisfy the fundamental equation of stellar dynamics, the steady-state collisionless Boltzmann equation (CBE; Binney & Tremaine 1988). With the axial symmetry assumption (in 3 dimensions), multiplication of the CBE with powers of the velocity moment, and subsequent integration over velocity space, we can reduce the Boltzmann equation into the Jeans equations (Jeans 1922), written in terms of the cylindrical coordinates (R, z, ϕ)

ν v R 2 ¯ − ν v ϕ 2 ¯ R + ∂ ( ν v R 2 ¯ ) ∂ R + ∂ ( ν v R ¯ v z ¯ ) ∂ z = − ν ∂ Φ ∂ R , $$ \begin{aligned}&\frac{\nu \overline{{v}^{2}_{R}}-\nu \overline{{v}^{2}_{\phi }}}{R} + \frac{\partial \left( \nu \overline{{v}^{2}_{R}} \right)}{\partial R} + \frac{ \partial \left( \nu \overline{{v}_{R}} \overline{{v}_{z}} \right)}{\partial z} = -\nu \frac{\partial \Phi }{\partial R}, \end{aligned} $$(8)

ν v R ¯ v z ¯ R + ∂ ( ν v z 2 ¯ ) ∂ z + ∂ ( ν v R ¯ v z ¯ ) ∂ R = − ν ∂ Φ ∂ z , $$ \begin{aligned}&\frac{\nu \overline{{v}_{R}} \overline{{v}_{z}}}{R} + \frac{\partial \left( \nu \overline{{v}^{2}_{z}} \right)}{\partial z} + \frac{ \partial \left(\nu \overline{{v}_{R}} \overline{{v}_{z}} \right) }{\partial R} = -\nu \frac{\partial \Phi }{\partial z}, \end{aligned} $$(9)

where Φ is the gravitational potential, ν the luminosity density and

ν v i ¯ v j ¯ ≡ ∫ v i v j f d 3 v . $$ \begin{aligned} \nu \overline{{v}_{i}} \overline{{v}_{j}} \equiv \int {{v}_{i} {v}_{j} f \mathrm{d}^{3} {\boldsymbol{v}}}. \end{aligned} $$(10)

The axisymmetry assumption seems to be valid, to first order, for most elliptical galaxies, unless photometric or kinematic evidence for bars or triaxiality is present. However, Eqs. (8) and (9) are still quite general and do not uniquely specify a solution (Cappellari 2008). By specifying the shape and orientation of the intersection of the velocity ellipsoid everywhere in the meridional plane, one can uniquely solve the Jeans equations numerically and recover the motion of stars in a gravitational potential Φ (Cappellari 2008). To solve the Jeans equations under the above assumptions, we use the Jeans Anisotropic MGE1 routine (Cappellari 2002, 2008; Cappellari & Copin 2003). This software makes the assumptions that: (1) the velocity ellipsoid is aligned with the cylindrical coordinate system (R, z, ϕ) and (2) the anisotropy is constant and quantified by v R 2 ¯ = β v z 2 ¯ $ \overline{{v}_{R}^{2}}= \beta\overline{{v}^{2}_{z}} $. Using these assumptions together with the boundary conditions that ν v ¯ z = 0 $ \nu \overline{{v}}_{z} =0 $ for z → ∞, one can finally solve the Jeans equations for v z 2 ¯ $ \overline{{v}_{z}^{2}} $ and v ϕ 2 ¯ $ \overline{{v}_{\phi}^{2}} $.

To derive solutions for the Jeans Eqs. (8) and (9), we first parametrise the source galaxy’s stellar surface brightness and surface mass density with a Multi-Gaussian Expansion (MGE). The MGE method was initially conceived by Bendinelli (1991) and then further developed in other works (Monnet et al. 1992; Emsellem et al. 1994, 1999; Cappellari 2002). It consists of a series of expansions of the galaxy’s image using 2D Gaussians functions

I ( R ∼ , θ ∼ ) = ∑ i = 1 N L i 2 π σ i 2 q i exp [ − 1 2 σ i 2 ( x 2 + y 2 q i 2 ) ] , $$ \begin{aligned} I(\tilde{R},\tilde{\theta }) = \displaystyle \sum _{i=1}^{N} \frac{L_{i}}{2\pi \sigma ^{2}_{i} q_{i}} \exp {\left[ - \frac{1}{2\sigma ^{2}_{i}} \left( x^{2} + \frac{{ y}^{2}}{q^{2}_{i}} \right) \right] }, \end{aligned} $$(11)

where N is the number of Gaussians with luminosity Li, axis ratio qi, and dispersion σi, (x, y) are a system of coordinates on the plane of the sky centred on the galaxy’s nucleus, and ( R ∼ , θ ∼ ) $ (\tilde{R},\tilde{\theta}) $ the relative polar coordinates ( x = R ∼ sin ( θ ∼ − ψ i ) , y = R ∼ cos ( θ ∼ − ψ i ) ) $ (x=\tilde{R}\sin({\tilde{\theta}-\psi_{i}}), \mathit{y}=\tilde{R}\cos({\tilde{\theta}-\psi_{i}})) $ with ψi the position angle measured counterclockwise from the y-axis to the major axis of the Gaussian. All Gaussians are assumed to have the same centre and position angle. This method allows for a straightforward and analytically convenient expression of an arbitrary surface brightness and surface mass density distribution (Cappellari 2002). Since galaxies have an unknown inclination i, one needs to deproject the MGE surface brightness and surface mass density profile to get the intrinsic tracer and mass density ν and Φ in Eqs. (8) and (9). Despite not being able to eliminate the intrinsic degeneracy of the deprojection, the MGE method provides realistic densities, resembling real galaxies, when projected at any angle (Cappellari 2008). With the deprojected MGE of the surface brightness and surface mass density profile, one can readily solve the Jeans equations and perform the integral along the line-of-sight to obtain the total observed first and second order velocity moment

v rms = v 2 + σ v 2 , $$ \begin{aligned} {v}_{\rm rms}=\sqrt{{v}^{2}+\sigma _{\rm v}^2}, \end{aligned} $$(12)

where v and σv are the source galaxy’s projected velocity and velocity dispersion respectively. As it is not clear how the ordered and random motions contribute to a particular vrms profile a priori, we stick to the prediction of projected second order moment in Eq. (12) for our modelling purposes and show an example of the predicted velocity moments (vrms) simulated with GLaD in Fig. 1.

thumbnail Fig. 1.

Simulated source stellar kinematics with GLaD. We produce an axisymmetric Jeans model of a mock source galaxy on the source plane (left panel) and we predict the lensed vrms profile in the image plane (right panel). The masked black region corresponds to regions in the image plane which have no correspondence on the source plane. The yellow cross marks the lens centroid position.

2.3. Gravitational lensing: GLEE

Our lensing analysis is carried out using GLEE, a software developed by A. Halkola and S. H. Suyu (Suyu & Halkola 2010; Suyu et al. 2012). This software uses parametrised mass profiles (discussed in Sect. 2.1) to describe the different lensing components, such as dark matter halos and galaxies, and allows us to compute the deflection angle for these profiles to map between the source and the image planes. The deflection angle is calculated as

α ̂ ( ξ ) = 4 G c 2 ∫ d 2 ξ ′ Σ ( ξ ′ ) ξ − ξ ′ | ξ − ξ ′ | 2 , $$ \begin{aligned} \hat{\alpha }(\xi ) = \frac{4G}{c^{2}} \int \mathrm{d}^{2} {\xi ^{\prime }} \Sigma (\boldsymbol{\xi }^{\prime }) \frac{\boldsymbol{\xi }-\boldsymbol{\xi }^{\prime }}{{|\boldsymbol{\xi }-\boldsymbol{\xi }^{\prime }|}^2}, \end{aligned} $$(13)

where G is the gravitational constant, ξ is the 2-dimensional impact vector, and Σ is the surface mass density of the galaxy, namely the mass density projected onto a plane perpendicular to the incoming light ray, defined as Σ = κΣcrit, where κ is the dimensionless surface mass density (convergence) and the critical surface mass density is Σ crit = c 2 4 π G D s D d D ds $ \Sigma_{\mathrm{crit}}= \frac{c^{2}}{4\pi G} \frac{D_{\mathrm{s}}}{D_{\mathrm{d}}D_{\mathrm{ds}}} $, which is a function of the angular diameter distances of lens and source.

For the development of GLaD we use only these features of GLEE, namely we use GLEE to (1) compute the deflection angles of analytic mass distributions, (2) model surface brightness profiles of galaxies, and (3) set up multi-lens-plane modelling. GLaD further builds upon this by incorporating dynamical modelling of the sources. Therefore, a key difference between GLEE and GLaD is on the model of the source: GLEE reconstructs the source surface brightness on a grid of pixels (Suyu et al. 2006), whereas GLaD uses analytic profiles for both the source surface brightness and mass distributions. An example of the use of GLEE to construct the source and lens surface brightnesses is shown in Fig. 2.

thumbnail Fig. 2.

Simulated strong lensing data with GLaD. We model the light of the source (left panel) and the mass and light of the lens and we predict how the surface brightness of the lensed source (central panel) and lensed source with lens (right panel) will look like. The units of the surface brightness are in counts.

2.4. Joint analysis

To sample the parameter space efficiently and obtain the posterior probability distributions we use EMCEE (Foreman-Mackey et al. 2013), a stable, well tested Python implementation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) proposed by Goodman & Weare (2010). The posterior probability of the model parameters, collectively denoted by η, is obtained using Bayes’ Theorem

P ( η | X obs ) ∝ L ( X obs | η ) ⏞ likelihood P ( η ) ⏞ prior , $$ \begin{aligned} P(\boldsymbol{\eta }|\boldsymbol{X}_{\rm obs}) \propto \overbrace{\mathcal{L} (\boldsymbol{X}_{\rm obs}|\boldsymbol{\eta })}^\mathrm{likelihood} \overbrace{P(\boldsymbol{\eta })}^\mathrm{prior}, \end{aligned} $$(14)

where P(η) is the prior probability on the model parameters, that we always assume to be uniform, and Xobs collectively denotes our observables (surface brightness Iobs and velocity map v rms obs $ \boldsymbol{v}^{\mathrm{obs}}_{\mathrm{rms}} $). To perform the joint analysis, we combine the lensing likelihood ℒlens with the dynamics likelihood ℒdyn such that the total likelihood function would be

L tot = L lens × L dyn . $$ \begin{aligned} \mathcal{L} _{\rm tot}= \mathcal{L} _{\rm lens} \times \mathcal{L} _{\rm dyn}. \end{aligned} $$(15)

The lensing likelihood is

L lens ( I obs | η ) ∝ exp [ − 1 2 ∑ i = 1 N pixels | I i obs − I i pred ( η ) | 2 σ i 2 ] , $$ \begin{aligned} \mathcal{L} _{\rm lens}(\boldsymbol{I}_{\rm obs}|\boldsymbol{\eta }) \propto \exp {\left[-\frac{1}{2} \displaystyle \sum _{i=1}^{N_{\rm pixels}} \frac{\vert {I}_{i}^\mathrm{obs}-{I}_{i}^\mathrm{pred}(\boldsymbol{\eta })\vert ^2}{\sigma _{i}^2} \right]}, \end{aligned} $$(16)

where Npixels is the number of pixels in the image, I i obs $ {I}_{i}^{\mathrm{obs}} $ is the observed surface brightness in a certain pixel i, I i pred ( η ) $ I_{i}^{\mathrm{pred}}(\boldsymbol{\eta}) $ is the modelled surface brightness of that same pixel and σi is the uncertainty on that pixel. Here, I and σ have the units of counts. We assume our uncertainty in every pixel σi to be composed of the background noise (with σback) and by the Poisson noise, namely

σ i = σ back 2 + I i obs . $$ \begin{aligned} \sigma _{i}= \sqrt{\sigma _{\rm back}^{2} + I^\mathrm{obs}_{i}}. \end{aligned} $$(17)

The lensing likelihood can be easily generalised to account for correlated noise among the pixels, by substituting a covariance matrix to σ in Eq. (16). The dynamics likelihood is

L dyn ( v rms obs | η ) ∝ exp [ − 1 2 ∑ i = 1 N bin | v rms , i obs − v rms , i pred ( η ) | 2 σ rms , i 2 ] , $$ \begin{aligned} \mathcal{L} _{\rm dyn}\left({\boldsymbol{v}}^\mathrm{obs}_{\rm rms} |{\boldsymbol{\eta }}\right) \propto \exp {\left[-\frac{1}{2}\displaystyle \sum _{i=1}^{N_{\rm bin}} \frac{\vert {{v}}_{\mathrm{rms},i}^\mathrm{obs}-{{v}}_{\mathrm{rms},i}^\mathrm{pred}(\boldsymbol{\eta })\vert ^2}{\sigma _{\mathrm{rms},i}^2} \right]}, \end{aligned} $$(18)

where v rms , i obs $ {{v}}_{\mathrm{rms},i}^{\mathrm{obs}} $ is the observed value of Eq. (12) in each bin, v rms , i pred $ {{v}}_{\mathrm{rms},i}^{\mathrm{pred}} $ our model prediction, and σrms, i the error on each bin. To obtain the value of the latter, we assume that the error in each IFU pixel scales inversely proportional to its S/N (see Emsellem et al. 2004). Since we do not mock up the spectroscopic data, we rely on the surface brightness information of the imaging data for assessing the S/N in each IFU spaxel. When binning the vrms, the signal to noise in each bin is

S N = ∑ j = 1 N pixels I j ∑ j = 1 N pixels ( σ back , j 2 + I j ) , $$ \begin{aligned} \frac{S}{N}= \frac{\sum _{j=1}^\mathrm{N_{\rm pixels}} I_{j}}{\sqrt{\sum _{j=1}^{N_{\rm pixels}} \left( \sigma _{\mathrm{back},j}^{2} + I_{j} \right)}}, \end{aligned} $$(19)

where Npixels is the number of pixels in each bin and I and σ here have the units of counts. The error on the vrms value in each bin is

σ rms = v rms S N · $$ \begin{aligned} \sigma _{\mathrm{rms}}= \frac{{{v}}_{\rm rms}}{\frac{S}{N}}\cdot \end{aligned} $$(20)

3. Demonstration: RX J1131−1231

RX J1131−1231 is a gravitational lensing system discovered by Sluse et al. (2003). This system is composed of four multiple images of a distant quasar and its host galaxy (zs = 0.654, Sluse et al. 2007), lensed by an intervening giant elliptical galaxy at redshift zd = 0.295 (Sluse et al. 2003, 2007), as shown in Fig. 3. The quadruply imaged quasar is surrounded by an Einstein ring of ∼3″ diameter. This system was used for studies on quasars and the region around black holes (e.g. Dai et al. 2010), on time-delay cosmography (e.g. Morgan et al. 2006; Suyu et al. 2013, 2014; Birrer et al. 2016; Chen et al. 2019), on dark matter substructures (e.g. Birrer et al. 2017, among others), on planet searches with microlensing (e.g. Dai & Guerras 2018), and on black hole and galaxy co-evolution (e.g. Ding et al. 2017). Thanks to the large Einstein ring and the high magnification, the lensed quasar host promises to be an interesting test bed for obtaining IFU stellar kinematic data with future IFU instruments such as the Near Infrared Spectrograph (NIRSpec; Bagnasco et al. 2007; Birkmann et al. 2016) on the James Webb Space Telescope (JWST). We estimate, through the Exposure Time Calculator2 (ETC), that we would need an on-source integration time of ∼6.5 h to obtain a signal-to-noise of ∼11 in the brightest pixel of the lensed arcs with NIRSpec. Through binning, we can easily increase the signal-to-noise ratio to 20, which is deemed necessary to properly extract the kinematic information from the data (Falcón-Barroso et al. 2017). This shows how measuring spatially resolved stellar kinematics of galaxies in 2D at z ≳ 0.6 would soon become feasible.

thumbnail Fig. 3.

HST ACS image of RX J1131−1231 in the F814W filter. The lens galaxy G is surrounded by an Einstein ring and the four images of the quasar labelled with letters (A, B, C, D). The satellite in the northern part of the lens galaxy is not included in our mock data. The image resolution is 0.05″ pixel−1.

3.1. RX J1131−1231 simulated data

We construct mock lensing and kinematic data for RX J1131−1231 starting from the best-fit lens mass model obtained by fitting both the quasar image positions as point sources as well as the full Einstein ring as an extended source (Yıldırım et al. 2020) to the Hubble Space Telescope (HST) Advanced Camera for Surveys (ACS) image of RX J1131−1231 in the F814W filter. The model of the lens mass distribution is a power law with an external shear (parametrised as described in Sect. 2.1). The power law profile was already proven to be a good fit for this lensing system in previous studies (among others Suyu et al. 2013, 2014) and, in general, power law profiles are a good representation of galaxies’ mass distributions according to studies based on X-ray observations (Humphrey & Buote 2010) and the Sloan Lens ACS survey (e.g. Koopmans et al. 2006, 2009; Gavazzi et al. 2007; Auger et al. 2010; Barnabè et al. 2011) and from studies on lens potential corrections (Suyu et al. 2009). The best-fit parameters of the lens mass model for simulating our data are presented in Table 1. To probe dependencies on the parametrisation of the lens mass model we also explore a composite mass distribution of baryons and dark matter in Sect. 3.2.4.

Table 1.

Input lens mass parameters for RX J1131−1231 obtained by fitting both the quasar image positions as point sources as well as the full Einstein ring as an extended source (Yıldırım et al. 2020).

3.1.1. Mock lensing data

Given the lens mass model described in Sect. 3.1, we predict the centroid position of the quasar host galaxy on the source plane using GLEE. We parametrise the source by a Sérsic in order to obtain a smooth surface brightness distribution, that we use for the dynamical analysis as well. We centre the source profile on the source centroid obtained with GLEE, and assume an axis ratio, orientation, effective radius and a Sérsic index. We lens the so obtained source surface brightness map to the lens plane. To simulate the lens light, we centre a Sérsic profile on the centroid coordinates of the lens mass profile described in Sect. 3.1, and we also assume the lens light orientation to be the same as that of the lens mass. We convolve the model with a Gaussian PSF with standard deviation of 4 pixels, which we estimated from measuring the FWHM of different stars in the field. Finally, we adjust the amplitudes of both light profiles, to obtain a surface brightness image that mimics the HST observations both in terms of the surface brightness distribution and the S/N. Since we mainly focus on the host galaxy’s arcs, we have excluded the quasar’s light from our model. We add simulated noise, that we obtain using Eq. (17). The resulting surface brightness model is shown in Fig. 4 and the input parameters are shown in Table 2. The data image has a total of 14641 pixels.

thumbnail Fig. 4.

Mock surface brightness image of RX J1131−1231 obtained, as explained in Sect. 3.1.1, using the parameters in Tables 1 and 2. The lens galaxy at the centre is surrounded by an Einstein ring and the four images of the host galaxy (quasar excluded). The image has a resolution of 0.05″ pixel−1. In the mock image we added simulated noise, as explained in Sect. 3.1.1.

Table 2.

Input Sérsic light parameters of the mock lens and source galaxies of RX J1131−1231.

3.1.2. Mock kinematic data

Once we have produced a mock lensing image to use as constraints, we produce mock kinematic data for the lensed arcs. To this end, we have to assume a mass profile for the source, in addition to the source (Sérsic) light profile described in Sect. 3.1.1. We adopt a pseudoisothermal elliptic mass distribution (PIEMD). We use a single profile for the combined dark and baryonic matter components given that the limited kinematic data in the arcs is most likely insufficient to effectively break the degeneracies between the two components. The parameters of this PIEMD profile are given in Table 3. We adopt the same centroid and same inclination for both the light and mass profiles of the source. We obtain the Einstein radius of the source mass distribution using the equation

θ E , ∞ = 4 π σ v 2 c 2 , $$ \begin{aligned} \theta _{\rm E,\infty } = \frac{4 \pi \sigma _{\rm v}^{2}}{c^2}, \end{aligned} $$(21)

Table 3.

Input source mass parameters of the mock RX J1131−1231.

where σv is the mean velocity dispersion that we assume to be 200 km s−1 for the source galaxy, a sensible estimate for a galaxy at redshift zs = 0.654 hosting a black hole, whose mass we fixed to around 108 M⊙ (Dai et al. 2010). Since we are studying the properties of the source, which does not act itself as a lens to any background galaxy, we compute all the quantities for an “artificial” source at redshift infinity, namely we fix the D ds D s = 1 $ \frac{D_{\mathrm{ds}}}{D_{\mathrm{s}}}=1 $, where Dds is the distance from our source at redshift zs = 0.654 to the artificial source at redshift infinity, and Ds is the distance to the artificial source. We show the source’s mass parameters in Table 3.

With the aforementioned parametrised profiles for the source mass and light, we produce an MGE model of the source light and mass on the source plane. We use the JAM routine to produce a model of the vrms of the source galaxy on the source plane, and lens the kinematic map to the image plane. We then convolve with a gaussian PSF with zero mean and standard deviation of 0.1″. Our final mock 2D kinematic map has a resolution of 0.1″, to resemble the JWST NIRSpec IFU instrument. As already mentioned in the introduction to this section, we estimated the signal-to-noise ratio in the arcs to be ∼11 from 6.5 h of observations with JWST NIRSpec according to the ETC. In order to obtain a more favourable balance between S/N and spatial resolution of our source kinematic data, while keeping the total integration time still reasonable, we scale the surface brightness of the source light profile such that the signal-to-noise ratio of the brightest pixel on the image plane becomes 15. To do that, we multiply the intensity by a factor Q such that

( S N ) bp = I bp Q σ back 2 + I bp Q = 15 , $$ \begin{aligned} \left( \frac{S}{N} \right)_{\rm bp}= \frac{I_{\rm bp} Q}{\sqrt{ \sigma _{\rm back}^{2} + I_{\rm bp} Q}}=15, \end{aligned} $$(22)

where “bp” denotes brightest pixel. We solve Eq. (20) to obtain the Q factor and then scale the intensity of each pixel Ii and the Poisson noise component of each pixel by multiplying them with the Q factor. We therefore obtain a new rescaled intensity map with a signal-to-noise of 15 in the brightest pixel.

We use Voronoi binning as described in Cappellari & Copin (2003), an adaptive spatial binning technique, to obtain the signal-to-noise ratio per bin of ∼20 (labelled SN20) and ∼30 (labelled SN30), necessary to obtain reliable measurements of the stellar kinematics from the spectroscopic data. The binned kinematic map is illustrated in Figs. 5 and 6 for the cases of SN20 and SN30, respectively. For simplicity, we assume that all bins contribute equally as constraints throughout the analysis. Realistically, however, the bins further away from the arc might not only have a lower constraining capacity than those closer to the arc, but it will be highly difficult to extract valuable kinematic information from the spectra, given the low S/N spaxels in these regions and the massive binning that would be needed to push the S/N above the threshold of 20 and 30, respectively. Therefore, before proceeding further with this assumption, we tested the case in which those bins have an artificially higher error on their kinematic moments, which effectively discards them from the fit. Since we noticed no significant discrepancy in the final parameter constraints, we assume all the bins to contribute equally throughout the analysis.

thumbnail Fig. 5.

Mock binned lensed kinematic data of the source galaxy of RX J1131−1231. We scale the signal of the source galaxy to obtain a signal-to-noise ratio in the brightest pixel of 15, to mimic a realistic scenario (see Sect. 3.1.2 for further details). The data are binned to obtain a signal-to-noise in each bin of 20. We obtain a total of 45 bins. The pixel resolution of the kinematic map is 0.1″. The masked black region corresponds to regions in the image plane which have no correspondence on the source plane. The yellow cross marks the lens centroid position.

thumbnail Fig. 6.

Mock lensed kinematic data of the source galaxy of RX J1131−1231. We scale the signal of the source galaxy to obtain a signal-to-noise ratio in the brightest pixel of 15, to mimic a realistic scenario (see Sect. 3.1.2 for further details). The data are binned to obtain a signal-to-noise in each bin of 30. We obtain a total of 19 bins. The pixel resolution of the kinematic map is 0.1″. The masked black region corresponds to regions in the image plane which have no correspondence on the source plane. The yellow cross marks the lens centroid position.

3.2. Mass models of mock data

Once we have simulated a set of mock lensing and kinematic data, we perform different tests to assess how well we are able to recover the lens and source parameters. We sample the parameter space by varying the lens and source mass and light distributions. To assess the constraining power of the combination of the lensing and kinematic data, we model all the parameters using both the lensing and dynamics data (labelled LD) and using the lensing data only (labelled L). We also model only the lensing parameters (i.e. fixing the source mass) by using only the lensing data as constraints. We test this for both the SN20 and SN30 cases, to assess the improvement we obtain with a higher amount of kinematic constraints (model SN20 has more than double the number of bins than model SN30). In addition, we explore a model that has the source’s mass distribution following its light (scaled by a M/L), which we denote as a mass-follows-light model (MFL), to assess systematic errors associated with the source mass parameterisation. Finally, we consider a different parametrisation of the lens mass distribution to determine the impact of imperfect lens mass model on the inference of the source properties.

3.2.1. LD: Lensing and dynamical models

Our first test consists of remodelling all the source’s and lens’s mass and light parameters using as constraints both the lensing and kinematic data, to check if GLaD is able to consistently recover the input parameters. We use both the kinematic data with SN20 and with SN30 (see Sect. 3.1.2). As shown in Table 4, for both SN20 and SN30 most of the parameters are recovered within the 1σ uncertainties, and all are recovered within the 2σ, as shown in Fig. 7; both of the models have χ red 2 ~1 $ \chi^2_{\rm red}\sim 1 $. We note a tight anticorrelation between the Sérsic amplitude Ie, s and the Sérsic index n of the source light, and between the source mass inclination and position angle. We also find a tight anticorrelation between the Sérsic amplitude Ie, s and the effective radius of the lens light profile, and between the lens Einstein radius and the slope of the power-law profile. The source mass parameters do not show correlation with other parameters, apart from a mild anticorrelation between the source axis ratio and Einstein radius. Moreover, the inclination is not well constrained. We do not find a significant improvement on the lens galaxy constraints when using different kinematic data quality (SN20 and SN30), probably because the amount of kinematic data points (19 in the case of SN30 and 45 in the case of SN20) is still subdominant with respect to the lensing constraints (14641 pixels). However, we note a factor ∼1.1−1.2 improvement on the source parameters such as Sérsic amplitude Ie, s, effective radius Reff, s, Sérsic index ns, source’s mass axis ratio q s piemd $ q^{\mathrm{piemd}}_{\mathrm{s}} $ and Einstein radius θ E , s piemd $ \theta^{\mathrm{piemd}}_{\mathrm{E,s}} $ in the SN20 case. We expect this improvement to increase with higher signal-to-noise per pixel on the kinematic data, and a consequently higher amount of bins. Therefore, for a given S/N in the kinematic data, it is preferable to spatially resolve as much as possible at the expense of a lower S/N per bin, provided that each bin has sufficient S/N to yield accurate vrms measurements.

thumbnail Fig. 7.

Joint 2D posterior probability distribution for the parameters of the power-law mock (presented in Sect. 3.1), with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin. The different contours in the 2D plots indicate, respectively, the 1σ, 2σ and 3σ CIs. Parameters shown are those where the improvement on the constraints coming from the combination of lensing and dynamics (LD in red contours) is more prominent as compared to lensing only (L in grey contours). These parameters are the source centroid xc, s, yc, s, the anisotropy β, the source and lens light axis ratio qs, ql, the lens Einstein radius θ E , ∞ , l spemd $ \theta^{\mathrm{spemd}}_{\mathrm{E,\infty,l}} $ and slope γ l spemd $ \gamma^{\mathrm{spemd}}_{\mathrm{l}} $, and the shear parameters γext, ∞ and ϕext. In the diagonal are shown the 1D histograms of the corresponding parameter on the x-axis. The input values are marked with a yellow cross.

Table 4.

Lens and source parameter constraints for the different configurations for our first set of mock data (presented in Sect. 3.1).

3.2.2. L: Lensing only models

In the second test we perform, we include only the lensing surface brightness data as constraints, in order to assess the changes when compared to the joint lensing and dynamical modelling in Sect. 3.2.1. As shown in Table 4 and in Fig. 7, parameters like the source mass and light centroid, the source Sérsic amplitude Ie, s and effective radius are better constrained (up to an order of magnitude tighter) in the combined lensing and dynamics case. Indeed, the dynamical analysis makes direct use of these quantities to predict the model of the lensed source kinematic map. Moreover, we note that also the lens amplitude, lens mass axis ratio and Einstein radius are constrained better up to a factor of 3 in the combined lensing and dynamics case. This shows that actually the addition of the kinematic data will put tighter constraints on the lens parameters as well. We also find that the addition of the kinematic data puts tighter constraints on the shear parameters (amplitude and orientation), with respect to lensing only case. As expected, when trying to vary the source mass parameters not including the kinematics analysis, we find them to be virtually unconstrained, as shown in Table 4 and in Fig. 8. Finally, we note a tight anticorrelation between the Sérsic amplitude Ie, s and the Sérsic index n of the source light, and between the source mass inclination and position angle, as already seen in the LD case. For this model, the input parameters are mostly recovered within the 1σ uncertainties, and are all recovered within the 2σ, with a χ red 2 ~1 $ \chi^2_{\rm red}\sim 1 $.

thumbnail Fig. 8.

Joint 2D posterior probability distribution for the source mass parameters of the power-law mock data presented in Sect. 3.1, with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin. The different contours in the 2D plots indicate the different σ CIs. Parameters shown are those where the improvement on the constraints coming from the addition of lensing and dynamics (LD in red contours) is more prominent as compared to lensing only (L in grey contours). These parameters are the anisotropy β, the inclination i, the source mass axis ratio q s piemd $ q^{\mathrm{piemd}}_{\mathrm{s}} $, the source Einstein radius θ E , ∞ , s piemd $ \theta^{\mathrm{piemd}}_{\mathrm{E,\infty, s}} $. In the diagonal are shown the 1D histograms of the corresponding parameter on the x-axis. The input values are marked with a yellow cross.

3.2.3. MFL: Mass-follows-light models

To test systematic modelling uncertainties, we model the source mass (originally simulated as an isothermal profile) using a different model. We use a mass-follows-light model, that is to say we scale the source light profile (Sérsic) using a mass-to-light ratio, and we allow the mass-to-light ratio to vary together with the other parameters. We fit to the lensing and kinematic data (those binned to obtain a signal-to-noise of 30 in each bin, SN30), and we test how well this model is able to reproduce the original data.

We find that assuming a different source mass profile biases some of the model parameters, which are not recovered within the uncertainties, as shown in Table 4. In this particular case, the MFL model cannot recover most of the source light, mass and dynamical parameters within the 1σ uncertainties, but most of them are recovered within the 2σ. The only parameter which is overestimated by ∼4σ is the anisotropy. The source inclination, Sérsic index, and light intensity seem to be more robust parameters and are less affected by systematics in this particular case, as they are mainly anchored by the lensing data. For the lens, we find that the parameters affected by systematics are the Sérsic intensity and the mass Einstein radius, ellipticity and slope, which are all tightly correlated to each other, as already mentioned for the LD case in Sect. 3.2.1. Finally, the external shear parameters are not recovered within the 3σ uncertainties. This might be the cause of the shift of the source centroid prediction, which is sensitive to the shear strength value, and explain the misfit of the kinematic data.

This model has a slightly higher χ red 2 $ \chi^2_{\rm red} $ than the other models ( χ red 2 ~1.07 $ \chi^2_{\rm red}\sim1.07 $), showing that it does not fit as well to the data. In particular, the dynamics χ dyn,red 2 $ \chi^2_{\rm dyn,red} $ is ∼1.63, while the lensing χ len,red 2 $ \chi^2_{\rm len,red} $ is ∼1.07, which tells us that the kinematic data are poorly fitted as compared to the lensing data. Indeed, as shown in Fig. 9, the vrms is underestimated, especially in the bins around the arc region. The surface mass density Σ(< R), shown in Fig. 10, is in good agreement with the data for radii smaller than 0.2 Reff, and underestimated otherwise, by up to a factor 2. We find that the size of the region within 0.2 Reff from the centre of the source on the source plane maps exactly to the location of the brighest bins (those with high signal-to-noise) on the image plane. This indicates that our constraints are robust in the region of high signal-to-noise, and instead model dependent on the regions with low signal-to noise. Therefore, we conclude that the addition of dynamics can help better distinguish between models, whereas the lensing analysis alone provides no constraints on the source mass at all.

thumbnail Fig. 9.

Mock lensed kinematic data of the source galaxy of RX J1131−1231 obtained using an isothermal mass profile (left panel), reconstructed kinematic map from our best-fit MFL model (central panel), and normalised residuals (right panel). The pixel resolution of the kinematic map is 0.1″. The masked black region is the region in the image plane which has no correspondence on the source plane. The yellow cross marks the lens centroid position. If we compare to the data, we see that our model predicted source’s vrms is slightly underestimated.

thumbnail Fig. 10.

Predicted averaged circularised surface mass density Σ(< R) (left) and mass enclosed M(< R) (right) within one effective radius Reff = 0.5″ (3.5 kpc) of the source galaxy, with the corresponding CI (shaded regions). Both the convergence and the enclosed mass of the source of the MFL model (see Sect. 3.2.3) are underestimated in the outskirts, as compared to the input and the LD model. This is also shown in the vrms value, which is slightly underestimated, as shown in Fig. 9.

3.2.4. Composite models

As a final test, we assess the improvement of the constraints on the lens mass model, particularly focussing on the dark matter component. To do so, we re-simulate both the lensing and kinematics data, as done in Sect. 3.1, assuming a composite mass model for the lens. In particular, to represent the baryonic mass we use a chameleon profile (discussed in Sect. 2.1) that we scale with a mass-to-light ratio, and for the dark matter we assume a NFW profile. We still assume an external shear component. For this model, the lens Sérsic light and the lens mass are decoupled. We show the parameters of this model in Table 5. For the kinematic data, we imposed a signal-to-noise of 30 in each bin (SN30), which allows us to obtain 22 kinematic data points (shown in Fig. 11). We remodel these new simulated data with both the full lensing and dynamics configuration and the lensing only configuration, to assess the improvement due to the addition of the kinematic constraints. Then we remodel the lens mass using a power-law model, to assess systematic uncertainties associated with the lens mass parameterisation. As shown from Table 5 and Fig. 12, we find improved constraints on the source centroid when we add the kinematic constraints, consistent with our findings for the power-law model. Moreover, the dark matter (NFW profile) orientation and Einstein radius are better constrained when including the kinematic data, compared to the lensing-only analysis. Both the LD composite and the L composite model have a χ red 2 ~1 $ \chi^2_{\rm red}\sim 1 $. When testing systematics errors for this set of mock data with the composite model, namely when modelling the lens composite mass with a single power-law model, we find that the addition of the kinematic data allows to better discern between models. Indeed, as shown in Table 5, the LD power-law model has a comparable total χ red 2 $ \chi_{\rm red}^2 $( χ red 2 =1.03 $ \chi_{\rm red}^2=1.03 $) than the lensing only L power-law ( χ red 2 ~1.02 $ \chi_{\rm red}^2\sim 1.02 $). However, the LD power-law has a higher misfit of the kinematic data χ dyn,red 2 $ \chi^2_{\rm dyn,red} $ ∼ 1.25. Moreover, in the LD power-law we find a better constrained lens Einstein radius by a factor 2 as compared to the L power-law. This is consistent with the previous set of mock data. Finally, as already noted previously, the source parameters are constrained when including kinematic data, in both the input set-up and the systematic test. In general, systematic uncertainties, from profile mismatch, dominate statistical uncertainties. Encouragingly, most of the mass parameters of the sources are recovered within 1σ uncertainty in the LD power-law model, showing that the source mass distribution could still be inferred from the kinematics despite differences in the lens mass profile between power-law and composite profiles.

thumbnail Fig. 11.

Mock lensed kinematic data of the source galaxy of RX J1131−1231 for the composite mass model (left panel), our prediction obtained by modelling the data with a power law profile (central panel) and the normalised residuals (right panel). The data are binned to obtain a signal-to-noise in each bin of 30, for a total of 22 bins. The pixel resolution of the kinematic maps is 0.1″. The masked black region is the region in the image plane which has no correspondence on the source plane. The yellow cross marks the lens centroid position.

thumbnail Fig. 12.

Joint 2D posterior probability distribution for the parameters of the composite mock (presented in Sect. 3.2.4), with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin (shown in Fig. 11). The different contours in the 2D plots indicate, respectively, the 1σ, 2σ and 3σ ranges. Parameters shown are those where the improvement on the constraints coming from the combination of lensing and dynamics (LD in red contours) is more prominent as compared to lensing only (L in grey contours). These parameters are the source centroid xc, s, yc, s, the anisotropy β, the position angle of the lens dark matter profile θ l nfw $ \theta^{\mathrm{nfw}}_{\mathrm{l}} $, its Einstein radius θ E , ∞ , l nfw $ \theta^{\mathrm{nfw}}_{\mathrm{E,\infty,l}} $, and the shear parameters γext, ∞ and ϕext. In the diagonal are shown the 1D histograms of the corresponding parameter on the x-axis. The input values are marked with a yellow cross.

Table 5.

Modelled lens and source parameter values for the different configurations of the composite mock data (presented in Sect. 3.2.4), with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin (shown in Fig. 11).

To show more in detail the differences among these models, we compare their average surface mass density3, namely

Σ ¯ ( < R ) = ∫ 0 R Σ ( R ′ ) 2 π R ′ d R ′ π R 2 · $$ \begin{aligned} \bar{\Sigma }({ < }R) = \frac{\int _{0}^{R} \! \Sigma (R^\prime ) 2 \pi R^\prime \, \mathrm{d}R^\prime }{\pi R^2}\cdot \end{aligned} $$(23)

where Σ(R) is the surface mass density. The Σ ¯ ( < R ) $ \bar{\Sigma}({ < }R) $ for the input, LD composite and L composite, is shown in the left panel of Fig. 13, while the Σ ¯ ( < R ) $ \bar{\Sigma}({ < }R) $ for the LD power-law and L power-law models is shown in the right panel. From Fig. 13 we see that the LD composite and L composite models are very consistent with the input in terms of average surface mass density. When we model with the power-law lens mass model, we still see a strong consistency with the input for the Σ ¯ ( < R ) $ \bar{\Sigma} ({ < }R) $. We note mild inconsistencies at R <  0.3″, where the LD power-law and L power-law models get peakier than the input.

thumbnail Fig. 13.

Average surface mass density Σ ¯ ( R ) $ \bar{\Sigma}(R) $ for the LD composite and the L composite models (left panel), and for the power-law models LD power-law and L power-law (right panel) as compared to the input composite mock, plotted with the relative error bars. All models recover the input closely at the value of the Einstein radius (∼1.6″ for source at reshift 0.654). On the left panel, the L composite model gets peakier towards the centre, but is still able to recover the input within the error bars (light blue shaded region). This does not appear to be the case for the LD power-law and L power-law (right panel), which differ from the input for R <  0.3″. However, at larger radii, we cannot easily distinguish between power-law and composite models from the average surface mass density.

Therefore we conclude that, for this system, the composite and power law models are not distinguishable in terms of Σ ¯ ( < R ) $ \bar{\Sigma} ({ < }R) $ and of the lensing analysis, but are thanks to the use of the dynamics analysis, which shows the higher misfit when a different lensing model than the input is used.

4. Summary

GLaD is a software that performs a joint gravitational lensing and dynamical modelling analysis, with the goal of studying properties of source galaxies and to assess the improvement on the constraints on the lens mass model when high-spatially resolved kinematic observations of the lensed source are avail- able. GLaD assumes parametrised mass and light profiles for both the source and the lens and produces a kinematic map (employing Jeans Anisotropic MGE) and a surface brightness map on the source plane. Then, by computing the deflection produced by the lens mass profile, these quantities are mapped to the image plane. We therefore obtain a lensed kinematic map (binned, if binning is performed on the input data) and a surface brightness map that we then compare to the input kinematic and surface brightness data. We show that this method allows one to infer the mass properties of the source and also to improve constraints on the lens mass parameters.

We test this method by simulating two different sets of mock data of the RX J1131−1231 system, and we remodel all the parameters with different assumptions. The first set of mock data assumes a power-law profile for the lens mass plus external shear, including both the baryonic and dark matter. The second set assumes a composite mass model for the lens mass, namely uses a chameleon light profile scaled by a mass-to-light ratio to describe the baryonic mass and a NFW profile for the dark matter component. With the first mock data set we test the improvement on the constraints on both lens and source due to the addition of the kinematic data as compared to a lensing-only analysis. We compare different quality of kinematic data and we test systematics on the source mass model by modelling it with different profiles. Finally we perform the previous analysis using the second mock data set. In this case, we also test systematic errors on the lens mass model, by remodelling the mass with a single power-law profile. Our findings are summarised below:

  • Our software is consistently able to recover the input parameter of the model within the 1−2σ uncertainties when modelling with the input configuration set-up, both in the case of the combined lensing and dynamics analysis and in the lensing-only analysis. This is true for both sets of mock data presented.

  • The addition of the source kinematics data allows us to put constraints on the source mass, that is otherwise unconstrained. Moreover, the joint lensing and dynamics analysis tightens constraints on the source centroid and the source light axis ratio (up to a factor ∼20). Furthermore, the statistical uncertainties on the lens mass axis ratio, Einstein radius and slope are improved by up to a factor ∼3. Tighter constraints are also observed on the shear parameters (up to a factor ∼5).

  • The addition of dynamics allows us to identify degeneracies between lens mass parameters and source kinematic properties. In particular, we note a tight correlation between the source centroid and the lens slope and Einstein radius.

  • Different kinematic data binning do not show, for this particular case, prominent differences on the lens mass parameters, probably due to the subdominant amount of data points with respect to the lensing constraints. Despite this, we do note a factor ∼1.1−1.2 improvement on some source parameters, such as Sérsic amplitude Ie, s, effective radius Reff, s, Sérsic index ns, source’s mass axis ratio q s piemd $ q^{\mathrm{piemd}}_{\mathrm{s}} $ and Einstein radius θ E , s piemd $ \theta^{\mathrm{piemd}}_{\mathrm{E,s}} $. The slight improvement of SN20 shows that, to tighten parameters’ constraints, it is better to obtain more spatially resolved elements, at the expense of higher S/N per bin.

  • When testing systematics, i.e. modelling the source mass with a mass-follows-light profile, we find that this biases the model parameters, in particular the predicted source inclination and the lens’s Einstein radius, slope and axis ratio. It also mildly underestimates the source vrms, leading to a higher χ dyn 2 $ \chi^2_{\rm dyn} $ for this model. The bias is particularly noticeable on the reconstructed surface mass density of the source galaxy which, for this particular system, is underestimated with respect to the data in the regions larger than 0.2 Reff. Since the region within 0.2 Reff corresponds to the higher signal-to-noise regions on the image plane, this shows that the surface mass density of the source is very sensitive to the choice of the model in regions that have lower signal-to-noise data. We conclude that the addition of the kinematic analysis can help better distinguish between models than the lensing analysis alone.

  • The analysis carried out with the mock data from the composite mass distribution for the lens mass shows similar results as the previous set of mock data. We still find improved constraints on the source centroid, and also on the dark matter (NFW profile) orientation and Einstein radius, when including the kinematic data compared to the lensing-only analysis. Furthermore, comparing the lensing and dynamics and lensing-only power-law model we find, consistent with the previous set of mock data, that the lens Einstein radius is better constrained by a factor 2 when including dynamics. Another interesting observation from this test is that if we model the composite lens mass with a single power-law profile, the combined lensing and dynamics analysis has a higher χ2 than the lensing only. This increased value comes from the misfit of the kinematic data, showing that the addition of these constraints can actually discern better between models.

Our analysis shows that the combination of the lensing and the dynamical analysis as complementary probes grants an interesting gain of information on both the source and lens galaxies, and it also allows us to break some of the degeneracies relative to each method. With the improvement in data quality brought by instruments as those on JWST and the Extremely Large Telescope (ELT), this method allows a detailed analysis of the source galaxies, giving us an extensive look into the high-redshift Universe.

1

Jampy and MGEfit, online available via https://www-astro.physics.ox.ac.uk/~mxc/software/

3

Average circularised, since we integrate all the mass distributions over circles.

Acknowledgments

We thank Michele Cappellari for making JamPy and MGEfit routines publicly available. GC thanks Stefan Schuldt and Elisa Ritondale for their help and the useful discussions. GC, AY and SHS thank the Max Planck Society for support through the Max Planck Research Group for SHS.

References

  1. Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bacon, R., Adam, G., Baranne, A., et al. 1995, A&AS, 113, 347 [NASA ADS] [Google Scholar]
  3. Bagnasco, G., Kolm, M., Ferruit, P., et al. 2007, in Cryogenic Optical Systems and Instruments XII, SPIE Conf. Ser., 6692, 66920M [CrossRef] [Google Scholar]
  4. Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681 [NASA ADS] [CrossRef] [Google Scholar]
  5. Barkana, R. 1998, ApJ, 502, 531 [NASA ADS] [CrossRef] [Google Scholar]
  6. Barnabè, M., & Koopmans, L. V. E. 2007, ApJ, 666, 726 [NASA ADS] [CrossRef] [Google Scholar]
  7. Barnabè, M., Czoske, O., Koopmans, L. V. E., Treu, T., & Bolton, A. S. 2011, MNRAS, 415, 2215 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barnabè, M., Dutton, A. A., Marshall, P. J., et al. 2012, MNRAS, 423, 1073 [NASA ADS] [CrossRef] [Google Scholar]
  9. Barth, A. J., Sarzi, M., Rix, H.-W., et al. 2001, ApJ, 555, 685 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bendinelli, O. 1991, ApJ, 366, 599 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bezanson, R., van Dokkum, P. G., Franx, M., et al. 2011, ApJ, 737, L31 [NASA ADS] [CrossRef] [Google Scholar]
  13. Binney, J., & Tremaine, S. 1988, Sky Telesc., 76, 45 [Google Scholar]
  14. Birkmann, S. M., Ferruit, P., Rawle, T., et al. 2016, in Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, SPIE Conf. Ser., 9904, 99040B [Google Scholar]
  15. Birrer, S., Amara, A., & Refregier, A. 2016, JCAP, 8, 020 [NASA ADS] [CrossRef] [Google Scholar]
  16. Birrer, S., Amara, A., & Refregier, A. 2017, JCAP, 5, 037 [CrossRef] [Google Scholar]
  17. Cappellari, M. 2002, MNRAS, 333, 400 [NASA ADS] [CrossRef] [Google Scholar]
  18. Cappellari, M. 2008, MNRAS, 390, 71 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cappellari, M. 2016, ARA&A, 54, 597 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345 [NASA ADS] [CrossRef] [Google Scholar]
  21. Cappellari, M., Verolme, E. K., van der Marel, R. P., et al. 2002, ApJ, 578, 787 [NASA ADS] [CrossRef] [Google Scholar]
  22. Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709 [NASA ADS] [CrossRef] [Google Scholar]
  23. Cappellari, M., Romanowsky, A. J., Brodie, J. P., et al. 2015, ApJ, 804, L21 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013, ApJ, 773, 149 [CrossRef] [Google Scholar]
  25. Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743 [Google Scholar]
  26. Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Daddi, E., Renzini, A., Pirzkal, N., et al. 2005, ApJ, 626, 680 [NASA ADS] [CrossRef] [Google Scholar]
  28. Dai, X., & Guerras, E. 2018, ApJ, 853, L27 [NASA ADS] [CrossRef] [Google Scholar]
  29. Dai, X., Kochanek, C. S., Chartas, G., et al. 2010, ApJ, 709, 278 [NASA ADS] [CrossRef] [Google Scholar]
  30. de Zeeuw, P. T., Bureau, M., Emsellem, E., et al. 2002, MNRAS, 329, 513 [NASA ADS] [CrossRef] [Google Scholar]
  31. Di Teodoro, E. M., Grillo, C., Fraternali, F., et al. 2018, MNRAS, 476, 804 [CrossRef] [Google Scholar]
  32. Ding, X., Treu, T., Suyu, S. H., et al. 2017, MNRAS, 472, 90 [CrossRef] [Google Scholar]
  33. Dutton, A. A., Brewer, B. J., Marshall, P. J., et al. 2011, MNRAS, 417, 1621 [NASA ADS] [CrossRef] [Google Scholar]
  34. Dye, S., Furlanetto, C., Swinbank, A. M., et al. 2015, MNRAS, 452, 2258 [NASA ADS] [CrossRef] [Google Scholar]
  35. Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285, 723 [NASA ADS] [Google Scholar]
  36. Emsellem, E., Dejonghe, H., & Bacon, R. 1999, MNRAS, 303, 495 [NASA ADS] [CrossRef] [Google Scholar]
  37. Emsellem, E., Cappellari, M., Peletier, R. F., et al. 2004, MNRAS, 352, 721 [NASA ADS] [CrossRef] [Google Scholar]
  38. Emsellem, E., Cappellari, M., Krajnović, D., et al. 2007, MNRAS, 379, 401 [Google Scholar]
  39. Falcón-Barroso, J., Lyubenova, M., van de Ven, G., et al. 2017, A&A, 597, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
  41. Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518 [NASA ADS] [CrossRef] [Google Scholar]
  42. Gavazzi, R., Treu, T., Rhodes, J. D., et al. 2007, ApJ, 667, 176 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13 [NASA ADS] [CrossRef] [Google Scholar]
  44. Girard, M., Dessauges-Zavadsky, M., Schaerer, D., et al. 2018, A&A, 613, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
  46. Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009a, ApJ, 695, 1577 [NASA ADS] [CrossRef] [Google Scholar]
  47. Gültekin, K., Cackett, E. M., Miller, J. M., et al. 2009b, ApJ, 706, 404 [NASA ADS] [CrossRef] [Google Scholar]
  48. Gültekin, K., Cackett, E. M., Miller, J. M., et al. 2012, ApJ, 749, 129 [NASA ADS] [CrossRef] [Google Scholar]
  49. Hogg, D. W., Blanton, M., Strateva, I., et al. 2002, AJ, 124, 646 [NASA ADS] [CrossRef] [Google Scholar]
  50. Humphrey, P. J., & Buote, D. A. 2010, MNRAS, 403, 2143 [NASA ADS] [CrossRef] [Google Scholar]
  51. Jeans, J. H. 1922, MNRAS, 82, 122 [NASA ADS] [CrossRef] [Google Scholar]
  52. Jones, T. A., Swinbank, A. M., Ellis, R. S., Richard, J., & Stark, D. P. 2010, MNRAS, 404, 1247 [NASA ADS] [Google Scholar]
  53. Kassiola, A., & Kovner, I. 1993, ApJ, 417, 450 [NASA ADS] [CrossRef] [Google Scholar]
  54. Koopmans, L. V. E., Treu, T., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 649, 599 [NASA ADS] [CrossRef] [Google Scholar]
  55. Koopmans, L. V. E., Bolton, A., Treu, T., et al. 2009, ApJ, 703, L51 [NASA ADS] [CrossRef] [Google Scholar]
  56. Kriek, M., van Dokkum, P. G., Franx, M., et al. 2006, ApJ, 649, L71 [NASA ADS] [CrossRef] [Google Scholar]
  57. Kriek, M., van der Wel, A., van Dokkum, P. G., Franx, M., & Illingworth, G. D. 2008, ApJ, 682, 896 [NASA ADS] [CrossRef] [Google Scholar]
  58. Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221 [NASA ADS] [CrossRef] [Google Scholar]
  59. Kronawitter, A., Saglia, R. P., Gerhard, O., & Bender, R. 2000, A&AS, 144, 53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Maller, A. H., Simard, L., Guhathakurta, P., et al. 2000, ApJ, 533, 194 [NASA ADS] [CrossRef] [Google Scholar]
  61. McConnell, N. J., Ma, C.-P., Murphy, J. D., et al. 2012, ApJ, 756, 179 [NASA ADS] [CrossRef] [Google Scholar]
  62. Monnet, G., Bacon, R., & Emsellem, E. 1992, A&A, 253, 366 [NASA ADS] [Google Scholar]
  63. Morgan, N. D., Kochanek, C. S., Falco, E. E., & Dai, X. 2006, in American Astronomical Society Meeting Abstracts, BAAS, 38, 927 [Google Scholar]
  64. Murphy, J. D., Gebhardt, K., & Adams, J. J. 2011, ApJ, 729, 129 [NASA ADS] [CrossRef] [Google Scholar]
  65. Naab, T., Oser, L., Emsellem, E., et al. 2014, MNRAS, 444, 3357 [NASA ADS] [CrossRef] [Google Scholar]
  66. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
  67. Newman, A. B., Smith, R. J., Conroy, C., Villaume, A., & van Dokkum, P. 2017, ApJ, 845, 157 [NASA ADS] [CrossRef] [Google Scholar]
  68. Newman, A. B., Belli, S., Ellis, R. S., & Patel, S. G. 2018, ApJ, 862, 126 [CrossRef] [Google Scholar]
  69. Oldham, L., Auger, M. W., Fassnacht, C. D., et al. 2017, MNRAS, 465, 3185 [CrossRef] [Google Scholar]
  70. Oser, L., Ostriker, J. P., Naab, T., Johansson, P. H., & Burkert, A. 2010, ApJ, 725, 2312 [NASA ADS] [CrossRef] [Google Scholar]
  71. Patrício, V., Richard, J., Carton, D., et al. 2018, MNRAS, 477, 18 [NASA ADS] [CrossRef] [Google Scholar]
  72. Rizzo, F., Vegetti, S., Fraternali, F., & Di Teodoro, E. 2018, MNRAS, 481, 5606 [CrossRef] [Google Scholar]
  73. Rybak, M., Vegetti, S., McKean, J. P., Andreani, P., & White, S. D. M. 2015, MNRAS, 453, L26 [NASA ADS] [CrossRef] [Google Scholar]
  74. Sarzi, M., Rix, H.-W., Shields, J. C., et al. 2001, ApJ, 550, 65 [NASA ADS] [CrossRef] [Google Scholar]
  75. Sluse, D., Surdej, J., Claeskens, J.-F., et al. 2003, A&A, 406, L43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Sluse, D., Claeskens, J. F., Hutsemékers, D., & Surdej, J. 2007, A&A, 468, 885 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Suyu, S. H., & Halkola, A. 2010, A&A, 524, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Suyu, S. H., Marshall, P. J., Hobson, M. P., & Blandford, R. D. 2006, MNRAS, 371, 983 [NASA ADS] [CrossRef] [Google Scholar]
  79. Suyu, S. H., Marshall, P. J., Blandford, R. D., et al. 2009, ApJ, 691, 277 [Google Scholar]
  80. Suyu, S. H., Hensel, S. W., McKean, J. P., et al. 2012, ApJ, 750, 10 [NASA ADS] [CrossRef] [Google Scholar]
  81. Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70 [NASA ADS] [CrossRef] [Google Scholar]
  82. Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35 [NASA ADS] [CrossRef] [Google Scholar]
  83. Swinbank, A. M., Dye, S., Nightingale, J. W., et al. 2015, ApJ, 806, L17 [Google Scholar]
  84. Szomoru, D., Franx, M., van Dokkum, P. G., et al. 2010, ApJ, 714, L244 [NASA ADS] [CrossRef] [Google Scholar]
  85. Szomoru, D., Franx, M., & van Dokkum, P. G. 2012, ApJ, 749, 121 [NASA ADS] [CrossRef] [Google Scholar]
  86. Thomas, J., Ma, C.-P., McConnell, N. J., et al. 2016, Nature, 532, 340 [NASA ADS] [CrossRef] [Google Scholar]
  87. Toft, S., van Dokkum, P., Franx, M., et al. 2005, ApJ, 624, L9 [NASA ADS] [CrossRef] [Google Scholar]
  88. Toft, S., van Dokkum, P., Franx, M., et al. 2007, ApJ, 671, 285 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  89. Toft, S., Gallazzi, A., Zirm, A., et al. 2012, ApJ, 754, 3 [NASA ADS] [CrossRef] [Google Scholar]
  90. Treu, T., & Koopmans, L. V. E. 2002, ApJ, 575, 87 [NASA ADS] [CrossRef] [Google Scholar]
  91. Treu, T., & Koopmans, L. V. E. 2004, ApJ, 611, 739 [NASA ADS] [CrossRef] [Google Scholar]
  92. Trujillo, I., Feulner, G., Goranova, Y., et al. 2006, MNRAS, 373, L36 [NASA ADS] [CrossRef] [Google Scholar]
  93. van de Sande, J., Kriek, M., Franx, M., et al. 2013, ApJ, 771, 85 [NASA ADS] [CrossRef] [Google Scholar]
  94. van de Ven, G., Falcón-Barroso, J., McDermid, R. M., et al. 2010, ApJ, 719, 1481 [NASA ADS] [CrossRef] [Google Scholar]
  95. van den Bosch, R. C. E., & de Zeeuw, P. T. 2010, MNRAS, 401, 1770 [NASA ADS] [CrossRef] [Google Scholar]
  96. van den Bosch, R. C. E., Gebhardt, K., Gültekin, K., et al. 2012, Nature, 491, 729 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  97. van der Wel, A., Holden, B. P., Zirm, A. W., et al. 2008, ApJ, 688, 48 [NASA ADS] [CrossRef] [Google Scholar]
  98. van der Wel, A., Rix, H.-W., Wuyts, S., et al. 2011, ApJ, 730, 38 [NASA ADS] [CrossRef] [Google Scholar]
  99. van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28 [NASA ADS] [CrossRef] [Google Scholar]
  100. van Dokkum, P. G., Franx, M., Kriek, M., et al. 2008, ApJ, 677, L5 [NASA ADS] [CrossRef] [Google Scholar]
  101. van Dokkum, P. G., Kriek, M., & Franx, M. 2009, Nature, 460, 717 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  102. van Dokkum, P. G., Whitaker, K. E., Brammer, G., et al. 2010, ApJ, 709, 1018 [NASA ADS] [CrossRef] [Google Scholar]
  103. Walsh, J. L., Barth, A. J., Ho, L. C., & Sarzi, M. 2013, ApJ, 770, 86 [NASA ADS] [CrossRef] [Google Scholar]
  104. Weijmans, A.-M., Cappellari, M., Bacon, R., et al. 2009, MNRAS, 398, 561 [NASA ADS] [CrossRef] [Google Scholar]
  105. Wellons, S., Torrey, P., Ma, C.-P., et al. 2016, MNRAS, 456, 1030 [NASA ADS] [CrossRef] [Google Scholar]
  106. Yıldırım, A., van den Bosch, R. C. E., van de Ven, G., et al. 2017, MNRAS, 468, 4216 [NASA ADS] [CrossRef] [Google Scholar]
  107. Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, MNRAS, 493, 4783 [CrossRef] [Google Scholar]
  108. Zirm, A. W., van der Wel, A., Franx, M., et al. 2007, ApJ, 656, 66 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Input lens mass parameters for RX J1131−1231 obtained by fitting both the quasar image positions as point sources as well as the full Einstein ring as an extended source (Yıldırım et al. 2020).

In the text
Table 2.

Input Sérsic light parameters of the mock lens and source galaxies of RX J1131−1231.

In the text
Table 3.

Input source mass parameters of the mock RX J1131−1231.

In the text
Table 4.

Lens and source parameter constraints for the different configurations for our first set of mock data (presented in Sect. 3.1).

In the text
Table 5.

Modelled lens and source parameter values for the different configurations of the composite mock data (presented in Sect. 3.2.4), with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin (shown in Fig. 11).

In the text

All Figures

thumbnail Fig. 1.

Simulated source stellar kinematics with GLaD. We produce an axisymmetric Jeans model of a mock source galaxy on the source plane (left panel) and we predict the lensed vrms profile in the image plane (right panel). The masked black region corresponds to regions in the image plane which have no correspondence on the source plane. The yellow cross marks the lens centroid position.
In the text
thumbnail Fig. 2.

Simulated strong lensing data with GLaD. We model the light of the source (left panel) and the mass and light of the lens and we predict how the surface brightness of the lensed source (central panel) and lensed source with lens (right panel) will look like. The units of the surface brightness are in counts.
In the text
thumbnail Fig. 3.

HST ACS image of RX J1131−1231 in the F814W filter. The lens galaxy G is surrounded by an Einstein ring and the four images of the quasar labelled with letters (A, B, C, D). The satellite in the northern part of the lens galaxy is not included in our mock data. The image resolution is 0.05″ pixel−1.
In the text
thumbnail Fig. 4.

Mock surface brightness image of RX J1131−1231 obtained, as explained in Sect. 3.1.1, using the parameters in Tables 1 and 2. The lens galaxy at the centre is surrounded by an Einstein ring and the four images of the host galaxy (quasar excluded). The image has a resolution of 0.05″ pixel−1. In the mock image we added simulated noise, as explained in Sect. 3.1.1.
In the text
thumbnail Fig. 5.

Mock binned lensed kinematic data of the source galaxy of RX J1131−1231. We scale the signal of the source galaxy to obtain a signal-to-noise ratio in the brightest pixel of 15, to mimic a realistic scenario (see Sect. 3.1.2 for further details). The data are binned to obtain a signal-to-noise in each bin of 20. We obtain a total of 45 bins. The pixel resolution of the kinematic map is 0.1″. The masked black region corresponds to regions in the image plane which have no correspondence on the source plane. The yellow cross marks the lens centroid position.
In the text
thumbnail Fig. 6.

Mock lensed kinematic data of the source galaxy of RX J1131−1231. We scale the signal of the source galaxy to obtain a signal-to-noise ratio in the brightest pixel of 15, to mimic a realistic scenario (see Sect. 3.1.2 for further details). The data are binned to obtain a signal-to-noise in each bin of 30. We obtain a total of 19 bins. The pixel resolution of the kinematic map is 0.1″. The masked black region corresponds to regions in the image plane which have no correspondence on the source plane. The yellow cross marks the lens centroid position.
In the text
thumbnail Fig. 7.

Joint 2D posterior probability distribution for the parameters of the power-law mock (presented in Sect. 3.1), with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin. The different contours in the 2D plots indicate, respectively, the 1σ, 2σ and 3σ CIs. Parameters shown are those where the improvement on the constraints coming from the combination of lensing and dynamics (LD in red contours) is more prominent as compared to lensing only (L in grey contours). These parameters are the source centroid xc, s, yc, s, the anisotropy β, the source and lens light axis ratio qs, ql, the lens Einstein radius θ E , ∞ , l spemd $ \theta^{\mathrm{spemd}}_{\mathrm{E,\infty,l}} $ and slope γ l spemd $ \gamma^{\mathrm{spemd}}_{\mathrm{l}} $, and the shear parameters γext, ∞ and ϕext. In the diagonal are shown the 1D histograms of the corresponding parameter on the x-axis. The input values are marked with a yellow cross.
In the text
thumbnail Fig. 8.

Joint 2D posterior probability distribution for the source mass parameters of the power-law mock data presented in Sect. 3.1, with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin. The different contours in the 2D plots indicate the different σ CIs. Parameters shown are those where the improvement on the constraints coming from the addition of lensing and dynamics (LD in red contours) is more prominent as compared to lensing only (L in grey contours). These parameters are the anisotropy β, the inclination i, the source mass axis ratio q s piemd $ q^{\mathrm{piemd}}_{\mathrm{s}} $, the source Einstein radius θ E , ∞ , s piemd $ \theta^{\mathrm{piemd}}_{\mathrm{E,\infty, s}} $. In the diagonal are shown the 1D histograms of the corresponding parameter on the x-axis. The input values are marked with a yellow cross.
In the text
thumbnail Fig. 9.

Mock lensed kinematic data of the source galaxy of RX J1131−1231 obtained using an isothermal mass profile (left panel), reconstructed kinematic map from our best-fit MFL model (central panel), and normalised residuals (right panel). The pixel resolution of the kinematic map is 0.1″. The masked black region is the region in the image plane which has no correspondence on the source plane. The yellow cross marks the lens centroid position. If we compare to the data, we see that our model predicted source’s vrms is slightly underestimated.
In the text
thumbnail Fig. 10.

Predicted averaged circularised surface mass density Σ(< R) (left) and mass enclosed M(< R) (right) within one effective radius Reff = 0.5″ (3.5 kpc) of the source galaxy, with the corresponding CI (shaded regions). Both the convergence and the enclosed mass of the source of the MFL model (see Sect. 3.2.3) are underestimated in the outskirts, as compared to the input and the LD model. This is also shown in the vrms value, which is slightly underestimated, as shown in Fig. 9.
In the text
thumbnail Fig. 11.

Mock lensed kinematic data of the source galaxy of RX J1131−1231 for the composite mass model (left panel), our prediction obtained by modelling the data with a power law profile (central panel) and the normalised residuals (right panel). The data are binned to obtain a signal-to-noise in each bin of 30, for a total of 22 bins. The pixel resolution of the kinematic maps is 0.1″. The masked black region is the region in the image plane which has no correspondence on the source plane. The yellow cross marks the lens centroid position.
In the text
thumbnail Fig. 12.

Joint 2D posterior probability distribution for the parameters of the composite mock (presented in Sect. 3.2.4), with kinematic data having a signal-to-noise of 15 in the brightest pixel and binned to have a signal-to-noise of 30 in each bin (shown in Fig. 11). The different contours in the 2D plots indicate, respectively, the 1σ, 2σ and 3σ ranges. Parameters shown are those where the improvement on the constraints coming from the combination of lensing and dynamics (LD in red contours) is more prominent as compared to lensing only (L in grey contours). These parameters are the source centroid xc, s, yc, s, the anisotropy β, the position angle of the lens dark matter profile θ l nfw $ \theta^{\mathrm{nfw}}_{\mathrm{l}} $, its Einstein radius θ E , ∞ , l nfw $ \theta^{\mathrm{nfw}}_{\mathrm{E,\infty,l}} $, and the shear parameters γext, ∞ and ϕext. In the diagonal are shown the 1D histograms of the corresponding parameter on the x-axis. The input values are marked with a yellow cross.
In the text
thumbnail Fig. 13.

Average surface mass density Σ ¯ ( R ) $ \bar{\Sigma}(R) $ for the LD composite and the L composite models (left panel), and for the power-law models LD power-law and L power-law (right panel) as compared to the input composite mock, plotted with the relative error bars. All models recover the input closely at the value of the Einstein radius (∼1.6″ for source at reshift 0.654). On the left panel, the L composite model gets peakier towards the centre, but is still able to recover the input within the error bars (light blue shaded region). This does not appear to be the case for the LD power-law and L power-law (right panel), which differ from the input for R <  0.3″. However, at larger radii, we cannot easily distinguish between power-law and composite models from the average surface mass density.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.