Gravitational Lensing and Dynamics (GLaD): combined analysis to unveil properties of high-redshift galaxies

The dynamical modelling of integral ﬁeld unit (IFU) stellar kinematics is a powerful tool to unveil the dynamical structure and mass build-up of galaxies in the local Universe, while gravitational lensing is nature’s cosmic telescope to explore the properties of galaxies beyond the local Universe. We present a new approach, which uniﬁes dynamical modelling of galaxies with the magniﬁcation power of strong gravitational lensing, to reconstruct the structural and dynamical properties of high-redshift galaxies. By means of axisymmetric Jeans modelling, we create a dynamical model of the source galaxy, assuming a surface brightness and surface mass density proﬁle. We then predict how the source’s surface brightness and kinematics would look when lensed by the foreground mass distribution and compare with the mock observed arcs of strong gravitational lensing systems. For demonstration purposes, we created and also analysed mock data of the strong lensing system RXJ1131 − 1231. By modelling both the lens and source, we recover the dynamical mass within the e ﬀ ective radius of strongly lensed high-redshift sources within 5% uncertainty, and we improve the constraints on the lens mass parameters by up to 50%. This machinery is particularly well-suited for future observations from large segmented-mirror telescopes, such as the James Webb Space Telescope, which will yield high sensitivity and angular-resolution IFU data for studies on distant and faint galaxies.


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. 2008van der Wel et al. , 2014van Dokkum et al. 2008;Szomoru et al. 2010Szomoru et al. , 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 Toft et al. 2007;Cimatti et al. 2008;van Dokkum et al. 2010), small Sérsic indices, and high stellar velocity dispersion 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 twophase 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 nonvisible (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. 2009aGültekin et al. ,b, 2012van 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. 2013Cappellari et al. , 2015Yı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 2002Barnabè & Koopmans 2007;Koopmans et al. 2009;van de Ven et al. 2010; Barnabè et al. 2011Barnabè et al. , 2012, 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. 2017Newman et al. , 2018Di 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 starforming 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 H 0 = 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.

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 2 LOS 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).

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 where (x, y) are the coordinates in the galaxy's plane (along the semi-major and semi-minor axes), R em is the elliptical mass radius, A135, page 2 of 16 e is the ellipticity e = 1−q 1+q with q the axis ratio, θ E,∞ is the lens Einstein radius for source at redshift infinity and r c 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 (x c , y c ), its axis ratio q, its position angle θ, its Einstein radius θ E,∞ , and its core radius r c .
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 where q is the axis ratio, r c 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 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 where ρ s is the characteristic overdensity and r s is the scale radius.
To represent both the lens and source light profiles, we use a Sérsic profile, whose intensity is given by where k is approximately 2n − 1 3 , R q = x 2 + y 2 /q 2 and n is the Sérsic index which, for most galaxies, spans values between 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 x c,s and y c,s , the axis ratio q, the position angle θ, the Sérsic amplitude I e , the effective radius R eff 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 where q l is the axis ratio, w t and w c are parameters of the profile (with w t > w c to keep L > 0). To represent the baryonic mass we scale the chameleon light profile by a constant mass-to-light ratio.

Dynamical modelling: axisymmetric Jeans modelling
The distribution function (DF) f (x, u) describing the positions x = (x, y, z) and velocities u = (v x , v y , v z ) 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 coor- where Φ is the gravitational potential, ν the luminosity density and 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 MGE 1 routine (Cappellari 2002(Cappellari , 2008Cappellari & 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 2 R = βv 2 z . Using these assumptions together with the boundary conditions that νv z = 0 for z → ∞, one can finally solve the Jeans equations for v 2 z and v 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. 1994Emsellem et al. , 1999Cappellari 2002). It consists of a series of expansions of the galaxy's image using 2D Gaussians functions where N is the number of Gaussians with luminosity L i , axis ratio q i , and dispersion σ i , (x, y) are a system of coordinates on the plane of the sky centred on the galaxy's nucleus, and (R,θ) the relative polar coordinates (x =R sin(θ − ψ i ), y = R cos(θ − ψ 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-ofsight to obtain the total observed first and second order velocity moment 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 v rms 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 (v rms ) simulated with GLaD in Fig. 1.

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 aŝ 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 , 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.

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 where P(η) is the prior probability on the model parameters, that we always assume to be uniform, and X obs collectively denotes our observables (surface brightness I obs and velocity map u obs rms ). To perform the joint analysis, we combine the lensing likelihood L lens with the dynamics likelihood L dyn such that the total likelihood function would be The lensing likelihood is where N pixels is the number of pixels in the image, I obs i is the observed surface brightness in a certain pixel i, I pred i (η) 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 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 where v obs rms,i is the observed value of Eq. (12) in each bin, v pred rms,i 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 A135, page 4 of 16  spaxel. When binning the v rms , the signal to noise in each bin is where N pixels is the number of pixels in each bin and I and σ here have the units of counts. The error on the v rms value in each bin is

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 (z s = 0.654, Sluse et al. 2007), lensed by an intervening giant elliptical galaxy at redshift z d = 0.295 (Sluse et al. 2003(Sluse et al. , 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  , 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.

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. 2013Suyu et al. , 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. 2006Koopmans et al. , 2009Gavazzi 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.

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.

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 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 z s = 0.654 hosting a black hole, whose mass we fixed to around 10 8 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, where D ds is the distance from our source at redshift z s = 0.654 to the artificial source at redshift infinity, and D s 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 v rms 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 where "bp" denotes brightest pixel. We solve Eq. (22)   Notes. The lens mass distribution is constituted by a power law profile (SPEMD) plus an external shear, whose parameters are, respectively, the centroid position x spemd c , y spemd c , the axis ratio q spemd , the position angle θ spemd , the Einstein radius for source at redshift infinity θ spemd E,∞ , the core radius r spemd core and the slope γ spemd for the power law, the external shear strength for source at redshift infinity γ ext,∞ and its orientation φ ext . 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 Notes. The Sérsic profile parameters are the centroid position x c and y c , the axis ratio q s , the position angle θ s , the Sérsic amplitude I e,s , the effective radius R eff,s and the Sérsic index n s . Notes. The source mass distribution is constituted by a pseudoisothermal elliptical mass distribution (PIEMD), whose parameters are the centroid position x piemd c , y piemd c , the axis ratio q piemd , the position angle θ piemd , the Einstein radius for source at redshift infinity θ piemd E,∞ and the core radius r piemd core . We also show the values of the anisotropy β and inclination i we use to mock the kinematic data.
in the final parameter constraints, we assume all the bins to contribute equally throughout the analysis.

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.
A135, page 7 of 16 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. 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.

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 χ 2 red ∼ 1. We note a tight anticorrelation between the Sérsic amplitude I e,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 I e,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 I e,s , effective radius R eff,s , Sérsic index n s , source's mass axis ratio q piemd s and Einstein radius θ piemd 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 v rms measurements.

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 I e,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 I e,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 χ 2 red ∼ 1.

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-tolight 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.
A135, page 8 of 16 Notes. The columns show the best-fitting values with the corresponding 1σ uncertainties for the joint lensing and dynamical models (LD) with different signal-to-noise ratios for the binning (SN20 and SN30), for the lensing-only model (L), and finally for models with a mass-follows-light profile for the source mass distribution (with SN30) to test for systematic uncertainties related to the source mass parameterisation.
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 χ 2 red than the other models (χ 2 red ∼ 1.07), showing that it does not fit as well to the data. In particular, the dynamics χ 2 dyn,red is ∼1.63, while the lensing χ 2 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 v rms 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 R eff , and underestimated otherwise, by up to a factor 2. We find that the size of the region within 0.2 R eff 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.

Composite models
As a final test, we assess the improvement of the constraints on the lens mass model, particularly focussing on the dark A135, page 9 of 16 . 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 x c,s , y c,s , the anisotropy β, the source and lens light axis ratio q s , q l , the lens Einstein radius θ spemd E,∞,l and slope γ spemd 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. 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-tonoise 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 powerlaw 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 χ 2 red ∼ 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 χ 2 red (χ 2 red = 1.03) than the lensing only L power-law (χ 2 red ∼ 1.02). However, the LD power-law has a higher misfit of the kinematic data χ 2 dyn,red ∼ 1.25. Moreover, in the LD powerlaw 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  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 piemd s , the source Einstein radius θ piemd E,∞,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. could still be inferred from the kinematics despite differences in the lens mass profile between power-law and composite profiles.
To show more in detail the differences among these models, we compare their average surface mass density 3 , namelȳ where Σ(R) is the surface mass density. TheΣ(<R) for the input, LD composite and L composite, is shown in the left panel of Fig. 13, while theΣ(<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). We note mild inconsistencies at R < 0.3 , where the LD power-law and L power-law models get peakier than the input.
Therefore we conclude that, for this system, the composite and power law models are not distinguishable in terms of Σ(<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.

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 available. 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 A135, page 11 of 16 3) are underestimated in the outskirts, as compared to the input and the LD model. This is also shown in the v rms value, which is slightly underestimated, as shown in Fig. 9.
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.
A135, page 12 of 16 Notes. The model presented are that including both the lensing and dynamics constraints (LD composite), that including only the lensing ones (L composite), and the two power-law models, with lensing and dynamics constraints (LD power-law) and with those from lensing only (L power-law). Parameters are presented with the 1σ uncertainties.
-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.   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 x c,s , y c,s , the anisotropy β, the position angle of the lens dark matter profile θ nfw l , its Einstein radius θ nfw E,∞,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.
Sérsic index n s , source's mass axis ratio q piemd s and Einstein radius θ piemd 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 incli-nation and the lens's Einstein radius, slope and axis ratio. It also mildly underestimates the source v rms , leading to a higher χ 2 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 R eff . Since the region within 0.2 R eff corresponds to the higher signal-to-noise regions on the image plane, this shows that Inpu composi e LD power⊙law L power⊙law Fig. 13. Average surface mass densityΣ(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.
the surface mass density of the source is very sensitive to the choice of the model in regions that have lower signal-tonoise 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.