Time-delay cosmographic forecasts with strong lensing and JWST stellar kinematics

TL;DR

This paper introduces a joint strong lensing and stellar dynamics framework using JWST-like data to precisely measure cosmological distances and the Hubble constant from individual lens systems, enhancing future cosmological constraints.

Contribution

It develops a novel method combining pixelated source reconstruction with axisymmetric Jeans equations, demonstrating high-precision $H_0$ and distance measurements with upcoming telescope data.

Findings

Time-delay distance constrained to 3.1% uncertainty.

Potential to measure $H_0$ with < 2.0% uncertainty from multiple lenses.

Distance to lens constrained with 2.4% precision.

Abstract

We present a joint strong lensing and stellar dynamical framework for future time-delay cosmography purposes. Based on a pixelated source reconstruction and the axisymmetric Jeans equations, we are capable of constraining cosmological distances and hence the current expansion rate of the Universe ($H_0$) to the few percent level per lens, when high signal-to-noise integral field unit (IFU) observations from the next generation of telescopes become available. For illustrating the power of this method, we mock up IFU stellar kinematic data of the prominent lens system RXJ1131-1231, given the specifications of the James Webb Space Telescope. Our analysis shows that the time-delay distance ($D_{\Delta t}$) can be constrained with 3.1% uncertainty at best, if future IFU stellar kinematics are included in the fit and if the set of candidate model parameterisations contains the true lens potential. These constraints would translate to a 3.2% precision measurement on $H_0$ in flat $\Lambda$CDM cosmology from the single lens RXJ1131-1231, and can be expected to yield an $H_0$ measure with < 2.0% uncertainty, if similar gains in precision can be reached for two additional lens systems. Moreover, the angular diameter distance ($D_\mathrm{d}$) to RXJ1131-1231 can be constrained with 2.4% precision, providing two distance measurements from a single lens system, which is extremely powerful to further constrain the matter density ($\Omega_\mathrm{m}$). The measurement accuracy of $D_\mathrm{d}$, however, is highly sensitive to any systematics in the measurement of the stellar kinematics. For both distance measurements, we strongly advise to probe a large set of physically motivated lens potentials in the future, to minimise the systematic errors associated with the lens mass parameterisation.

Figures

Equations

$$t(\vec{\theta },\vec{\beta })=\frac{(1+z_{\mathrm{d}})}{c}\frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{ds}}}\varphi (\vec{\theta },\vec{\beta })$$

$$\phi(\vec{\theta},\vec{\beta})=\Big{[}\frac{(\vec{\theta}\;-\;\vec{\beta})^{2}}{2}-\psi_{\rm L}({\vec{\theta}})\Big{]}$$

$$D_{\Delta t}\equiv (1+z_{\mathrm{d}})\frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{ds}}}.$$

$$t(\vec{\theta },\vec{\beta })\propto \frac{1}{H_0}\varphi (\vec{\theta },\vec{\beta }),$$

$$\Delta t_{ij}=t_{i}-t_{j}\propto\frac{1}{H_{0}}\>\Big{[}\phi(\vec{\theta_{i}},\vec{\beta})\;-\;\phi(\vec{\theta_{j}},\vec{\beta})\Big{]}$$

$${\psi }_{\mathrm{L},\lambda }(\vec{\theta })=\frac{\lambda }{2}|\vec{\theta }{|}^2+\vec{s}\vec{\theta }+c_0+(1-\lambda ){\psi }_{\mathrm{L}}(\vec{\theta }),$$

$${\nabla }^2{\psi }_{\mathrm{L}}=4\pi G{\rho }_{2D}=2\kappa ,$$

$$\kappa =\frac{\Sigma (D_{\mathrm{d}}\vec{\theta })}{{\Sigma }_{\mathrm{crit}}}$$

$${\Sigma }_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{d}}{D}_{\mathrm{ds}}}$$

$${\beta }_{\lambda }=(1-\lambda )\beta -\vec{s}$$

$${\kappa }_{\lambda }=\lambda +(1-\lambda )\kappa (\vec{\theta }).$$

$$D_{\Delta t,\lambda }=\frac{D_{\Delta t}^{\mathrm{model}}}{(1-\lambda )},$$

$$D_{\Delta t}=\frac{D_{\Delta t}^{\mathrm{model}}}{(1-{\kappa }_{\mathrm{ext}})},$$

$$\frac{\partial f}{\partial t}+\sum _{i=1}^3{v}_i\frac{\partial f}{\partial x_i}-\frac{\partial {\psi }_{\mathrm{D}}}{\partial x_i}\frac{\partial f}{\partial v_i}=0,$$

$$\frac{\nu \overline{v_R^2}-\nu \overline{v_{\varphi }^2}}{R}+\frac{\partial (\nu \overline{v_R^2})}{\partial R}+\frac{\partial (\nu \overline{v_R{v}_z})}{\partial z}=-\nu \frac{\partial {\psi }_{\mathrm{D}}}{\partial R}$$

$$\frac{\nu \overline{v_R{v}_z}}{R}+\frac{\partial (\nu \overline{v_z^2})}{\partial z}+\frac{\partial (\nu \overline{v_R{v}_z})}{\partial R}=-\nu \frac{\partial {\psi }_{\mathrm{D}}}{\partial z},$$

$$\frac{{\beta }_z\nu \overline{v_z^2}-\nu \overline{v_{\varphi }^2}}{R}+\frac{\partial ({\beta }_z\nu \overline{v_z^2})}{\partial R}=-\nu \frac{\partial {\psi }_{\mathrm{D}}}{\partial R}$$

$$\frac{\partial (\nu \overline{v_z^2})}{\partial z}=-\nu \frac{\partial {\psi }_{\mathrm{D}}}{\partial z}.$$

$$\begin{split}\overline{v_{\mathrm{LOS}}^{2}}=&\frac{1}{\mu(x^{\prime},y^{\prime})}\;\int_{-\infty}^{\infty}\;\nu\Big{[}(\overline{v_{R}^{2}}\sin^{2}\phi+\overline{v_{\phi}^{2}}\cos^{2}\phi)\sin^{2}i\\ &+\overline{v_{z}^{2}}\cos^{2}i-\overline{v_{R}v_{z}}\sin\phi\sin(2i)\Big{]}dz^{\prime}\;\equiv\;v^{2}+\sigma^{2}.\end{split}$$

$$\mu(x^{\prime},y^{\prime})=\sum_{j=1}^{N}\;\mu_{0,j}\;\exp\Big{[}-\frac{1}{2\sigma^{\prime 2}_{j}}\;(x^{\prime 2}+\frac{y^{\prime 2}}{q^{\prime 2}_{j}})\Big{]},$$

$$\nu(R,z)=\sum_{j=1}^{N}\frac{q^{\prime}_{j}\mu_{0,j}}{\sqrt{2\pi}\sigma^{\prime}_{j}q_{j}}\;\exp\Big{[}-\frac{1}{2\sigma^{\prime 2}_{j}}(R^{2}+\frac{z^{2}}{q^{2}_{j}})\Big{]},$$

$${\psi }_{\mathrm{D},j}(R,z)=-\frac{2G{\rm Y}\nu }{\sqrt{2\pi }{\sigma }_j^{\prime }}{\int }_0^1{\mathcal{F}}_j(u)du$$

$$\mathcal{F}_{j}(u)=\exp\Big{[}-\frac{u^{2}}{2\sigma^{\prime 2}_{j}}\Big{(}R^{2}+\frac{z^{2}}{\mathcal{Q}^{2}_{j}(u)}\Big{)}\Big{]}\frac{1}{\mathcal{Q}_{j}(u)}$$

$$\vec{\theta }=\vec{\beta }-\vec{\nabla }{\psi }_{\mathrm{L}}(\vec{\theta })=\vec{\beta }-\vec{\alpha }(\vec{\theta }),$$

$$d=\mathsf{f}s+n.$$

$$\begin{split}\chi_{\mathrm{L}}^{2}\propto\;\exp\Bigg{[}-\frac{1}{2}&\Big{[}\sum_{i}\frac{({d_{i}}-d_{i,{\rm m}})^{2}}{\sigma_{{\rm d},i}^{2}}+\sum_{j}\frac{(a_{j}-a_{j,{\rm m}})^{2}}{\sigma_{{\rm a},j}^{2}}\\ &+\sum_{k}\frac{(\Delta t_{k}-\Delta t_{k,{\rm m}})^{2}}{\sigma_{\Delta t,k}^{2}}\Big{]}\Bigg{]}.\end{split}$$

$$\kappa=\Big{[}\frac{\zeta^{2}+\zeta_{\rm c}^{2}}{E^{2}}\Big{]}^{\eta/2-1},$$

$$P_{\mathrm{L}}({\tau }_{\mathrm{L}}|d_{\mathrm{L}})\propto \stackrel{\text{likelihood}}{\overbrace{P_{\mathrm{L}}(d_{\mathrm{L}}|{\tau }_{\mathrm{L}})}}\underset{\text{prior}}{\underbrace{P_{\mathrm{L}}({\tau }_{\mathrm{L}})}},$$

$$P_{\mathrm{D}}({\tau }_{\mathrm{D}}|d_{\mathrm{D}})\propto P_{\mathrm{D}}(d_{\mathrm{D}}|{\tau }_{\mathrm{D}})P_{\mathrm{D}}({\tau }_{\mathrm{D}})$$

$$\chi_{\mathrm{D}}^{2}=\exp\Big{[}-\frac{1}{2}\sum_{l}\frac{\left(V_{\mathrm{rms},l}-\sqrt{\overline{v_{\mathrm{LOS},l,m}^{2}}}\right)^{2}}{\sigma^{2}_{V_{\mathrm{rms},l}}}\Big{]}.$$

Full text

1 Introduction

According to our standard cosmological model, we live in a flat, cold, dark matter and dark energy dominated Universe ($\Lambda$CDM). While little is yet known about the nature of dark matter and dark energy, our standard cosmological model - anchored mainly through measurements of anisotropies in the Cosmic Microwave Background (CMB; de Bernardis et al., 2000; Hanany et al., 2000; Bennett et al., 2003) - has been well established and provides an accurate description of e.g. the large scale structure formation and distribution (see; Springel et al., 2006, for a review) and the abundance of various elements in the Universe. As powerful as this model is, however, it has been constantly facing challenges. On small scales, the well-known "core-cusp" issue (Moore, 1994; McGaugh & de Blok, 1998) as well as the "Missing Satellites Problem" (Kauffmann et al., 1993; Klypin et al., 1999; Moore et al., 1999; Boylan-Kolchin et al., 2011) are stubbornly defying predictions from cosmological N-body simulations within the $\Lambda$CDM framework. Certainly, some of the disagreements can be attributed to observational effects (van den Bosch et al., 2000), underlying modelling assumptions (Evans et al., 2009) and our ignorance of the small-scale physics and baryonic feedback processes and interactions (Oh et al., 2011). But, while more recent studies claim to solve some of the small scale discrepancies, by thoroughly accounting for e.g. star formation, supernovae feedback (e.g. Zolotov et al., 2012; Madau et al., 2014; Tollet et al., 2016) and environmental impact such as ram pressure and tidal stripping (e.g. Brooks et al., 2013; Dutton et al., 2016; Sawala et al., 2016), similar discrepancies between theory and observations also arise at large scales. After making assumptions for a handful of parameters, such as spatial flatness and a constant dark energy equation of state of $-1$ (corresponding to the cosmological constant $\Lambda$), the standard cosmological model provides stringent constraints for the expansion rate of the Universe $H_{\mathrm{0}}$ through observations of the CMB (Planck Collaboration et al., 2018), which appears to be at odds with local measurements based on Cepheids and Type Ia Supernovae (Riess et al., 2018a, b, 2019). Especially the latter has been of particular interest lately. Given the significant 4.4$\sigma$ (i.e. $\sim$9.8%) discrepancy between the most recent measurements from both the Planck Collaboration and the Cepheid distance ladder, this result is either interpreted as corroborating evidence for a non-standard cosmological model (e.g. Lusso et al., 2019) or claimed to be part of unknown systematic effects that are not properly accounted for (e.g. Shanks et al., 2019). Naturally, relaxing our assumption about spatial flatness, a constant dark energy equation of state (i.e., not fixed to $-1$ that corresponds to the cosmological constant $\Lambda$) or increasing the number of relativistic species would allow us to reconcile both measurements. But, before such drastic conclusions are drawn, independent measurements of $H_{\mathrm{0}}$ should be carried out, if feasible, to test for possibly unknown systematics in any single method and to assess the need for physics beyond the standard model.

Time-delay cosmography (TDC; Refsdal, 1964) provides a methodologically independent tool for measuring $H_{\mathrm{0}}$ to the percent level (see recent reviews by, e.g., Treu & Marshall, 2016; Suyu et al., 2018). By means of a multiply imaged, time-variable background source and an accurate description of the foreground lens mass distribution, the time-delay distance ($D_{\Delta t}$) can be inferred, which is inversely proportional to the Hubble-Lemaître constant $H_{\mathrm{0}}$. The technique has long been plagued by poor time-delay measurements, invalid assumptions about the lens mass profile and systematic errors. However, it has been demonstrated that exhaustive studies of lensed quasars with exquisite light curves allow the measurement of $H_{\mathrm{0}}$ for a single system with an accuracy of $\sim$ 7% (Suyu et al., 2010, 2013). In addition, it was shown that TDC leads to tight constraints on other cosmological parameters, competing with those from contemporary Baryon Acoustic Peak studies, when each probe is combined with the CMB (Suyu et al., 2013).

In light of the aforementioned discrepancy between the current best cosmological probes, TDC has gained momentum and the H0LiCOW111$H_{\mathrm{0}}$ Lenses in COSMOGRAIL’s Wellspring; http://www.h0licow.org/ ($H_{\mathrm{0}}$ Lenses in COSMOGRAIL’s Wellspring) program has been initiated, which aims at measuring $H_{\mathrm{0}}$ with better than 3.5% precision and accuracy (Suyu et al., 2017). As part of these efforts, complementary data sets consisting of i) high-cadence and long-baseline monitoring of quasar light curves, mostly through COSMOGRAIL222COSmological MOnitoring of GRAvItational Lenses; https://cosmograil.epfl.ch, ii) high-spatially resolved photometric observations of the foreground lenses, lensed background quasars and quasar hosts, iii) wide-field photometric and spectroscopic observations of the lens’ environments and iv) stellar kinematic data of the foreground lenses have been obtained. Each of these ingredients are crucial to break the inherent modelling degeneracies, i.e. the mass-sheet degeneracy (MSD; Falco et al., 1985), in strong lensing studies and to reliably pin down the time-delay distance and hence the Hubble-Lemaître constant in a single system. As of now, H0LiCOW reports a 3.0% measurement of $H_{\mathrm{0}}$ in flat LCDM, based on a joint analysis of 6 strong lensing systems (Wong et al., 2019; Chen et al., 2019; Rusu et al., 2019; Birrer et al., 2019; Bonvin et al., 2017; Wong et al., 2017; Suyu et al., 2014, 2010). Yet, it is noteworthy that an $H_{\mathrm{0}}$ measurement which is comparable in precision with the best available probes (i.e. $\sim$ 2%) would still require the combination of almost a dozen lenses (Shajib et al., 2018). Accordingly, a measurement with 1% precision - a value that is considered as being highly beneficial for any Stage III and IV cosmological study to further constrain the dark energy equation of state (Weinberg et al., 2013) - would not be available until 40 such measurements have been carried out with similar precision.

Given these numbers and forecasts, a truly competitive TDC probe would greatly benefit from a much improved accuracy and precision for each lens study. In particular, three sources of uncertainty have been identified as the biggest contributors to the total time-delay distance error budget, i) the time delays, ii) the mass along the line of sight (LOS), and iii) the lens mass parameterisation. Assuming that future time delays can be measured to the percent level - based on long-baseline optical monitoring campaigns and new curve-shifting algorithms (Tewes et al., 2013a) -, any time-delay cosmological probe that aims at obtaining an accuracy and precision of 1-2 percent in the near future, will have to drastically improve their estimate of the external convergence ($\kappa_{\rm ext}$) associated with LOS mass distributions and lift the modelling degeneracy due to different lens mass parameterisations.

Interestingly, McCully et al. (2017) developed a new framework to model LOS mass distributions efficiently and quantified the environmental effects through realistic simulations of lens fields. By reconstructing the three-dimensional mass distribution of strong-lens sightlines, they obtain constraints on $\kappa_{\rm ext}$ that are consistent with those from statistical approaches of combining galaxy number density observations with N-body simulations (Hilbert et al., 2007, 2009; Collett et al., 2013; Suyu et al., 2014), but with a 4 times narrower distribution which yields much stronger priors on $\kappa_{\rm ext}$ (which affects $H_{\mathrm{0}}$ linearly). Progress in reducing the uncertainty in the lens mass parameterisation, on the other hand, has been moderate. Stellar kinematic data of the foreground lens are now commonly employed to break the MSD and to align the time-delay distance measurements when e.g. a Navarro-Frank-White (NFW; Navarro et al., 1996, 1997) or power law profile are adopted for the lens mass model. But, follow-up kinematic observations from adaptive-optics (AO) assisted ground-based facilities struggle to go beyond a single aperture-averaged velocity dispersion measurement, due to the difficulty in separating the bright quasar light from the foreground lens galaxy and the faintness of the lens itself. As a consequence, the final precision on $H_{\mathrm{0}}$ is currently limited to $\sim$ 7% from a single lens system. Moreover, time-delay studies employing kinematic data assume a spherical mass model for recovering the aperture-averaged stellar velocity dispersion, whereas the strong lensing mass model is of elliptical or even triaxial nature, and are thus not fully self-consistent. Whether this assumption introduces a bias in the inferred time-delay distances also remains to be seen.

To drastically improve the precision and accuracy of a single lens study, more flexible dynamical models along with both high-spatially resolved observations that map the 2D stellar kinematics in great detail and sufficient signal-to-noise (S/N) to reliably extract the kinematic moments across the entire field-of-view (FOV) are necessary. The next generation of telescopes - such as the James Webb Space Telescope (JWST), the Thirty Meter Telescope (TMT) and the European Extremely Large Telescope (E-ELT) - will provide the required improvement in sensitivity and resolution. The aim of this paper is to present a fully self-consistent, physically motivated modelling machinery for TDC, that will be capable of exploiting this data set to its full potential. Given the specifications of JWST, we will create mock stellar kinematics and, based on a joint strong lensing & stellar dynamical analysis, forecast the cosmological constraints from future space- and ground-based telescope observations.

The paper is organised as follows. In Section 2, we cover the strong lensing and stellar dynamical theory and formalism. Section 3 will be used to present the already available HST observations, time delays and mock future stellar kinematics of the prominent strong lens configuration RXJ1131$-$1231. We model the data of RXJ1131$-$1231 within a Bayesian framework in Section 4, show the probability density function (PDF) for its time-delay distance and lens distance ($D_{\mathrm{d}}$) and discuss possible sources of uncertainty. The inference of the cosmological parameters is carried out in Section 5 and finally followed by a summary in Section 6.

Throughout this paper, we adopt a standard cosmological model with $H_{\mathrm{0}}$ $=82.5$ km s*-1* Mpc*-1*, a matter density of $\Omega_{\mathrm{m}}=0.27$ and a dark energy density of $\Omega_{\mathrm{\Lambda}}=0.73$, where our particular choice for $H_{\mathrm{0}}$ $=82.5$ km s*-1* Mpc*-1*is driven by the time-delay distance measurements of RXJ1131$-$1231 in Suyu et al. (2014).

2 Theory

2.1 Historical context

Strong gravitational lensing and stellar dynamics are powerful tracers of their underlying gravitational potential. Being subject to the mass-sheet degeneracy (in lensing) and mass-anisotropy degeneracy (in stellar dynamics), however, it has been quickly realised that a combination of both would be capable of lifting their inherent modelling degeneracies, while providing even tighter constraints for the respective mass models of any given system (Grogin & Narayan, 1996; Romanowsky & Kochanek, 1999; Koopmans, 2004). Early implementations of joint strong lensing and stellar dynamical models can also be found in Koopmans & Treu (2002) and Treu & Koopmans (2002a), where strong lensing and stellar kinematic data have been utilised to infer the internal mass distribution and dark matter content of intermediate-redshift ($z\leq 1$) galaxies (Koopmans & Treu, 2003; Treu & Koopmans, 2004), culminating in the SLACS survey (Bolton et al., 2006). These early studies mainly relied on HST imaging and (aperture averaged) spectroscopic data within the effective radius, while adopting spherical Jeans models. Ultimately, this approach has also been commonly applied for cosmological purposes (Treu & Koopmans, 2002b; Koopmans et al., 2003), to break the mass profile degeneracies in the lensing-only models.

With the advent of integral field spectroscopy and further refinement in lens (e.g. Blandford et al., 2001; Warren & Dye, 2003; Koopmans, 2005) and dynamical modelling machineries (e.g. Cretton et al., 1999; Cappellari, 2008), the joint analysis has been expanded to cover valuable information from the extended Einstein rings and 2D kinematics, while employing pixelated source reconstruction and more sophisticated two-integral (Barnabè & Koopmans, 2007; Barnabè et al., 2009, 2011) and three-integral dynamical models (van de Ven et al., 2010; Barnabè et al., 2012). Whereas the earliest implementations treated the subject in an inconsistent manner, fully decoupling the strong lensing and stellar kinematic data by, e.g. adopting elliptical and spherical mass models respectively and making simplistic assumptions about the orbital anisotropy profile (Osipkov, 1979; Merritt, 1985a, b), later models have been capable to remedy most of these shortcomings by treating the subject within a fully self-consistent (van de Ven et al., 2010) and statistically meaningful (Barnabè et al., 2012) framework. Here, we build on the work by van de Ven et al. (2010) and Barnabè et al. (2012), by extending the machinery to include time-delay data for cosmological purposes. Moreover, even with the next generation of ground- and space-based stellar kinematic data, the likelihood functions will be swamped by the lensing information, which is why we make use of the Bayesian Information Criterion (BIC; Schwarz, 1978) as a statistical tool to further break degeneracies in future TDC studies.

Our joint strong lensing and stellar dynamical modelling machinery relies on a pixelated source fitting algorithm and the solutions of the Jeans equations in axisymmetric lens geometry, which are embedded in a Bayesian framework. For brevity, we refer the reader to Suyu et al. (2006, 2010, 2013) and Cappellari (2008), which cover in detail the theory and application of each to real observational data. Here, we confine ourselves to a description of the main formalisms, following the strong lensing and stellar dynamical framework developed and formulated in Schneider (1985); Schneider et al. (1992); Blandford & Narayan (1986) and Binney & Tremaine (1987).

2.2 Time-delay strong lensing

In any strong lens configuration with a time-variable background source, an excess time delay

[TABLE]

will be observed. Here, $\vec{\theta}$ is the angular image position, $\vec{\beta}$ the corresponding source position, $z_{\mathrm{d}}$ the lens redshift, $D_{\rm d}$, $D_{\rm s}$ & $D_{\rm ds}$ the angular diameter distance to the lens, the source and between the lens and source respectively, and

[TABLE]

the Fermat potential. The difference in the light propagation time at image position $\vec{\theta}$ with respect to the non-lensed case can therefore be attributed to the first and second term in the Fermat potential, which represent the geometric excess path length and the gravitational time-delay of the lens potential $\psi_{\rm L}(\vec{\theta})$ respectively, and a combination of cosmological distances which are generally referred to as the time-delay distance

[TABLE]

With $D_{\Delta t}$ being inversely proportional to $H_{\mathrm{0}}$, Eq.1 can be rewritten as

[TABLE]

i.e., the excess time delay can be used as a cosmological probe, if the form of the lens potential is sufficiently well known. However, since $t(\vec{\theta},\vec{\beta})$ itself cannot be measured, we rely on relative time delays

[TABLE]

between multiple images $i$ and $j$, with e.g. quadruply imaged systems naturally providing more constraints than doubly lensed sources.

Given the observables $\vec{\theta_{i}}$ and $\Delta t_{ij}$, the lens potential $\psi_{\mathrm{L}}$ and the source position $\vec{\beta}$ need to be modelled accurately to infer $H_{\mathrm{0}}$. A major drawback of this inference, though, is the MSD333A special case of the Source-Position Transformation (Schneider & Sluse, 2014; Unruh et al., 2017). (Falco et al., 1985). For illustration purposes, we assume a transformation of the lens potential $\psi_{\mathrm{L}}$ of the form

[TABLE]

where $\lambda$, $c_{0}$ and $\vec{s}$ are constant scalars and vectors respectively. Moreover, the projected matter density $\rho_{\mathrm{2D}}$ is related to its gravitational potential via Poisson’s equation

[TABLE]

where

[TABLE]

is the projected dimensionless surface mass density (SMD) and

[TABLE]

the critical SMD, that is used to discriminate between the weak ($\kappa\lesssim 1$) and strong lensing regime ($\kappa\gtrsim 1$). According to Eqs. 6 and 7, any transformation $\psi_{\mathrm{L},\lambda}(\vec{\theta})$ of $\psi_{\mathrm{L}}(\vec{\theta})$ will translate into a transformed source position

[TABLE]

and a dimensionless SMD of the form

[TABLE]

That is, as $\lambda$, $c_{0}$ and $\vec{s}$ only change the position and scaling of the source, which by itself is not directly observable, the above transformation essentially implies that any transformation of the lens potential (and hence of the dimensionless lens SMD) will be compensated by a corresponding scaling in the source plane, leaving many observables invariant under the transformation. Unfortunately, however, $D_{\Delta t}$ does not belong to this set of invariants. Being highly susceptible to the MSD, $D_{\Delta t}$ will also be scaled as follows

[TABLE]

with $D_{\Delta t}^{\rm model}$ being the model time-delay distance (without accounting for the mass-sheet-transformation parameter $\lambda$ in Eq. 11), and any cosmological inference based on strong lensing alone is therefore fundamentally limited by our ignorance of $\lambda$.

Since lensing is sensitive to all mass along the LOS, including small and large scale structures in the projected vicinity which can contribute to the SMD at the lens location, the MSD is inherently linked to this external convergence ($\kappa_{\mathrm{ext}}$). In fact, the MSD stems from the degeneracy between $\kappa_{\mathrm{ext}}$ and the normalisation of the lens potential (but see Schneider & Sluse, 2013, for a critical discussion). Consequently, Eq. 12 reduces to

[TABLE]

where $D_{\Delta t}$ is the true time-delay distance to the specific sightline of the lens, after accounting for the MSD. In contrast to $\lambda$, though, $\kappa_{\mathrm{ext}}$ has the benefit of not simply being an arbitrary scaling of the lens potential, but being observationally and/or numerically assessable via photometric and spectroscopic observations of the lens environment (Fassnacht et al., 2006; Momcheva et al., 2006; Rusu et al., 2017; Sluse et al., 2017) as well as ray-tracing methods through e.g. the Millennium Simulations (Hilbert et al., 2007, 2009; Greene et al., 2013; Collett et al., 2013) and weak lensing (Tihhonova et al., 2018). Whereas early attempts to quantify $\kappa_{\mathrm{ext}}$ have been only moderately successful, yielding external convergences that can affect the final measurement of $H_{\mathrm{0}}$ by 5% and more, more recent studies indicate that the distribution and impact of $\kappa_{\mathrm{ext}}$ can be drastically reduced when e.g. sightlines are not significantly overdense (Greene et al., 2013; Rusu et al., 2017) or individual lens fields are modelled (McCully et al., 2017).

Nonetheless, other means are needed to effectively break the MSD, and to reliably measure $D_{\Delta t}$. This is particularly evident from Eq. 4 and 7. Assuming, for instance, a simple power law profile for the 3D density distribution (i.e. $\rho_{\mathrm{3D}}(\vec{r})\propto r^{-\gamma}$), the lens potential becomes $\psi_{\mathrm{L}}\propto r^{2-\gamma}$. Any uncertainty in the slope of the mass density will thus translate into an uncertainty on the inferred time-delay distance ($D_{\Delta t}\propto\frac{1}{\gamma-1}$). To constrain $\gamma$, methods have been developed to make use of the spatially extended lensed images of the source galaxy, i.e. the Active Galactic Nuclei (AGN) host galaxy in the case of time-delay lenses (e.g. Warren & Dye, 2003; Suyu et al., 2006; Dye et al., 2008; Birrer et al., 2015). In this work, we follow the work of Suyu et al. (2006) and Suyu et al. (2013) for the lens modelling by describing the AGN host galaxy surface brightness on a grid of pixels, and the lens mass distribution with parameterised profiles. In addition, stellar kinematics of the lens galaxy are employed, which provide an independent assessment of the lens potential at different radii, further breaking lens mass model degeneracies for constraining $D_{\Delta t}$.

2.3 Axisymmetric Jeans modelling

The dynamical state of a system of particles is fully described by its distribution function (DF) $f(\vec{x},\vec{v})\geq 0$, with particle positions $\vec{x}$ and velocities $\vec{v}$. In the case that these particles are collisionless, interact purely via gravitational forces and are embedded in background potential $\psi_{\rm D}$ that is smooth in time and space, the time evolution of the DF is subject to the Collisionless Boltzmann Equation (CBE) (Binney & Tremaine, 1987)

[TABLE]

which basically postulates a conservation of the phase-space density. Yet, as the phase-space distribution is not accessible for objects beyond our Galaxy, where only bulk motions and positions of stars along specific LOSs are observed, the CBE is impractical for real observational purposes. In fact, any real world application of the CBE would require a more practical formalism, incorporating kinematic moments which are more easily measurable via line profile shifts and widths. This can be achieved by multiplying the CBE with powers of the velocity moment and subsequent integration over velocity space. Further, rewriting the CBE in terms of the cylindrical coordinate system ($R,z,\phi$) and under the assumption of axial symmetry, we obtain the two Jeans equations (Jeans, 1922; Binney & Tremaine, 1987)

[TABLE]

[TABLE]

where $\nu(\vec{x})=\int fd^{3}\vec{v}$ is the zeroth velocity moment and tracer density of the gravitational potential $\psi_{\mathrm{D}}$444$\psi_{\mathrm{D}}$ is the 3D gravitational potential, in contrast to $\psi_{\mathrm{L}}$ which denotes the 2D lens potential..

Given the four unknown second-order velocity moments $\overline{v_{i}^{2}}=\frac{1}{\nu}\int fv_{i}^{2}d^{3}v_{i}$ and $\overline{v_{i}v_{j}}=\int v_{i}v_{j}fd^{3}v$, the Jeans equations do not have a unique solution. In practice, assumptions about the shape and alignment of the velocity ellipsoid are made to simplify Eq. 15 and 16. These usually include the alignment of the velocity ellipsoid with the cylindrical coordinate system (i.e. $\overline{v_{R}v_{z}}=0$) and a flattening in the meridional plane, i.e. $\beta_{z}=1-{\overline{v_{z}^{2}}}/{\overline{v_{R}^{2}}}$ (Binney, 1980; Binney & Mamon, 1982; Binney & Tremaine, 1987), which yield the more commonly seen form (Cappellari, 2008)

[TABLE]

[TABLE]

These equations now link a mass and tracer density to three intrinsic second-order velocity moments which, in turn, can be used to obtain a projected second-order velocity moment along the LOS

[TABLE]

Here, $x^{\prime}$ and $y^{\prime}$ are the cartesian coordinates on the plane of the sky, $z^{\prime}$ the coordinate along the LOS, $i$ the inclination angle, $\mu$ the observed surface brightness (SB) - in contrast to $\nu$, which represents the (deprojected) intrinsic luminosity density - $\cos\phi=x/R$ (where $x$ and $R^{2}=x^{2}+y^{2}$ denote the intrinsic coordinate axis and cylindrical radius) and $v$ and $\sigma$ the observed mean LOS velocity and velocity dispersion (with $V_{\mathrm{rms}}=\sqrt{v^{2}+\sigma^{2}}$).

The assumptions for the shape and alignment of the velocity ellipsoid have been found to be a good description of the internal dynamical structure of fast-rotating (Emsellem et al., 2007) early-type galaxies (ETGs; Cappellari et al., 2007). Hence, the axisymmetric Jeans equations provide a decent fit to the observed kinematics and are usually in agreement with constraints that are obtained via more sophisticated orbit-based dynamical models. However, the Jeans equations do not make use of the higher order kinematic moments, which contain valuable information regarding the intrinsic shapes of galaxies. This is most prominent for massive, slow-rotating and pressure supported ETGs, where generally worse fits are obtained as the assumption of axial symmetry also breaks down (Li et al., 2018). In principle, the Jeans equations can be extended to triaxial systems, consisting of three equations and six second-order moments, but the set of solutions still contains unphysical DFs ($f<0$) (van de Ven et al., 2003). We proceed with these caveats in mind and note that the validity of the axisymmetric Jeans equations needs to be evaluated on a case by case basis, especially for the purposes of precision cosmology, where we aim to constrain the mass profile and hence time-delay distances to the percent level.

When constructing axisymmetric Jeans models, the intrinsic luminosity density is obtained by deprojecting the observed SB distribution. For this, we make use of a Multi-Gaussian Expansion (MGE) (Monnet et al., 1992; Emsellem et al., 1994; Cappellari, 2002)555For the dynamical modelling part, we make use of publicly available Python implementations of JAM (jampy; Cappellari, 2008) and MGEfit (mgefit; Cappellari, 2002), online available via https://www-astro.physics.ox.ac.uk/~mxc/software/.. In brief, the distribution is parameterised by a set of two-dimensional Gaussians, such that the SB can be written as

[TABLE]

where $\mu_{0}$ is the peak SB, $\sigma^{\prime}$ the dispersion along the projected major axis and $q^{\prime}$ the apparent flattening of each Gaussian. In an oblate axisymmetric case, the inclination $i$ is the only free viewing angle required to perform the deprojection. The deprojection is not unique (Rybicki, 1987), unless the galaxy is viewed edge-on, but konus densities which project to zero SB have been found to be of little effect for SB distributions of realistic (elliptical) galaxies (van den Bosch, 1997). If a deprojectable (i.e. $\cos^{2}i<q^{\prime 2}_{\rm min}$, with $q_{\rm min}$ being the axis ratio of the flattest Gaussian in the fit) inclination $i$ has been chosen, the intrinsic luminosity density in cylindrical coordinates reads as

[TABLE]

where $\sigma_{j}=\sigma^{\prime}_{j}$ and $q_{j}=\frac{\sqrt{q_{j}^{{}^{\prime}2}-\cos^{2}i}}{\sin i}$ now denote the intrinsic dispersion and flattening of the Gaussians. Simple (mass-follows-light) models can easily be constructed by linking the tracer density $\nu$ to the mass density $\rho$ via a mass-to-light ratio (M/L) for the individual Gaussians ($\Upsilon$j). The MGE is particularly handy here, as the gravitational potential can then be obtained by means of a simple, one-dimensional integral

[TABLE]

with

[TABLE]

and $\mathcal{Q}^{2}_{j}(u)=1-(1-q^{2}_{j})\;u^{2}$.

2.4 Joint formalism & Bayesian Inference

We start our joint formalism with the strong lensing part. Given the image positions $\vec{\theta}$ of the AGN and AGN host, a lens potential $\psi_{\rm L}$ will be adopted, which relates the AGN and AGN host source positions to those in the image plane via the lens equation

[TABLE]

where $\vec{\alpha}(\vec{\theta)}$ is the scaled deflection angle. We describe the source intensity distribution on a grid of pixels with values $\mathbb{s}$ (vector with dimension $N_{\rm s}$, the number of source pixels), which is related to the observed intensity value of the image plane $\mathbb{d}$ via

[TABLE]

Here, $\mathbb{d}$ is a vector with length $N_{\rm d}$ (the number of image pixels) and $\mathbb{n}$ the noise in the data. $\mathsf{f}$ represents the lensing operator (a matrix of dimension $N_{\rm d}\times N_{\rm s}$), which contains information regarding the lens potential and observational effects - such as telescope point spread function (PSF) - and which is constrained by the extended image positions and intensities as well as the relative time delays between the individual AGN images. In general, lens potentials for which the deflection angles $\vec{\alpha}(\vec{\theta})$ (and hence the lensing operator $\mathsf{f}$) can be obtained analytically and/or with only moderate numerical effort are adopted, and the goodness-of-fit for a particular model is defined as

[TABLE]

In the above equation, $d_{i}$ represents the image pixel intensities, $d_{i,{\rm m}}$ the modelled image pixel intensities and $\sigma_{{\rm d},i}$ the corresponding pixel uncertainties. Whereas the first term in Eq. 26 represents the fit to the image intensity distribution, the second and third terms account for the $\chi^{2}$ contribution from fitting in particular to the AGN positions ($a_{j}$) and their relative time delays ($\Delta t_{k}$). A best-fitting model is usually quickly obtained by minimising the cost function in Eq. 26, when the parameters of interest are few. However, since we are interested in inferring credible confidence intervals for all parameters of interest, we perform an analysis within the framework of Bayesian statistics.

For simplicity, let us assume that the strong lens configuration is well parameterised by a softened power law elliptical mass distribution (SPEMD; Barkana, 1998) with SMD

[TABLE]

with $\zeta^{2}=x^{\prime 2}+y^{\prime 2}/q^{2}$, where $E$ is a normalisation factor, $\eta$ the power law index, $\zeta_{\rm c}$ the core radius and $q$ the observed flattening with the $x^{\prime}$-axis being aligned with the projected galaxy major axis. According to Bayes’ theorem, the posterior PDF for this lensing-only model - with the set of parameters $\tau_{\mathrm{L}}=\{E,\eta,\zeta_{\rm c},q,\omega\}$ and data sets $d_{\mathrm{L}}=\{{d_{i},a_{j},\Delta t_{k}}\}$ - is given by

[TABLE]

where the log likelihood corresponds to Eq. 26 ($\log P_{\rm L}\propto-\chi^{2}_{\rm L}/2$, after marginalising over the source intensity pixel parameters $\mathbb{s}$; Suyu & Halkola, 2010) and $\omega$ comprises a set of remaining variables, such as $D_{\Delta t}^{\rm model}$, $\kappa_{\mathrm{ext}}$ and the position angle (PA) of the projected SPEMD on the plane of the sky (we measure counter-clockwise from the $x^{\prime}$-axis to the $y^{\prime}$-axis). Once (non-)informative priors for $\tau_{\mathrm{L}}$ have been chosen, the marginalised posterior PDF for a parameter of interest can be obtained by integrating the joint (lensing-only) PDF over all nuisance parameters. Similarly, the dynamics-only posterior PDF

[TABLE]

can easily be obtained by means of the dynamics-only likelihood

[TABLE]

Note that $\tau_{\mathrm{L}}\neq\tau_{\mathrm{D}}$. Unlike the lensing data, the dynamics are insensitive to e.g. $\kappa_{\mathrm{ext}}$, while depending explicitly on additional parameters, such as $\beta_{z}$ and $D_{\mathrm{d}}$. In the case of a joint lensing & dynamics model, the joint prior is simply a union of both $\tau_{\mathrm{LD}}=\{\tau_{\mathrm{L}},\tau_{\mathrm{D}}\}$ and the lensing & dynamics posterior PDF reads as

[TABLE]

given the independence of both data sets.

The exploration of the parameter space and adequate sampling of the joint lensing & dynamics posterior PDF is carried out by means of the affine-invariant ensemble sampler emcee (Foreman-Mackey et al., 2013). Starting points are obtained by first carrying out a pre-annealing process, in order to avoid low probability modes of the multi-parameter space, and initialising the walkers of the sampler such that they sample well the prior probability distribution. At each step, i.e. at each parameter combination that is probed by the walkers, the lensing likelihood in Eq. 26 will be evaluated. Secondary products of this evaluation are the dimensionless SMD (Eq. 8) and SB distribution, which are then transformed into a physical mass and luminosity density profile, before being parameterised by a MGE (Eq. 20 and 21) to allow for a straightforward, analytical calculation of the lens potential according to Eq. 22. Here, the critical SMD in Eq. 9 has to be expressed in terms of $D_{\Delta t}$, $D_{\mathrm{d}}$ and the lens redshift $z_{\mathrm{d}}$

[TABLE]

The SMD profile (as obtained from the lensing part) is multiplied with ($1-\kappa_{\mathrm{ext}}$), in order to take into account any contribution from the external convergence (see Eq. 11). As a consequence, $\kappa_{\mathrm{ext}}$ cancels out in the mass density profile that is used in Eq. 17 and 18, and we are thus insensitive to $\kappa_{\mathrm{ext}}$ when using the kinematic data. Note, however, that the absolute scaling of the lens potential is fixed and the MSD broken when stellar kinematics are included, which provide an independent measurement of the lens potential. Moreover, as $D_{\mathrm{d}}$ itself is insensitive to the external convergence along the LOS when inferred from kinematic and lensing observations (Jee et al., 2015, and also noted above), we have

[TABLE]

After deprojection and adoption of an anisotropy parameter $\beta_{z}$, the likelihood function in Eq. 30 can finally be evaluated. In combination with $\tau_{\mathrm{LD}}$, this yields the joint lensing & dynamics posterior PDF in Eq. 31, where the marginalised distributions can then be visualised by histograms with the most probable model and the 1$\sigma$ uncertainties being approximated by the median and 16th and 84th percentiles of the distribution.

2.5 Bayesian Information Criterion

A large, flexible and physically motivated set of light and mass parameterisations is generally utilised to accurately model the foreground lens. This approach is largely motivated by our ignorance of the true underlying gravitational potential, which is tightly linked to the excess time-delay (Eq. 5). Even if lens galaxies are found to be well approximated by power law density distributions (Koopmans et al., 2006), the inferred time-delay distances can differ significantly when compared to density distributions that follow more closely e.g. a NFW profile. The discrepancy in $D_{\Delta t}$ between different lens mass models can be partially alleviated by including stellar kinematic information. However, the vast amount of data from both lensing and future IFU stellar kinematics should enable us to perform a model selection, by means of the differences in their likelihood functions. To this end we will make use of the BIC (Schwarz, 1978), which is an approximation to the Bayesian evidence666In practice, the BIC approximation performs well, despite an asymptotic error of O(1). In cases where the data is of lower quality (i.e. with little constraining power) the impact of the error can still be mitigated by choosing appropriate priors, such that it becomes O($n^{-1/2}$), where $n$ is the number of data points (Kass & Wasserman, 1995). (see e.g. Raftery, 1995, for a thorough derivation) via

[TABLE]

where $\mathcal{M}$ is a hypothesis (in our case a model with lensing and dynamical parameters $\tau_{\mathrm{LD}}$). Following Bayes’ theorem, the ratio of posterior probabilities of two competing hypothesis $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are given by

[TABLE]

Accordingly, an approximation to the Bayes factor can be obtained by means of Eq. 34 and 35. Here, in particular, we make use of the BIC differences ($\Delta\mathrm{BIC}=\mathrm{BIC_{1}}-\mathrm{BIC_{2}}$) between model $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ with model parameters $\tau_{\mathrm{LD,1}}$ and $\tau_{\mathrm{LD,2}}$, which yields

[TABLE]

That is, the approximation to the Bayes factor can readily be utilised to assess the posterior odds, if the prior odds are assumed to be equally probable. Given a larger set of models $\mathcal{M}_{j}$ with BIC values $\mathrm{BIC}_{j}$, the posterior probability for a particular model $i\in j$ follows as

[TABLE]

which in turn can be used as a weighting scheme for model averaging purposes, since according to the law of total probability, $\sum\limits_{1}^{j}P_{\mathrm{LD}}(\mathcal{M}_{j}|d_{\mathrm{LD}})=1$.

The BIC is given by

[TABLE]

In the above equation, $k$ represents the number of model parameters, $n$ the number of data points and $\hat{L}$ the corresponding model maximum likelihood. The BIC discerns between candidate models by penalising models of increased complexity (i.e. with higher degrees freedom). The model maximum likelihood is given by $P_{\mathrm{L}}(d_{\mathrm{L}}|\tau_{\mathrm{L}})$, $P_{\mathrm{D}}(d_{\mathrm{D}}|\tau_{\mathrm{D}})$ or $P_{\mathrm{LD}}(d_{\mathrm{LD}}|\tau_{\mathrm{LD}})$, depending on which data set is fitted. Yet, owing to the MSD and the minuscule differences in the likelihood function of different mass parameterisations when modelling RXJ1131$-$1231 (Suyu et al., 2014), the lensing-only likelihood cannot be used to discriminate between models. Even in a joint model, the likelihood would easily be dominated by the large number of pixel intensities that are fitted in the lensing part, which is why we will use the differences in the goodness of fit of the IFU kinematics from our joint strong lensing & dynamics run, to perform a proper weighting of our models according to the BIC. We follow the approach of Birrer et al. (2019), where a weighting scheme with respect to the minimal BIC is defined as

[TABLE]

This weighting scheme follows Eq. 37, after accounting for the fact that the calculation of the denominator is cumbersome in most realistic scenarios. Especially for modelling time-delay lenses, where we rely on a relatively small subset of models, which are physically motivated by e.g. hydrodynamical numerical simulations (Navarro et al., 1996, 1997). Given its purpose of being a normalising factor, we approximate the weighting via Eq. 39, while making sure that $f_{\mathrm{BIC}}$ is still bound by 1.

For a given model including the lens mass/light distribution, PSF, AGN light, and AGN host galaxy surface brightness, the BIC value of this model could be computed to rank it relative to other models. We are particularly interested in comparing the lens mass parametrisation, and thus the changes in BIC due to different mass parameterisation. However, the BIC depends on also other modelling choices/parameters, especially the number of surface brightness pixels used to describe the AGN host galaxy, which introduces an uncertainty on the BIC (see Suyu et al., 2013, where source pixelisation dominates the uncertainties in the BIC for a given form of lens mass parametrisation). Given finite computing resources and thus a finite number of source intensity grids that we could explore, we quantify the uncertainty on the BIC due to the source grid pixelisation effect by comparing the BIC values of a range of source grids and estimating the scatter $\sigma_{\rm BIC}$. To account for the uncertainty in the BIC in weighting models, we follow Birrer et al. (2019) and convolve $f_{\mathrm{BIC}}$ in Eq. 39 with a Gaussian of variance $\sigma_{\rm BIC}^{2}$, obtaining the new weights as

[TABLE]

where

[TABLE]

Carrying out the convolution integral, we find an analytic expression as follows

[TABLE]

where the Erf and Erfc functions are defined as

[TABLE]

and

[TABLE]

3 Data

In this section, we present real and mock observations of RXJ1131$-$1231, which are used to infer future cosmological constraints by means of our joint strong lensing & stellar dynamical analysis. We focus on RXJ1131$-$1231, in particular, given that it i) is the brightest lens galaxy among the H0LiCOW base sample, ii) has the most precise time-delay measurements, with only 1.3% uncertainty in the longest time delay, and iii) has plenty of ancillary data, which make it the most promising candidate for future integral field unit (IFU) observations.

Our data consist of HST imaging, precise time delay measurements and mock JWST stellar kinematics. Rather than mocking up JWST imaging along with JWST IFU data, we choose to rely on the literature HST data instead. Even if JWST will be able to acquire comparably high S/N observations with much shorter exposure times, JWST NIR imaging will provide only a marginal improvement over HST, given its slightly smaller nominal pixel size of 0.03″/pixel and expected PSF FWHM of $\sim$ 2 pixels. More importantly, both HST and JWST are capable of obtaining spatially-resolved imaging with sufficient S/N for bright lenses, as in the case of RXJ1131$-$1231, to reach the required precision of 0.02 for the total mass density power law slope, such that the uncertainty on the Fermat potential is already subdominant with respect to other sources of uncertainty (Meng et al., 2015). Mocking up ground-based AO imaging from future 30-40m class telescopes would have been a viable alternative for assessing the overall improvements from both, next generation imaging and spectroscopy. However, the move from HST to ground-based AO observations introduces other sources of uncertainty, such as the recovery of the complex PSF; given that current state-of-the-art analysis of AO images from 8-10m class telescopes yield distance constraints comparable to those from HST (Chen et al., 2019), we adopt HST imaging for our current study.

3.1 Imaging

HST Advanced Camera for Surveys (ACS) data of RXJ1131$-$1231 have been obtained as part of programme GO:9744 (PI:Kochanek). The data set comprises imaging in the F814W and F555W filters, where five exposures each have been taken with a total integration time of 1980s. For the analysis, we give preference to the F814W imaging, given the fact that the stellar mass-to-light ratios become a weaker function of the underlying stellar populations (Bell & de Jong, 2001; Cole et al., 2001). Moreover, the F814W data shows a clearer separation between the AGN and the spatially extended Einstein ring, whereas the F555W filter is more difficult to model due to diffraction spikes extending into the lensed arcs (but see Birrer et al. (2016) for a joint modelling of both bands).

The reduction and combination of the imaging data is performed via the standard Multidrizzle pipeline (Fruchter & Hook, 2002), with charge transfer inefficiencies properly taken into account by empirically tracing back the charge-coupled device (CCD) detector readout mechanism and thus the initial charge distribution (Anderson & Bedin, 2010). The images are sky subtracted, corrected for geometric and photometric distortions and cosmic ray cleaned, before drizzled onto a final science frame with 0.05″/pixel resolution. Flux uncertainties for each pixel are obtained by adding in quadrature Poisson noise from the source and background noise from the sky and detector readout.

The final science frame is displayed in Fig. 1, where the centrally located galaxy lenses the background AGN into a quadruple lens configuration (A, B, C & D). The background quasar host is a spiral galaxy (Claeskens et al., 2006), which forms the extended Einstein ring. Discovered by Sluse et al. (2003), spectroscopic measurements of the foreground lens and background source yield a redshift of $z_{\mathrm{d}}=0.295$ and $z_{\mathrm{s}}=0.654$ (Sluse et al., 2003, 2007), respectively. The foreground lens is further accompanied by a satellite galaxy (S), which is assumed to be a dwarf elliptical (Claeskens et al., 2006).

3.2 Time delays

Time-delay measurements of RXJ1131$-$1231 have been carried out by means of a dedicated optical monitoring campaign within the COSMOGRAIL program. Based on high-cadence (3 days) long-baseline (9+ years and more than 700 epochs) observations with meter-class telescopes, and new curve-shifting techniques, Tewes et al. (2013a) report time-delays relative to image B of $\Delta t_{\rm AB}=0.7\pm 1.2$ days, $\Delta t_{\rm CB}=-0.4\pm 1.5$ days and $\Delta t_{\rm DB}=91.4\pm 1.2$ days, with systematic errors already taken into account in the uncertainty estimates. In general, the long-baseline measurements result in time delays with $\sim$ 3% precision (Tewes et al., 2013b; Liao et al., 2015; Bonvin et al., 2016), and microlensing shifts in the time-delays (Tie & Kochanek, 2018) are found to be negligible, given the long time delay (Chen et al., 2018).

Given the above measurements, RXJ1131$-$1231 is particularly suitable for TDC purposes. The background AGN is not only quadruply imaged, providing three independent time-delay constraints, but the longest time delay yields an uncertainty of only 1.3% and forms a comparably low floor for any time-delay distance measurement. While percent level precision of $\Delta t$ is a necessary condition for any TDC probe that aims to measure $H_{\mathrm{0}}$ to the percent level in a single lens study, it is not sufficient. Even though the time-delay of RXJ1131$-$1231 is the smallest contributor to the error budget in $D_{\Delta t}$ (Suyu et al., 2013), the final precision of $\sim$ 7% is still substantial and mostly dominated by uncertainties in the LOS mass contribution and degeneracies in the lens mass model.

3.3 IFU stellar kinematics

Spectroscopic observations of RXJ1131$-$1231 with state-of-the-art instruments have, so far, only been able to yield a single stellar velocity dispersion measurement (Suyu et al., 2013), due to the faintness of the lens and difficulties in separating the bright quasar light from the galaxy. Yet, future observatories, such as JWST and E-ELT, will be capable of obtaining far more than a single aperture averaged measurement of the stellar kinematics, due to their improved sensitivity and resolution. In order to assess the improvements in constraining $D_{\Delta t}$, when IFU data from JWST become available, we have mocked up stellar kinematics based on the specifications of JWST’s near-infrared spectrograph (NIRSpec). To this end, we have carried out a lensing-only fit to the imaging data and time-delays of RXJ1131$-$1231, with a source resolution of $64\times 64$ pixels. The light model consists of four pseudo-isothermal elliptical profiles (PIEMDs), which are used to mimic a two-component Sérsic distribution (Dutton et al., 2011; Suyu et al., 2014). The mass model consists of a baryonic and non-baryonic component, where the former is obtained by multiplying the light profile with a constant stellar M/L and the latter is accounted for by a NFW halo. The best-fitting model of this fit is then employed to create a luminosity and mass density profile, and complemented by a minimal set of random dynamical parameters, to create a mock kinematic map of the second-order velocity moment according to Eq. 19. Whereas the mock time-delay distance $D_{\Delta t}^{\rm model}$ is based on the best-fitting lensing-only model, the mock input lens distance $D_{\mathrm{d}}^{\rm model}$ is obtained by adopting our standard cosmological model in Sec. 1 with $z_{\mathrm{d}}=0.295$ and $z_{\mathrm{s}}=0.654$.

To simulate observational effects and to assess maps of different quality, we spatially bin the map beforehand via an adaptive spatial binning procedure (Cappellari & Copin, 2003) to a target S/N of 20, 40 and 60 respectively. Especially the latter two are deemed more than sufficient to reliably extract the kinematics across the entire FOV (Falcón-Barroso et al., 2017). The binning is achieved by assuming a S/N in the central (brightest) pixel. Given the parameterised SB distribution, and assuming a Poisson noise dominated regime, the relative pixel intensities can then be translated into a relative 2D S/N map, which is then binned according to the above requirement (Fig. 2). The S/N in the central pixel is obtained by means of JWST’s exposure time calculator (ETC V1.3), where we aimed for a S/N that is both high enough to yield enough spatially resolved measurements and achievable with reasonable on-source integration times. Our final data sets are comprised of three different S/N combinations. A central S/N of 60 (or 100 or 30, respectively), with a target S/N of 40 (or 60 or 20, respectively) in each bin. The combination of the central and target S/N levels are denoted as 60/40 (100/60 & 30/20) hereafter, implying a binning scheme where the central S/N is e.g. chosen to be 60 with a target S/N of 40 across all bins.

The mock kinematics cover a 3″$\times$3″ FOV for JWST at 0.1″/pixel resolution (see Table 1). This mimics a single pointing with JWST, centred on the lens, where a small cycle dither pattern with subarcsecond shifts can be carried out to allow for identification and removal of cosmic rays and detector defects. The FOV contains $\sim$ 900 spaxels. In reality, however, the number of useful spaxels (and hence the final number of bins) will be smaller than the total number within the nominal FOV, due to contamination from B, D and S. These will be masked during the fitting of the stellar spectra and extraction of the stellar kinematics. The loss in spatial information, though, should be minimal given the small point spread function (PSF). Nonetheless, we will also probe a smaller FOV of 2″$\times$2″, to quantify the changes in our cosmological constraints when less data is available. This smaller FOV is a conservative assumption and results in a considerable loss of spatial information when compared to the nominal FOV of 3″$\times$3″, but can be considered as a worst case scenario, where we aim to predict the minimal gain in our cosmological inference. Note that these mock maps are created to harness the full power of JWST’s IFU spectrograph and are in contrast to previous studies of spatially resolved, but only radially averaged profiles of the velocity dispersion (Shajib et al., 2018).

The kinematic data is convolved with a single Gaussian PSF of size 2$\times$ the diffraction limit of JWST and has a FWHM of 0.08″. In reality, more complex shapes, that follow more closely e.g. a Moffat profile, can be expected. But, we are not worried about the actual shape of the real JWST observations, as any shape can also easily be adopted by similarly expanding the PSF with a set of multiple Gaussians, to cover e.g. the extended PSF wings. We add realistic errors to the mock kinematics, where the error in each bin is derived by drawing a random number from a Gaussian distribution with $\mu=0$ and $\sigma_{\rm stat}=(\overline{v_{\mathrm{LOS},l}^{2}})^{1/2}\times\frac{1}{\textrm{(S/N)}_{l}}$ (where ($\overline{v_{\mathrm{LOS},l}^{2}})^{1/2}$ is the mock $(\overline{v_{\mathrm{LOS}}^{2}})^{1/2}$ value at bin position $l$ and $\textrm{(S/N)}_{l}$ its corresponding S/N)777The signal and noise in an IFU spaxel is commonly defined as the median flux and flux standard deviation across the spectral range. Here, we make the assumption that the S/N scales inversely proportional with the errors of our kinematic measurements (Emsellem et al., 2011).. This standard deviation is employed as our true measurement error. Moreover, correlated and uncorrelated uncertainties of 2% each have been added on top of the random Gaussian noise, to account for observational errors that can arise due to e.g. stellar template mismatches. The former simulates a systematic floor in our mock data set, where we utilise the median $(\overline{v_{\mathrm{LOS}}^{2}})^{1/2}\times\ \frac{1}{50}$ across all bins to offset all measurements by a constant value; the latter follows the approach described above by adding again random Gaussian noise with $\mu=0$ and $\sigma_{\rm uncorr}=(\overline{v_{\mathrm{LOS,l}}^{2}})^{1/2}\times\frac{1}{50}$ to all bins across the entire FOV. In summary, we have for each bin

[TABLE]

where $\delta v_{\rm stat}=$ Gaussian[0,$\sigma_{\rm stat}$], $\delta v_{\rm corr}=0.02\,(\overline{v_{\mathrm{LOS}}^{2}})^{1/2}$, and $\delta v_{\rm uncorr}=$ Gaussian[0,$\sigma_{\rm uncorr}$]. The IFU kinematic maps are shown in Fig. 3, where the last column depicts our final mock data, which includes all sources of uncertainty and is employed as our reference data set throughout our joint analysis.

4 Analysis

We construct time-delay strong lensing and stellar dynamical models within the axisymmetric Jeans formalism, as described in Sec. 2, and make use of the high-resolution HST data, time-delays and mock IFU stellar kinematic maps in Sec. 3.

4.1 Setup

To reliably constrain cosmological distances to the percent level, we require flexible and accurate prescriptions of the underlying lens mass distribution. In our joint modelling, we therefore make use of observationally and theoretically motivated mass models. In the case of RXJ1131$-$1231, this includes a SPEMD (which accounts for both the dark and luminous mass) and a COMPOSITE mass model (Suyu et al., 2014), consisting of parameterised fits to the baryonic and non-baryonic matter distribution. Here, the baryonic mass is parameterised via four PIEMDs, which are used to mimic a two-component Sérsic contribution, and a NFW dark halo, and is therefore identical in nature with the mass model that has been used to mock up the IFU kinematics in Sec. 3.3. Both mass models have been shown to provide excellent fits to the strong lensing data (Suyu et al., 2014), but an aperture averaged stellar velocity dispersion measurement was essential to bring both $D_{\Delta t}$ distributions into agreement. While the discrepancy could be resolved by including stellar kinematic data, the precision and accuracy is still limited, and this mass model degeneracy is the main contributor to the error budget, which we aim to constrain further by modelling the mock IFU data set.

For simplicity, we neglect the satellite when mocking up the IFU map (Sec. 3) as well as during the modelling of the strong lensing and stellar kinematic data. The satellite is small enough to result in a loss of only a few spaxels, when being masked during the extraction of the IFU kinematics. More importantly, though, the satellite has a negligible effect for our mass model and cosmological inference (Suyu et al., 2013), as it contributes as little as 1% to the SMD at the lens location. Our analysis relies on two separate mass models, as described above. The modelling parameters for both mass parameterisations are presented in Table 2, along with the dynamical modelling parameters, which have also been used to mock up the kinematic data set. Note that the PIEMDs (i.e. the Sérsic profiles) are fixed during the fitting process. Fits to the SB distribution are carried out beforehand, and the SB distribution is translated into a SMD profile by means of a variable stellar M/L. Moreover, the PA and centroids of the dark and luminous matter distribution in the COMPOSITE mass model are fixed to the same value, to ensure that the projected SMD can be deprojected to an intrinsically axisymmetric mass distribution.

Our respective models include a total of 9 variable parameters with mostly uniform priors for both the SPEMD and COMPOSITE mass models. These priors are weakly informative in a sense that they are bound to ranges around the mock input value, where we have used observationally motivated minimum and maximum values for e.g. the anisotropy $\beta_{z}$ (e.g. Cappellari et al., 2007) and power law slope $\gamma^{\prime}$ (e.g. Auger et al., 2010). The prior ranges are sufficiently large to explore well the PDF, while still allowing for a fast convergence. Note, though, that flat priors are not always suitable for all parameters in the fit, as e.g. the intrinsic shape distribution (and consequently the inclination under the assumption of axisymmetry) is well known to be described by a Gaussian profile (e.g. Weijmans et al., 2014). However, the inclination of RXJ1131$-$1231 is severely limited by the high flattening of the second Sérsic in the fit to the SB distribution and, as a consequence, we have employed a flat prior for the possible range of deprojections. Moreover, the choice of our particular priors has been tested extensively and found to be insignificant for the inference presented in Sec. 4.2. Whereas this indicates that the data is indeed powerful enough to draw credible conclusions from the posterior distribution, irrespective of the prior choice, it cannot be considered as evidence for "truly" non-informative priors in the sense of e.g. Jeffreys priors. We therefore advise to probe the impact of such "naive" non-informative prior assumptions, as the posterior can be highly susceptible to the prior choice in less constraining cases.

While the COMPOSITE model has an additional variable parameter ($r_{\rm s}$), this parameter is not constrained at all, given that it lies well beyond the coverage of the lensing and kinematic data. Using a Gaussian prior instead (Gavazzi et al., 2007), the two mass models (SPEMD and COMPOSITE) have effectively the same number of free parameters. This number is significantly smaller than the total number of variables in the final lens model of RXJ1131$-$1231 in Suyu et al. (2013). However, most of those variables have a negligible effect on the key cosmological parameters, such as $D_{\Delta t}$ and $D_{\mathrm{d}}$, and we therefore adopt those optimised variables as fixed values during the modelling process. The only exception with respect to the optimised values in Suyu et al. (2014) is the PA of the NFW halo and the SPEMD, which are now both aligned with the SB distribution. As axial symmetry is a necessary condition for the construction of the Jeans models, a vastly different PA would violate our underlying modelling assumption. Source resolutions of varying sizes, on the other hand, will be probed, given that the parameter constraints show significant shifts depending on the pixelisation scheme for the AGN host galaxy surface brightness. Besides the uncertainty due to different mass parameterisations, this systematic uncertainty is, in fact, one of the biggest sources of uncertainty for a given mass model, resulting in a distribution that can be $2-3\times$ as wide as for a fixed source grid resolution (Suyu et al., 2013). For each mass model, we will therefore examine eight different source grid resolutions of $54\times 54$, $56\times 56$, $58\times 58$, $60\times 60$, $62\times 62$, $64\times 64$, $66\times 66$ and $68\times 68$ pixels. These source grid resolutions are usually sufficient to achieve an adequate $\chi^{2}$, while stabilising the final modelling constraints towards a common value. The source grid size is chosen such that it contains the entire source intensity distribution, with the outermost source grid pixels converging towards zero intensity values. Our final lensing-only models will then equally weight the constraints from models of different source grid resolutions and different lens mass parameterisations (i.e. COMPOSITE and SPEMD), by combining the individual Markov Chains. An equal weighting is applied in the lensing-only case, since we are incapable of differentiating between different models due to the MSD. In the case of a joint lensing & dynamics run, however, a weighting scheme according to the BIC will be applied, where the dynamical likelihood $P_{D}(d_{D}|\tau_{D})$ will be utilised to perform a model selection.

4.2 Modelling

We visualise the results of our joint strong lensing & stellar dynamical models in Fig. 4 and 5, where we show the marginalised 1D PDFs for our main parameters of interest, i.e. $D_{\Delta t}^{\rm model}$ and $D_{\mathrm{d}}^{\rm model}$  with a S/N of 60/40. The top panels display the constraints from fits to the strong lensing and kinematic data with statistical noise only (hereafter IDEAL, corresponding to the second panel of Fig. 3), whereas the bottom panels show the results from fits to the IFU stellar kinematics including various sources of uncertainty (hereafter FIDUCIAL, corresponding to the fourth panel of Fig. 3). The blue shaded region displays the PDF for models with a COMPOSITE mass distribution, the red shaded region for models with a SPEMD and the grey shaded region the combined PDF from both distributions when the IFU stellar kinematics are included in the fit and a weighting according to the BIC is performed. Note that the blue and red shaded region in Fig. 4 is the PDF from lensing-only models. This is in contrast to Fig. 5, where these also include the stellar kinematics, as the lensing-only models are insensitive to $D_{\mathrm{d}}^{\rm model}$.

Both panels in Fig. 4 clearly show the discrepancy in the time-delay distance, when different mass parameterisations are employed to model the foreground lens, with a double peaked distribution. Since these two models yield a comparable goodness of fit (a manifestation of the MSD), both are equally plausible and an equal weighting is applied in the lensing-only case. The combined COMPOSITE and SPEMD PDF therefore results in a distribution of $1835^{+89}_{-116}$ Mpc (here, the 50th and 16th/84th percentiles of the distribution rather than the mean and standard deviation of a single Gaussian). The discrepancy in the time-delay distance measurement can be resolved by including IFU stellar kinematics, where the combined and weighted distributions converge towards the common mock input value of 1823 Mpc. The 2D kinematics clearly distinguish between both mass models, where the systematic offset in $\overline{v_{\mathrm{LOS}}^{2}}$ is utilised to significantly downweight the contribution from the SPEMD (Fig. 6). A Gaussian fit to the combined and weighted COMPOSITE and SPEMD PDF with stellar kinematics included (i.e. grey shaded region) yields [$\mu_{\rm I},\sigma_{\rm I}$] = [1789 Mpc, 37 Mpc] in the IDEAL case and [$\mu_{\rm F},\sigma_{\rm F}$] = [1794 Mpc, 36 Mpc] for our FIDUCIAL data set respectively. This time-delay distance constraint is a considerable improvement in precision (2.0%) when compared to the lensing-only models (5.6%), even when systematic uncertainties in the stellar kinematics are generously taken into account.

The improvement can be traced back to three effects, in particular, i) a smaller width of the PDF for individual mass parameterisations with different source resolutions, ii) a shift of the mean of the distribution towards the true input time-delay distance and iii) a drastic downweighting of models with a significantly worse goodness of fit. In Fig 7, we illustrate the first two effects by showing the PDFs from lensing-only and joint strong lensing and stellar dynamical models. In both cases, we adopted a COMPOSITE mass distribution and modelled with six different source resolutions. The joint fit to the IFU stellar kinematics considerably reduces the width of the combined PDF from different source resolutions, effectively erasing the low probability wings from lensing-only models in Fig. 4, while shifting the whole distribution towards the input time-delay distance.

A joint fit with stellar kinematics of even higher S/N almost perfectly recovers the input time-delay distance while reaching a precision of 1.7%. In contrast, lower S/N kinematics yield not only worse precision (of 2.5%) but also worse accuracy, where the time-delay distance is only recovered within 1.2$\sigma$. The differences in the constraints for models of different S/N can be attributed to the aforementioned three effects, which are less prominent when the quality of the IFU kinematics degrades.

When it comes to the constraints for the lens distance $D_{\mathrm{d}}^{\rm model}$, we observe a much tighter distribution for our reference S/N of 60/40. Fits to the data with statistical noise only recover remarkably well the input lens distance of 775 Mpc. A Gaussian fit to the combined and weighted PDF from models with a COMPOSITE and SPEMD yields [$\mu_{\rm I},\sigma_{\rm I}$] = [769 Mpc, 14 Mpc]. This is a 1.8% precision measurement for the secondary cosmological distance we aim to infer. Yet, the distribution is clearly biased towards lower distances for fits to our FIDUCIAL data set with [$\mu_{\rm F},\sigma_{\rm F}$] = [734 Mpc, 13 Mpc], which is a consequence of the correlated noise and systematic floor we have added to $(\overline{v_{\mathrm{LOS}}^{2}})^{1/2}$ to mock realistic observational errors. The bias is especially prominent for models with a COMPOSITE mass distribution, where the joint PDF across all source resolutions (with equal weighting) yields a distribution with [$\mu,\sigma$] = [730 Mpc, 13 Mpc]. Keep in mind, however, that the mock data has been created with a source resolution of $64\times 64$ pixels, which is at the edge of the joint COMPOSITE PDF for both $D_{\Delta t}$ and $D_{\mathrm{d}}$ (see e.g. Figure 7). Picking a mock source resolution, and hence an input lens distance, which is closer to the median of the joint PDF in the first place would have partially alleviated this strong bias for our COMPOSITE mass models. Despite our mock source resolution, we observe a similar bias in the joint COMPOSITE and SPEMD PDF across all S/N. Nonetheless, we are able to recover the true lens distance from our joint COMPOSITE and SPEMD PDF within $\sim 3\sigma$, even for our highest S/N pick (with correlated and uncorrelated systematics included), where the model uncertainties are the tightest. We also observe that the $D_{\Delta t}$ constraints are overall preserved, irrespective of the various kinematic noise properties and realisations. This is attributable to the fact that $D_{\Delta t}$ is mainly anchored by the lensing data and time delay measurements. In contrast, information regarding $D_{\mathrm{d}}$ is directly embedded in our mock stellar kinematics. As a consequence, it is highly sensitive to any systematic changes in the data. In fact, $(\overline{v_{\mathrm{LOS}}^{2}})^{1/2}\propto 1/(D_{\mathrm{d}})^{1/2}$. The median $(\overline{v_{\mathrm{LOS}}^{2}})^{1/2}$ difference between the IDEAL and FIDUCIAL data across all bins amounts to $\sim$ 2.1%, which perfectly explains the $\sim$ 4% shift in the mean $D_{\mathrm{d}}$ between both distributions. We note that the inference of $H_{0}$ from the joint constraint of $D_{\Delta t}$ and $D_{\mathrm{d}}$ depends on the kinematic measurements, given the dependence of $D_{\mathrm{d}}$ on kinematics (as illustrated in Section 5).

It is worth mentioning, that correlations between the various parameters in the models will artificially inflate the model uncertainties, since the multidimensional least-squares plane will have more than one unique solution. In practice, this would imply an overestimation of the confidence intervals reported here. We can, however, rule out a high degree of correlation between the various model parameters, as this would manifest itself in unstable solutions with regard to small variations in the observations. The latter is essentially probed by our various realisations of the mock input data, which show that the solutions are generally stable and recover the mock input values regardless of e.g. the noise properties.

The final modelling results are summarised in Table 3, along with our set of complementary models, which we have constructed to assess uncertainties related to modelling i) a smaller FOV, ii) a miscalculation of the PSF size and iii) a single aperture measurement. Especially the latter has been carried out for direct comparison with literature measurements from H0LiCOW, which are based on a single aperture velocity dispersion.

4.3 Sources of uncertainty

To understand the impact of various sources of uncertainty, we complement our FIDUCIAL ("Full FOV") models by fitting to different data sets, by adopting different mass parameterisations and by accounting for PSF mismatches. In all cases, the data have been mocked up and modelled in the same manner as outlined in Sec. 3.3 and 4.1.

4.3.1 Field Of View

In a first test, we fit to a smaller $2\arcsec\times 2\arcsec$ FOV. As mentioned in Sec. 3.3 and 4.1, we have omitted the satellite and AGN images in the data construction and modelling phase. Clearly, both will result in a loss of spatial information, as the affected spaxels will have to be masked when measuring the kinematic moments. In addition to contamination from nearby sources, spatial information will also suffer due to the breakdown of a Poisson noise dominated regime in the remote regions. It is hard to quantify this loss beforehand, given that the final S/N in any given spaxel will be a complex function of the noise properties of the detector, but we try to mimic both effects by drastically reducing the FOV by almost 50%. Keep in mind, however, that this does not translate to a loss of 50% of spatial information. The final number of bins is still 56, compared to 78 for the nominal FOV, and a consequence of the low S/N spaxels beyond $\sim$ 1″ being discarded, which are otherwise massively binned to reach the target S/N.

The modelling results of this run are presented in Table 3, where fits to the IDEAL (i.e. with statistical noise only) and FIDUCIAL (i.e. with correlated and uncorrelated systematics) data are taken into account for a S/N configuration of 60/40. As expected, the constraints for $D_{\Delta t}^{\rm model}$ [$\mu,\sigma$]IDEAL = [1783 Mpc, 42 Mpc] & [$\mu,\sigma$]FIDUCIAL = [1791 Mpc, 42 Mpc] suffer and we obtain a precision of 2.5%, when compared to our reference "Full FOV" models with the same S/N, which has a precision of 2.1%. While we achieve a comparable precision for $D_{\mathrm{d}}^{\rm model}$ with this smaller FOV, it is noteworthy that the bias towards lower distances is slightly more pronounced [$\mu,\sigma$]FIDUCIAL = [732 Mpc, 13 Mpc], where the true lens distance can now only be recovered within $\sim$ 3.2$\sigma$. These findings urge us to aim for the deepest and highest S/N observations, as any loss in quality (i.e. S/N) and/or quantity (i.e. FOV) of the IFU stellar kinematics quickly diminishes any potential gain in the cosmological distance measurements. This is most evident, when we consider the extreme case of a single aperture measurement in Sec. 4.3.3.

4.3.2 Point Spread Function

For the creation of the mock IFU stellar kinematics as well as during the modelling phase, we convolve the predictions of the axisymmetric Jeans equations with a PSF of 0.08″ FWHM size. This PSF is twice as large as the diffraction limit of JWST, but a reasonable choice considering past applications and recent simulations for the next generation of telescopes (Tecza, 2011). In contrast to state-of-the-art AO assisted instruments, this remarkable angular separation is achieved by means of a significantly smaller, flux dominating PSF core. The PSF will, however, be undersampled, given NIRSpec’s pixel size, and a dithering strategy is vital to achieve the nominal spatial resolution. Nonetheless, slight mismatches with the true PSF can be expected when measuring the PSF size from real observations, and we account for this mismatch by convolving the model predictions with a PSF that is roughly 10% larger (i.e. $\sim$ 0.09″).

Our PSF mismatch modelling results are again summarised in Table 3. In light of the sub-pixel PSF size, the cosmological distance constraints are stable across both error assumptions (i.e. for statistical noise only and with correlated and uncorrelated errors on top), yielding almost identical precision on both $D_{\Delta t}^{\rm model}$ and $D_{\mathrm{d}}^{\rm model}$. We therefore omit probing models of different S/N and note that our constraints are insensitive to minor deviations from the true PSF size. Keep also in mind, that the correlation between individual bins is minimal. Our PSF is of the order of JWST NIRSpec’s nominal pixel size of 0.1″. Besides the most central bins, which consist of individual spaxels, most kinematic measurements are essentially independent, given that they are also considerably larger, even if we assume a PSF which is 3$\times$ as large as the diffraction limit.

4.3.3 Single aperture

Stellar kinematics are now commonly employed to break the inherent modelling degeneracies in (time-delay) strong lensing studies (e.g. Treu & Koopmans, 2002a; Koopmans et al., 2003; Treu & Koopmans, 2004; Koopmans et al., 2006). Yet, due to the faintness of the lens and difficulties in separating the bright quasar light from the galaxy, even with state-of-the-art facilities, the data is confined to a single aperture measurement of the stellar velocity dispersion. Moreover, currently employed techniques utilising this kinematic information in strong lensing studies are usually not self-consistent or physically too simple to capture the true complexity of realistic lens galaxies. For instance, most implementations in strong lensing studies assume elliptical lens mass models, but model predictions for the stellar kinematics are based upon a spherically symmetric mass distribution (e.g. Birrer et al., 2019; Jee et al., 2015; Suyu et al., 2013). Similarly, simple assumptions for the velocity anisotropy profile are made (Osipkov, 1979; Merritt, 1985a, b).

In order to assess the impact of using a single aperture measurement on $D_{\Delta t}^{\rm model}$and $D_{\mathrm{d}}^{\rm model}$, under the aforementioned modelling limitations, we mock up stellar kinematics of RXJ1131$-$1231 within a 0.8″$\times$ 0.8″ FOV at JWST resolution, according to the procedure outlined in Sec. 3.3. The kinematics are then luminosity weighted to simulate a single aperture measurement of $\overline{v_{\mathrm{LOS}}^{2}}$. In principle, we would have to split $\overline{v_{\mathrm{LOS}}^{2}}$ into two first order moments (Satoh, 1980), defining the contribution of its ordered vs. random motions, for a straightforward comparison with a stellar velocity dispersion measurement. For simplicity, however, we assume that this massive elliptical lens is dispersion dominated such that $\overline{v_{\mathrm{LOS}}^{2}}\approx\sigma^{2}$ within our FOV.

Contrary to our previous models in this section, we now break the self-consistency of our joint lensing and stellar dynamical models. We run the axisymmetric Jeans models in the spherical limit, by fixing the projected short- vs. long-axis ratio of the luminous and dark matter distribution to $q^{\prime}=0.99$. This is in line with literature studies, where spherical Jeans models have been employed for fitting the stellar velocity dispersion. For both the COMPOSITE and SPEMD, the models can easily recover the single aperture $\overline{v_{\mathrm{LOS}}^{2}}$ measurement of 325$\pm$12 km s*-1*, yielding very similar goodness of fit values across all source resolutions. Ergo, the BIC is not capable of discerning between the two different mass parameterisations, leaving the final distribution still double peaked with D_{\Delta t}^{\rm model}$$=1836^{+89}_{-117} Mpc for our FIDUCIAL models (median and 16th & 84th percentiles). The consequences of anchoring the cosmological distance measurements in RXJ1131$-$1231 on a single aperture velocity dispersion are most noticeable for $D_{\mathrm{d}}^{\rm model}$, where the 1D PDF is only loosely constrained, implying a precision $\geq 10\%$. This is smaller than the 18% found in Jee et al. (2019), based on spherical Jeans models that fit the literature stellar velocity dispersion of 323$\pm$20 km s*-1*, but the difference is likely to be attributed to the smaller errors in our mock data. Even though the precision on the distance measurements degrades substantially with only a single aperture averaged second order velocity moment instead of a 2D kinematic map, the input distances are recovered well within 1$\sigma$, without obvious signs of bias despite the spherical symmetry assumption employed here. This seems to be counterintuitive when compared to our elliptical mass models (see Table 3), which fit the 2D kinematics and can exhibit strong biases (especially in $D_{\mathrm{d}}$), but can be explained by two effects. First, the single aperture measurement is simply a tracer of the enclosed mass, in contrast to the 2D kinematics, which traces the spatially resolved mass distribution in detail. Both the COMPOSITE and SPEMD are capable of recovering the enclosed mass within the errors, with no preference for either model. Given the relatively loose constraints from the single anchor measurement in the single aperture case, the corresponding errors in the cosmological distance inferences from a single lens are sufficiently large to cover any potential bias. With more lenses and thus a reduction in the uncertainty of global parameters (e.g. $H_{0}$), such a bias could then become significant relative to the uncertainty and lifting the spherical assumption for the kinematics (as we have done previously) becomes important. Secondly, RXJ1131$-$1231 is quite roundish with a mean flattening for the light profile and NFW halo of $\sim$0.85 and $\sim$0.75 respectively. Recovering the input kinematics with axisymmetric Jeans models in the spherically symmetric case is therefore feasible without introducing a significant bias. This might change for highly flattened gravitational lenses, which (however) are less likely to be found, given the fact that the lensing cross section increases with galaxy mass and the most massive galaxies being quite roundish in projection.

We also emphasise that our assumption of a dispersion dominated lens within the single aperture FOV is unlikely to introduce a significant bias in our inference. With any net streaming motions present, our mock $V_{\mathrm{rms}}$ data cannot be attributed to $\sigma$ only. As a consequence, the true velocity dispersion must be lower. Assuming, however, net streaming motions of e.g. 50 km s*-1* within the 0.8″$\times$ 0.8″ FOV888We assume low net streaming motions, since i) the aperture FOV only covers half the half-light radius of RXJ1131$-$1231, ii) the contributions of $v$ are luminosity weighted, such that the contribution of the higher velocity wings will be downweighted, and iii) the increase of the lensing cross section increases with galaxy mass and the fact that the most massive galaxies are found to be slow-rotating objects, with low $v/\sigma$ profiles (e.g. Cappellari et al., 2007)., results in a change of only 4 km s*-1* for $\sigma$. Considering our mock $V_{\mathrm{rms}}$ observations of 325$\pm$12 km s*-1*, the change in $\sigma$ is thus well within the error bars. We therefore conclude that our findings for models which employ single aperture kinematics are likely to hold, unless the lens system exhibits significant rotation or is highly flattened. It is worth noting, however, that a single aperture measurement is less powerful in lifting the MSD as a result, due to the minuscule differences also in the kinematic $\chi^{2}$.

4.3.4 Generalised NFW

In a last effort to quantify the systematic uncertainties in our models, we adopt a generalised NFW profile (gNFW, Zhao, 1996), where the halo follows a density distribution according to

[TABLE]

with $\rho_{0}=\delta_{c}\;\rho_{c}$ being a product of the characteristic density $\delta_{c}$ and the critical density $\rho_{c}=3H^{2}/8\pi G$ at the time of halo formation, halo scale radius $r_{s}$ and density slope $\gamma$. The use of a gNFW halo is physically motivated by dissipational cosmological simulations, where the dark halo reacts to an accumulation of the central baryonic component via contraction (Blumenthal et al., 1986). More importantly, though, a mass model with a halo of gNFW form will allow us to better understand the systematics associated with a mass model which is comparably close to the true lens mass distribution. This is particularly interesting, given our general ignorance of the true underlying mass distribution. To this end, we make use of a gNFW halo with a Gaussian prior on the halo slope $\gamma$ with mean and standard deviation $1.0\pm 0.1$ (Mantz et al., 2016), while adopting the same priors for the remaining variables.

When employing the above density parameterisation for the dark halo, we obtain strong biases for the time-delay distance, whereas the lens distance again is susceptible to any systematics in the kinematic measurements. With $D_{\Delta t}^{\rm model}$ [$\mu,\sigma$]IDEAL = [1890 Mpc, 39 Mpc] & [$\mu,\sigma$]FIDUCIAL = [1891 Mpc, 39 Mpc] and $D_{\mathrm{d}}^{\rm model}$ [$\mu,\sigma$]IDEAL = [780 Mpc, 15 Mpc] & [$\mu,\sigma$]FIDUCIAL = [761 Mpc, 14 Mpc], we achieve a 2.1% and 1.9% precision measurement for $D_{\Delta t}^{\rm model}$ and $D_{\mathrm{d}}^{\rm model}$ respectively. However, the accuracy for the time-delay distance measurement suffers significantly, with $D_{\Delta t}^{\rm model}$ only being recovered within $2\sigma$. This simple toy model is another display of the MSD, where comparably good fits to the lensing data are obtained while yielding significantly different distance constraints. Only by means of the BIC due to i) slightly worse kinematic likelihoods of these models and ii) an increase in the model degrees of freedom999Despite our Gaussian prior with mean $\gamma=1.0$, the models strongly converge towards a much shallower density slope of $\gamma\sim 0.6$. Models with a flat prior of $0.5\leq\gamma\leq 1.5$ yielded identical constraints. In contrast to the scale radius $r_{s}$, where we used a Gaussian prior as it is not constrained by data, the halo slope is constrained and therefore constitutes an additional degree of freedom in the gNFW models., can we break this degeneracy. Even if the contribution of the gNFW models to the final distance constraints is non-negligible, in contrast to the SPEMD mass models which have vanishing posterior weights (see Appendix A), the final cosmological inference (see Sec. 5.2) is still vastly improved when compared to the literature single aperture findings.

4.3.5 Other sources of uncertainty

Optimising the likelihood function in Eq. 28, implicitly assumes uncorrelated measurements in the lensing and kinematic data. While this is a sensible assumption for the surface brightness and line-of-sight velocity moments, the time delays between individual pairs of images are correlated to some degree. To test for the impact of covariances between the image pairs, we adopt and include the covariance matrix of the time-delay errors into our models, as measured and reported in (Tewes et al., 2013b). From a single test, based on our COMPOSITE mass model with a source resolution of $64\times 64$ pixels and fitting to our FIDUCIAL data, we find no differences ($<1\%$) for the inferred cosmological distances $D_{\Delta t}$ and $D_{\mathrm{d}}$, which is not surprising given the weak correlations in the measurement errors between individual image pairs (Tewes et al., 2013b, Fig. 8).

When constructing mass models of increased complexity, we have knowingly omitted to probe models with e.g. a radially varying stellar $M/L$. Despite mounting evidence for such an IMF induced change (see e.g. van Dokkum & Conroy, 2010; Martín-Navarro et al., 2015; Conroy et al., 2017), our dark halo already already incorporates radial variations in the total $M/L$, albeit with a different slope than the baryonic component. That is, the increased model degrees of freedom in the gNFW profile allow us to assess the implications for the cosmological distances, if the $M/L$ profile of the model differs from the data and vice versa. Considering the findings in Sec. 4.3.4, even the slightest mismatches in the $M/L$ could have adverse consequences for the cosmological distance measurements, while providing acceptable fits to the lensing-only data. However, assuming that the family of mass models include the true lens model (as is the case in this study), these will be properly downweighted by the BIC. In any case, this convincingly demonstrates the paramount importance of adopting a wide range of plausible mass parameterisations for TDC purposes, as strong offsets from the true time-delay distance will be measured otherwise, even if high-quality IFU kinematics are included in the fit.

5 Cosmological Forecast & Discussion

We describe the cosmological constraints in flat $\Lambda$CDM using the forecasted distance measurements from the last section for both the IDEAL and FIDUCIAL lensing and dynamics models with a S/N of 60/40. Our final inference is based on the "Full FOV", "PSF" and "gNFW" models, given that these make use of the same data sets and allow for a proper evaluation and comparison within the BIC. While the former two models have been probed by both the COMPOSITE and SPEMD mass models, the latter employs a generalised NFW model, which allows us to include a wide range of systematic effects related also to the lens mass parameterisation. We also compare the constraints from using the joint $D_{\Delta t}$-$D_{\mathrm{d}}$ measurement, with that from using only the marginalised $D_{\Delta t}$ measurement.

5.1 Importance sampling with the forecasted distances

In order to obtain a cosmographic forecast, the first step is to get the posterior probability distribution of the cosmological distance measurements, accounting for systematic uncertainties. From the lensing and dynamical modelling detailed in the previous section, we have Markov chains containing the sampled $D_{\Delta t}^{\rm model}$ and $D_{\mathrm{d}}^{\rm model}$ parameters, for various models and set ups. For each data set, we weight the various models (with different lensing mass parametrisation and lensing source grid resolutions) using their BIC values following Eq. 40, where we have estimated $\sigma_{\rm BIC}$ through the scatter in BIC values from models that differ only in lensing source resolutions (given that the source resolutions have a dominant effect on the scatter). We then combine the weighted chains/models, and fit the marginalised $D_{\Delta t}^{\rm model}$-$D_{\mathrm{d}}^{\rm model}$ distribution with a multivariate Gaussian to obtain $P(D_{\Delta t}^{\rm model},D_{\mathrm{d}}^{\rm model}|d_{\mathrm{L}},d_{\mathrm{D}})$.

To further account for the uncertainty due to the external convergence from mass structures along the LOS, we use the $\kappa_{\mathrm{ext}}$ distribution from McCully et al. (2017), which is obtained through a 3-dimensional reconstruction of the mass structures in the field of RXJ1131$-$1231. By reconstructing the mass distribution specific to the RXJ1131$-$1231 sightline, the resulting $\kappa_{\mathrm{ext}}$ distribution is substantially more precise compared to that obtained statistically through galaxy number counts and numerical simulations in Suyu et al. (2014). While the galaxy-number counts approach is thoroughly tested and robust (e.g., Greene et al., 2013; Rusu et al., 2017), there is often a multitudes of sightlines with similar galaxy-number counts with different $\kappa_{\mathrm{ext}}$, giving rise to the scatter/uncertainty in $\kappa_{\mathrm{ext}}$ from this approach. By focusing on the environment of specific lens, the new method developed by McCully et al. (2017) (see also Collett et al., 2013) help reduce the scatter and uncertainty in $\kappa_{\mathrm{ext}}$, and we expect further developments/applications of this new method to yield better constrained $\kappa_{\mathrm{ext}}$ for lens systems in the future. In this paper, we adopt and approximate the $\kappa_{\mathrm{ext}}$ from McCully et al. (2017) as a Gaussian distribution with mean 0.05 and standard deviation of 0.01 (McCully, private comm.), and we also consider the $\kappa_{\mathrm{ext}}$ from the number counts (Suyu et al., 2014) for comparison.

With the posterior probability distribution $P(D_{\Delta t}^{\rm model},D_{\mathrm{d}}^{\rm model}|d_{\mathrm{L}},d_{\mathrm{D}})$, we can relate this to constraints on the cosmological parameters in any background cosmology through importance sampling (e.g., Lewis & Bridle, 2002; Suyu et al., 2010). As an illustration, we consider the constraints on $H_{0}$ specifically for the flat $\Lambda$CDM cosmology, where we adopt uniform priors on $H_{\mathrm{0}}$ between [50,120] km s*-1* Mpc*-1*and on the matter density parameter $\Omega_{\mathrm{m}}$ between $[0.05,0.5]$. We draw $10^{7}$ samples in $\{H_{0},\Omega_{\mathrm{m}}\}$, and compute the corresponding $D_{\Delta t}$ and $D_{\mathrm{d}}$ values given the lens and source redshifts in flat $\Lambda$CDM. For each of these samples, we also draw a value of $\kappa_{\mathrm{ext}}$ from the $\kappa_{\mathrm{ext}}$ distribution, and scale the distances according to Eqs. 13 and 33 to obtain $D_{\Delta t}^{\rm model}$ and $D_{\mathrm{d}}^{\rm model}$. We finally weight the sample by $P(D_{\Delta t}^{\rm model},D_{\mathrm{d}}^{\rm model}|d_{\mathrm{L}},d_{\mathrm{D}}$). From the distribution of the weighted samples, we obtain constraints on $H_{0}$ and $\Omega_{\mathrm{m}}$.

5.2 Forecasted $H_{0}$ constraint in flat $\Lambda$CDM

We show in Fig. 8 the cosmographic constraints for the FIDUCIAL lensing and dynamical models with S/N of 60/40 (magenta). The constraints with statistical errors only (blue) are comparable, though slightly shifted in $D_{\mathrm{d}}$ towards the mock input value of 775 Mpc. After including all sources of uncertainty, we expect to achieve a measurement of $H_{\mathrm{0}}$ $=85.4^{+2.8}_{-2.6}$ km s*-1* Mpc*-1*, with 3.2% precision (defined by the 50th, 16th and 84th percentile), by having high-quality spatially-resolved kinematic data from JWST. The marginalised $D_{\Delta t}$ constraint is $1913^{+59}_{-61}{\rm Mpc}$, which is of similar precision as $H_{0}$. Compared to the 6.6% uncertainty in $D_{\Delta t}$ without spatially resolved kinematic data (Suyu et al., 2014), we are reducing the systematic uncertainty by a factor of $\sim$ 2. Even when accounting for a more conservative $\kappa_{\mathrm{ext}}$ distribution (Suyu et al., 2014), by ray tracing through overdense LOSs in the Millennium Simulations (Springel et al., 2005), a measurement of $H_{\mathrm{0}}$ $=85.1^{+2.9}_{-2.6}$ km s*-1* Mpc*-1* (i.e. with 3.3% precision) is within reach, partly constrained by the assumption of flat $\Lambda$CDM, which restricts the range of plausible $D_{\Delta t}$ and $D_{\mathrm{d}}$ values. That is, most of our improved constraints (a factor of $\sim$2 in precision or 3.3% respectively) stem from our mock 2D JWST kinematics, when compared to the lensing-only case, with an additional improvement of 0.1%, if the smaller scatter in the LOS convergence (McCully et al., 2017) is taken into account.

With spatially resolved kinematics, we would also constrain $D_{\mathrm{d}}$ to $740^{+17}_{-18}$, with 2.4% uncertainty. This measurement is substantially better than the $\sim 18\%$ from the single-aperture average velocity dispersion measurement (Jee et al., 2015, 2019), yet biased towards lower values due to the systematic floor we have added to the FIDUCIAL mock kinematics. Nonetheless, we recover the mock input lens distance within $\sim 2\sigma$. To see whether $D_{\mathrm{d}}$ helps to further constrain $H_{0}$, we repeat the cosmographic forecast above using the marginalised probability distribution of $D_{\Delta t}$, i.e., $P(D_{\Delta t}|d_{\mathrm{L}},d_{\mathrm{D}},d_{\mathrm{env}})=\int{\rm d}{D_{\mathrm{d}}}P(D_{\Delta t},D_{\mathrm{d}}|d_{\mathrm{L}},d_{\mathrm{D}},d_{\mathrm{env}})$. With only the marginalised $D_{\Delta t}$, the constraint on $H_{0}$ degrades slightly to $85.5^{+3.5}_{-3.3}$ (i.e. 4.0% uncertainty) for the FIDUCIAL model. Note, however, that the marginalised $D_{\Delta t}$ measurement is not totally ignoring the distance information in $D_{\mathrm{d}}$, as it has been used to break the MSD in the first place. Moreover, the combination of two distance measurements from a single lens provides not only tight constraints on $H_{\mathrm{0}}$, but adds significant constraining power also for the matter density $\Omega_{\mathrm{m}}$ (Fig. 10).

Given the low redshift of the lens galaxy in RXJ1131$-$1231, $H_{0}$ is primarily constrained by $D_{\Delta t}$ in this case, but the tight constraint on $D_{\mathrm{d}}$ would provide substantial constraints on cosmological models beyond flat $\Lambda$CDM. For illustration, we show the cosmographic constraints for flat $w$CDM cosmology in Fig. 10, where we focus on our FIDUCIAL models with improved (McCully et al., 2017) and conservative (Suyu et al., 2014) $\kappa_{\mathrm{ext}}$ distributions. We note that the prior range on $D_{\Delta t}$ and $D_{\mathrm{d}}$ in the more general $w$CDM model is substantially broader compared to flat $\Lambda$CDM; the conservative $\kappa_{\mathrm{ext}}$ distribution thus leads to a wider $D_{\Delta t}$ distribution, in comparison to the case of $\Lambda$CDM in Fig. 8. With $w=-1.35^{+0.61}_{-0.77}$, the time-independent dark energy is only loosely constrained but, nonetheless, comparable to the combined constraints from 3 single lenses without 2D kinematic data. (Bonvin et al., 2017). Also, such a $D_{\mathrm{d}}$ measurement would serve as a stringent anchor for the inverse distance ladder approach for inferring $H_{0}$ (Jee et al., 2019).

Our forecasted constraints are essentially limited by the S/N of the spatially resolved kinematic maps for breaking lens model degeneracies. As shown in the previous section, higher S/N helps to discriminate between the different mass parameterisations better and hence provide tighter constraints on $D_{\Delta t}$, as the differences in the goodness of fit (and thus in the BIC weighting) become more prominent. For the case of S/N of 60/40 of the FIDUCIAL model, the resulting $H_{0}$ would have an uncertainty of 2.4%. While this S/N can be achieved with reasonable observation times ($\sim$ 6h), higher quality data to constrain $H_{\mathrm{0}}$ further would be difficult to obtain with JWST given the long integration time needed. Future giant segmented mirror telescopes like the E-ELT and TMT, however, could achieve a S/N $\geq$ 100/60 within the same time, owing to their $\sim 5-6\times$ larger aperture.

6 Summary & Outlook

In this paper, we presented a self-consistent joint strong lensing & stellar dynamical modelling machinery for TDC purposes, which employs a pixelated source reconstruction model and the solutions of the Jeans equations in axial symmetry. Our analysis is carried out within the framework of Bayesian statistics and suited, especially, for the study of strong lens configurations for which IFU stellar kinematic data will become available in the near future, by means of the next generation of ground- and space-based telescopes. To assess the performance of the machinery and the expected gain in the inference of cosmological distances and parameters, we mocked up IFU observations of the prominent lens system RXJ1131$-$1231, at JWST NIRSpec resolution. RXJ1131$-$1231 was a particularly natural choice for this study as it is the brightest known lens galaxy for which precise time-delay measurements are already available. The mock stellar kinematic map was based on the best-fitting lensing-only mass model, with a dark and luminous matter contribution, while making random assumptions about the orbital anisotropy and viewing orientation. The mock lens distance has been obtained by assuming a standard cosmological model with $H_{\mathrm{0}}$ $=82.5$ km s*-1* Mpc*-1*, $\Omega_{\mathrm{m}}=0.27$, $\Omega_{\mathrm{\Lambda}}=0.73$ and a lens and source redshift of $z_{\mathrm{d}}=0.295$ and $z_{\mathrm{s}}=0.654$, respectively.

With this suite of data, consisting of deep HST imaging, precise time-delays and mock IFU stellar kinematics of various levels of quality and including various sources of uncertainty, we constructed joint strong lensing & stellar dynamical models. Our models relied on two different mass parameterisations (COMPOSITE and SPEMD), which have been shown to yield significantly different time-delay distances when lensing-only fits are carried out (Suyu et al., 2014). Given the vast amount of information from the spatially resolved kinematics, we utilised the systematic differences in the predicted second-order LOS velocities between the different models, to apply a model selection according to the BIC. The main results of our study can be summarised as follows:

• •

The models recover remarkably well our input time-delay distance $D_{\Delta t}$ ($\leq 1\sigma$), when high-quality IFU stellar kinematics (S/N $\geq$ 60/40) are available. This result is irrespective of the IFU stellar kinematic errors (i.e. assuming purely statistical errors or with correlated and uncorrelated systematics of 2% each included).

• •

The time-delay distance can only be recovered within $1.2\sigma$ or worse, when the S/N of the IFU kinematics degrades below that of our reference 60/40 model.

• •

The lens distance $D_{\mathrm{d}}$ is recovered in all cases within $1\sigma$, when purely statistical errors for the stellar kinematics are assumed. But, a strong offset from the mock lens distance is observed, when the stellar kinematics are systematically biased towards higher or lower values. In these instances, the true lens distance can only be recovered within $3\sigma$ or worse, depending on the systematic offset of the data. Controlling the systematics in the measurement of the stellar kinematics is therefore key for a reliable inference of $D_{\mathrm{d}}$.

• •

The aforementioned results are valid for a 2D map, that covers the LOS velocity distribution within a 3″$\times$3″ FOV (e.g. JWST NIRSpec nominal FOV). Modelling a smaller 2″$\times$2″ FOV, to account for loss of spatial information due to contamination from nearby objects, yields similar accuracy but is less precise ($\Delta D_{\Delta t}^{\rm model}/D_{\Delta t}^{\rm model}=2.5\%$). This highlights the importance of deep and high-quality IFU data, as the gain in the cosmological inference is easily diminished when less spatial information is available.

• •

Small mismatches with the true kinematic PSF size have a negligible impact on the final modelling constraints.

• •

A single aperture stellar velocity dispersion is not very effective in breaking the MSD in RXJ1131$-$1231, yielding a marginal improvement in precision for $D_{\Delta t}^{\rm model}$ over lensing-only models. The constraints for $D_{\mathrm{d}}^{\rm model}$ suffer the most, with a precision $\geq 11\%$.

• •

We achieve a 2.0% precision measurement on $D_{\Delta t}^{\rm model}$, for our FIDUCIAL models with a S/N of 60/40 and including various sources of uncertainty while mocking up the IFU stellar kinematics. Accounting for a wide range of additional sources of systematic uncertainties, as sampled by our test models in Sec. 4.2, this translates to a 3.1% and 3.2% precision measurement on $D_{\Delta t}$ and $H_{\mathrm{0}}$, respectively, in flat $\Lambda$CDM.

• •

A 2.4% precision measurement can be achieved for $D_{\mathrm{d}}$. Yet, this measurement is sensitive to the aforementioned systematics in the stellar kinematics, since the constraints are mainly anchored by the IFU data.

• •

The constraints for $D_{\Delta t}$, $D_{\mathrm{d}}$ and hence $H_{\mathrm{0}}$ improve by a factor of $\sim$ 2, when high-quality IFU stellar kinematics are incorporated in the fit. The improvement can be traced back to three effects in particular, i) a smaller width of the PDF for individual mass models with different source resolutions, ii) a shift of the mean of the distribution towards the true time-delay distance and iii) a drastic downweighting of models with a significantly worse goodness of fit, which is otherwise not feasible due to the MSD.

The increased flexibility of our models allows for a more realistic modelling approach, while circumventing many of the assumptions and limitations of literature time-delay studies. Yet, as the lensing cross section increases with mass, gravitational lenses are likely to be massive elliptical galaxies, which have grown through numerous violent minor and major merger encounters (Wellons et al., 2016). As a consequence, lens galaxies are neither spherical nor elliptical. In fact, recent studies strongly indicate that the most massive galaxies are triaxial (Li et al., 2018), and modelling within an axisymmetric framework might be equally inadequate. It is beyond the scope of this paper to quantify the systematic uncertainties that can be traced back to the violation of axial symmetry, but literature studies show that the reconstructed dynamical masses can be underestimated by as much as 50%, depending on the viewing orientation (Thomas et al., 2007). Taking into account the link between the gravitational potential and the excess time delays, we advise against a modelling within this framework if strong signatures of triaxiality, such as isophotal twists or kinematically decoupled components, are present.

The modelling machinery presented in this paper, along with high-quality IFU data from future space- and ground-based telescopes, provides a promising outlook for constraining cosmological parameters to the few percent level from axisymmetric lenses. Given our forecast for the single lens system RXJ1131$-$1231, an $H_{\mathrm{0}}$ measurement of $\leq$ 2.0% precision could be within reach, if similar gains in precision can be obtained for a total of three lens systems. This would be an important boost in precision when compared to the combination of three such lens systems without 2D kinematic data (Bonvin et al., 2017) and comparable to the current best cosmological probes. It is worth noting, however, that the stellar kinematics have been mocked up by means of a COMPOSITE mass model, which in turn was used to model the suite of data. In reality, though, the set of candidate models is unlikely to contain the true form of the lens potential. As the BIC only applies a relative weighting scheme between all available models, the final accuracy and precision of the cosmological inference heavily relies on an adequate description of the true lens potential. In fact, even the slightest deviations from the true lens potential, as demonstrated by our gNFW models, can result in a biased inference of $H_{\mathrm{0}}$, which is in agreement with Sonnenfeld (2018), where a simple power law model was found to be insufficient to provide an unbiased measurement of $H_{\mathrm{0}}$ in most cases. As a consequence, the study presented here can only be regarded as a best-case scenario, and we strongly encourage to probe a large set of flexible and physically motivated lens mass parameterisations with sufficient degrees of freedom in the radial density profile, to minimise the systematic errors associated with it.

Acknowledgements

We thank G. Chirivì, T. Treu, E. Komatsu and S. Birrer for helpful discussions and comments on the manuscript, and C. McCully for sharing with us the $\kappa_{\rm ext}$ distribution from McCully et al. (2017). AY and SHS thank the Max Planck Society for support through the Max Planck Research Group for SHS. This research was supported in part by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science.

Appendix A Summary of models, BIC values and posterior weights