11institutetext: Max Planck Institute for Astrophysics, Karl-Schwarzschild-
Str. 1, 85748 Garching, Germany 11email<EMAIL_ADDRESS>22institutetext: Technische Universität München, Physik-Department, James-
Franck-Str. 1, 85748 Garching, Germany 33institutetext: Institute of Astronomy
and Astrophysics, Academia Sinica, 11F of ASMAB, No.1, Section 4, Roosevelt
Road, Taipei 10617, Taiwan 44institutetext: Department of Physics and
Astronomy, University of California, Los Angeles, CA 90095, USA
55institutetext: Kavli Institute for the Physics and Mathematics of the
Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Chiba, 277-8583,
Japan

# TDCOSMO VIII: Cosmological distance measurements in light of the mass-sheet
degeneracy - forecasts from strong lensing and IFU stellar kinematics

A. Yıldırım 11 S. H. Suyu  112233 G. C.-F. Chen  E. Komatsu 441155

Time-delay distance measurements of strongly lensed quasars have provided a
powerful and independent probe of the current expansion rate of the Universe
($H_{\mathrm{0}}$). However, in light of the discrepancies between early and
late-time cosmological studies, current efforts revolve around the
characterisation of systematic uncertainties in the methods. In this work, we
focus on the mass-sheet degeneracy (MSD), which is commonly considered as
being a significant source of systematics in time-delay strong lensing
studies, and aim to assess the constraining power provided by integral field
unit (IFU) stellar kinematics. To this end, we approximate the MSD with a
cored, two-parameter extension to the adopted lensing mass profiles (with core
radius $r_{\rm c}$ and mass-sheet parameter $\lambda_{\mathrm{int}}$), which
introduces a full degeneracy between $\lambda_{\mathrm{int}}$ and
$H_{\mathrm{0}}$ from lensing data alone. In addition, we utilise spatially
resolved mock IFU stellar kinematics of time-delay strong lenses, given the
prospects of obtaining such high-quality data with the James Webb Space
Telescope (JWST) in the near future. We construct joint strong lensing and
generalised two-integral axisymmetric Jeans models, where the time delays,
mock imaging and IFU observations are used as input to constrain the mass
profile of lens galaxies at the individual galaxy level, and consequently
yield joint constraints on the time-delay distance ($D_{\Delta t}$) and
angular diameter distance ($D_{\mathrm{d}}$) to the lens. We find that mock
JWST-like stellar kinematics constrain the amount of internal mass sheet that
is physically associated with the lens galaxy and limit its contribution to
the uncertainties of $D_{\Delta t}$ and $D_{\mathrm{d}}$, each at the $\leq$
4% level, without assumptions on the background cosmological model.
Incorporating additional uncertainties due to external mass sheets associated
with mass structures along the lens line-of-sight, these distance constraints
would translate to a $\lesssim$ 4% precision measurement on $H_{\mathrm{0}}$
in flat $\Lambda$CDM cosmology for a single lens. Our study shows that future
IFU stellar kinematics of time-delay lenses will be key in lifting the MSD on
a per lens basis, assuming reasonable and physically motivated core sizes.
However, even in the limit of infinite $r_{\rm c}$, where $D_{\Delta t}$ is
fully degenerate with $\lambda_{\mathrm{int}}$ and is thus not constrained,
stellar kinematics of the deflector, time delays and imaging data will provide
powerful constraints on $D_{\mathrm{d}}$, which becomes the dominant source of
information in the cosmological inference.

###### Key Words.:

distance scale — galaxies: individual — galaxies: kinematics & dynamics —
gravitational lensing: strong — stellar dynamics — methods: data and analysis

##  1 Introduction

The current expansion rate of the Universe ($H_{\mathrm{0}}$) is a fiercely
debated issue. By means of many independent probes, and efforts spanning
multiple decades, a picture has now emerged in which there appears to be a
clear discrepancy between early and late-time cosmological studies (Planck
Collaboration et al., 2020; Riess et al., 2019; Wong et al., 2020; Abbott et
al., 2018; Pesce et al., 2020; Verde et al., 2019). At the heart of this
debate is nothing less than our understanding of our standard cosmological
model and, at face value, this discrepancy can only be resolved by yet unknown
and unaccounted physics (e.g., Poulin et al., 2019; Agrawal et al., 2019; Knox
& Millea, 2020) or by our lack of accounting for systematic uncertainties in
the methods (e.g., Freedman et al., 2019; Shanks et al., 2019).

Given the possibly drastic implications, the assessment of the systematics has
been of particular interest lately. In the context of time-delay strong
lensing (TDSL; Refsdal, 1964; Treu & Marshall, 2016; Suyu et al., 2018), the
prominent mass-sheet degeneracy111A special case of the Source-Position
Transformation (Schneider & Sluse, 2014; Unruh et al., 2017; Wertz et al.,
2018) (MSD; Falco et al., 1985) is commonly quoted as being a fundamental
source of uncertainty (Schneider & Sluse, 2013; Kochanek, 2020). In brief, the
MSD manifests itself through a degeneracy between the lens potential and the
observationally inaccessible source position and size. That is, a
transformation of the former can be compensated by a translation and scaling
of the latter (see Eq. 6-12 in Yıldırım et al., 2020, and references therein).
This mass-sheet transformation (MST) leaves many observables invariant and
thus recovers e.g. the image positions and time delays equally well. However,
with the time-delay distance, and hence $H_{\mathrm{0}}$, being sensitive to
the (functional/radial form of the) lens potential (Sonnenfeld, 2018),
accounting for the MSD is crucial for obtaining an unbiased measurement of
$H_{\mathrm{0}}$ with realistic uncertainties.

To reduce the complexity of the MSD, we follow the classification introduced
in Suyu et al. (2010) and differentiate between two kinds of mass sheets: i)
an internal mass sheet that is physically associated with the strong lens
galaxy and hence affects its stellar kinematics, and ii) an external mass
sheet that is attributed to mass structures along the strong-lens line-of-
sight (LOS) which are not physically associated with the strong lens galaxy
and thus do not affect its stellar kinematics. The H0LiCOW ($H_{\mathrm{0}}$
Lenses in COSMOGRAIL’s222COSmological MOnitoring of GRAvItational Lenses;
https://cosmograil.epfl.ch Wellspring; Suyu et al., 2017) and
TDCOSMO333http://www.tdcosmo.org/ (Millon et al., 2020b) collaborations
account for external mass sheets by explicitly modelling e.g. the nearest and
brightest galaxies (in projection) along the LOS and comparing observed galaxy
counts in the vicinity of the strong-lens system to counts from ray tracing
through simulations. These efforts show that a final precision of 2.4% on
$H_{\mathrm{0}}$ is achieved from an ensemble of 6 time-delay strong lenses,
when well motivated mass models for the lens galaxies are adopted (Wong et
al., 2020). However, if these mass model assumptions are relaxed and a model
with internal mass sheet – that is maximally degenerate with $H_{\mathrm{0}}$
– is used, then the constraint on $H_{\mathrm{0}}$ degrades to 5% or 8%,
respectively, depending on whether or not external priors (obtained from non-
time-delay lenses) for the orbital anisotropy are adopted (Birrer et al.,
2020).

In order to gain back the precision lost by this more general and inclusive
approach, several avenues can be explored. These range from acquiring even
more stellar kinematic data of non-time-delay lenses (to further constrain the
orbital anisotropy and mass profile of lens galaxies at the galaxy population
level) to high-quality, partially resolved stellar kinematics (i.e.,
kinematics within azimuthally averaged rings from the ground) of individual
time-delay lenses (Birrer & Treu, 2021).

Given our previous work on forecasting cosmological distance measurements with
high-spatially resolved stellar kinematics from the next generation
telescopes, such as the James Webb Space Telescope (JWST) or the planned
European Extremely Large Telescope (E-ELT) (Yıldırım et al., 2020, hereafter
Y20), this study aims to reassess the constraining power of joint strong
lensing and stellar dynamical models for cosmological purposes, while
incorporating now a more general family of mass models that closely mimics the
contribution of an internal mass-sheet component. The mass-sheet component
used in this work resembles the parameterisation presented in Blum et al.
(2020, herafter B20) of a cored mass-density distribution, which has been
tailored to be insensitive to the lensing data and is therefore maximally
degenerate with $H_{\mathrm{0}}$. Adopting this cored density distribution
hence allows us to fully explore the uncertainties associated with the MSD and
the prospects of limiting its impact effectively by means of integral field
unit (IFU) data, which is otherwise not possible unless, e.g., aided by
additional distance indicators, the assumption of a cosmological world model
(Chen et al., 2020) or the aforementioned use of prior information (Birrer et
al., 2020).

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_{m}=0.27$,
and a dark energy density of $\Omega_{\Lambda}=0.73$, where our particular
choice for $H_{\mathrm{0}}$ is driven by the time-delay distance measurement
of RXJ1131$-$1231 in Suyu et al. (2014)

##  2 Data

As a test-bed for our study, we will make use of RXJ1131$-$1231 that was
discovered by Sluse et al. (2003). Being part of the H0LiCOW base sample (Suyu
et al., 2017), RXJ1131$-$1231 represents one of the most prominent TDCOSMO
systems with exquisite data, consisting of Hubble Space Telescope (HST)
imaging (Fig. 1) and highly precise time-delay measurements (Tewes et al.,
2013a; Millon et al., 2020a). In addition, an aperture averaged stellar
velocity dispersion measurement for the lensing galaxy of $\sigma=323\pm 20$
km s-1 is available (Suyu et al., 2013).

Following the procedure outlined in Y20, we make use of a hybrid data
strategy, where the archival HST data (GO: 9744; PI: Kochanek) and literature
time-delays (Tewes et al., 2013b) are aided by mock IFU stellar kinematics. In
detail, the HST observations have been used to remock the imaging data. This
is achieved by means of the best-fitting COMPOSITE (i.e. stars & NFW halo)
lensing model achieved on the HST lens image (Fig. 1), as obtained in Suyu et
al. (2014). In this model, the extended source has been reconstructed on a
source grid with a resolution of $64\times 64$ pixels. Unlike the COMPOSITE
lens model in Suyu et al. (2014), we excluded the nearby satellite (S in Fig.
1) in our mock HST image, to ensure consistency with the mock IFU stellar
kinematics, where the satellite has also been omitted (see Sec. 3 for further
explanation). The stellar kinematics of RXJ1131$-$1231 have been generated
based on this COMPOSITE lens mass model and using arbitrary dynamical model
parameters (see Sec. 3). Utilising the axisymmetric Jeans formalism for
mocking up the stellar kinematics (Binney & Tremaine, 1987; Cappellari, 2008;
Barnabè et al., 2012), we obtain IFU data at JWST resolution (i.e.
0.1″/pixel), covering a 3″$\times$3″ field-of-view (FOV). The mock IFU data
closely resemble the Near-Infrared Spectrograph (NIRSPEC) properties and
characteristics, assuming that the kinematics have been measured with a single
pointing centred on RXJ1131$-$1231 (Fig. 1).

The mock kinematics are further convolved with a single Gaussian PSF of 0.08″
FWHM, which is already two times larger than the diffraction limit of JWST.
Difficulties in properly measuring the PSF size have been shown to be of
minimal concern (Y20), which is why this PSF is employed for both the mock and
modelling stage.

Since the signal-to-noise (S/N) in each spatial resolution element is usually
insufficient to reliably measure the stellar kinematics across the entire FOV,
we perform a spatial binning of the mock data (Cappellari & Copin, 2003). To
this end, the surface brightness (SB) distribution at JWST resolution is
obtained from parameterised fits to the HST imaging data. This SB distribution
of RXJ1131$-$1231 is then used to achieve a desired S/N level, by means of
JWST’s exposure time calculator (ETC V1.3) Within the 3″$\times$3″ FOV, this
S/N map is, in turn, used to spatially bin the individual pixels to a target
S/N of $\sim$ 40 per bin. Our procedure results in $\sim$ 80 spatially
resolved measurements of the line-of-sight velocity distribution (LOSVD)
across the entire FOV, achievable with reasonable on-source integration times
of 8-10h.

To imitate realistic observations, various sources of statistic and systematic
uncertainties are added to the mock input kinematics
$\bigg{(}\sqrt{\overline{v_{\mathrm{LOS,\,input}}^{2}}}\bigg{)}$. Here, we
focus on two individual data sets. Following the notation in Y20, the IDEAL
data incorporate random Gaussian errors only ($\delta v_{\rm{stat}}$). More
specifically, the error in each bin $l$ is drawn from a Gaussian distribution
with a mean of zero and a standard deviation
($\sigma_{\rm{stat,\,\textit{l}}}$) that is inversely proportional to its
signal-to-noise (S/N). Given that all bins reach the target S/N of 40, this
approach yields $\sigma_{\rm{stat,\,\textit{l}}}\approx 0.025\times\
\sqrt{\overline{v_{\mathrm{LOS,\,input,\,\textit{l}}}^{2}}}$. The IDEAL data
hence reflect a best-case scenario for breaking the MSD, where the results
will be unaffected by e.g. systematic effects in the measurement of the
stellar kinematics. In the second data set, which we label as FIDUCIAL, we
incorporate systematic errors in addition to the random Gaussian noise in the
IDEAL data set. That is, the FIDUCIAL data set contains a 2% systematic shift
downward of all mock
$\sqrt{\overline{v_{\mathrm{LOS,\,input,\,\textit{l}}}^{2}}}$ values
$\bigg{(}\delta
v_{\rm{sys}}=-0.02\times\sqrt{\overline{v_{\mathrm{LOS,\,med}}^{2}}}\bigg{)}$,
where $\bigg{(}\sqrt{\overline{v_{\mathrm{LOS,\,med}}^{2}}}\bigg{)}$ is the
median of $\sqrt{\overline{v_{\mathrm{LOS,\,input,\,\textit{l}}}^{2}}}$ from
all bins. This systematic shift to the downside is introduced to account for
possible systematics in the measurement of the stellar kinematics (e.g. due to
stellar template mismatches). The mock kinematics can be quickly summarised as
follows:

$\sqrt{\overline{v_{\mathrm{LOS,\,data,\,\textit{l}}}^{2}}}=\underbrace{\overbrace{\sqrt{\overline{v_{\mathrm{LOS,\,input,\,\textit{l}}}^{2}}}+\delta
v_{\rm{stat}}}^{\rm{{\color[rgb]{1,.5,0}IDEAL}}}+\delta
v_{\rm{sys}}}_{\rm{{\color[rgb]{0,0,1}FIDUCIAL}}}$ (1)

To check for consistency with real observations, we have collapsed the mock
IFU data within an 0.8″ wide aperture for the IDEAL case, which yields a mock
single second-order velocity moment along the LOS of
$\sqrt{\overline{v_{\mathrm{LOS,ap}}^{2}}}=325\pm 12$ km s-1. This mock value
is in excellent agreement with the literature stellar velocity dispersion
measurement, assuming that rotation is minimal, and confirms the validity of
our approach.

The mock HST observations, time-delay measurements and mock kinematics are
then employed as constraints for our joint modelling machinery and the
modelling results for both mock data sets (i.e. IDEAL and FIDUCIAL) are
presented and discussed in Sec. 3 and 4.

Figure 1: HST ACS F814W imaging cutout of RXJ1131$-$1231, illustrating the
prominent lens configuration with a quadruply imaged background quasar (A, B,
C & D) and a nearby satellite (S). Spectroscopic measurements locate the lens
and quasar at redshifts $z_{\mathrm{d}}=0.295$ and $z_{\mathrm{s}}=0.654$,
respectively. We generated mock imaging data, closely resembling the archival
HST data shown here, based on a best-fitting lensing model and excluding the
satellite. Overlaid is the JWST NIRSpec nominal FOV of 3″$\times$3″, within
which we create mock stellar kinematics of the foreground lens at 0.1″/pixel
resolution. The FOV is oriented such that the x-axis is aligned with the
galaxy major axis.

##  3 Theoretical Framework

###  3.1 Strong lensing under the internal MST

We adopt the COMPOSITE mass model of RXJ1131$-$1231 as our baseline model
(Y20). This model consists of a baryonic and non-baryonic component. The
former is parameterised by multiple pseudo-isothermal elliptical profiles
(PIEMDs; Kassiola & Kovner, 1993), to resemble a two-component Sérsic
distribution, multiplied with a spatially constant stellar mass-to-light ratio
($\Upsilon$), and the latter by a NFW profile (Navarro et al., 1996, 1997).
Moreover, we add a cored density distribution to the COMPOSITE model, here
also parameterised as a PIEMD, with a core radius of $\theta_{\rm
c}=20\arcsec$, which is well beyond the Einstein radius ($\theta_{\rm
E}\approx 1.6\arcsec$)444We refer to $\theta_{\rm c}$ as the angular core
radius. The physical core radius $r_{\rm c}$ in units of kpc is obtained via
$r_{\rm c}=\theta_{\rm c}\,D_{\mathrm{d}}$. See also Sec. 22. As introduced in
B20, such a cored density distribution is built to mimic an internal mass
sheet with constant convergence $\kappa_{\mathrm{int}}$. In concordance with
earlier notations, we define

$\lambda_{\mathrm{int}}\equiv 1-\kappa_{\mathrm{int}},$ (2)

with

$\lambda_{\mathrm{int}}\in[0.5,1.5].$ (3)

By scaling the total convergence profile according to

$\kappa_{\rm
mst,int}(\theta)=(1-\lambda_{\mathrm{int}})+\lambda_{\mathrm{int}}\;\kappa_{\rm
composite}(\theta),$ (4)

where $\kappa_{\rm composite}(\theta)$ is the original COMPOSITE profile, this
MST leaves all observables invariant, provided the time-delay distance
($D_{\Delta t}$)

$D_{\Delta
t}\equiv(1+z_{\mathrm{d}})\;\frac{D_{\mathrm{d}}D_{\mathrm{s}}}{D_{\mathrm{ds}}}\propto
H_{0}^{-1},$ (5)

scales as

$D_{{\Delta t},{\rm int}}=\frac{{D_{\Delta t,{\rm
composite}}}}{\lambda_{\mathrm{int}}}.$ (6)

Here $z_{\mathrm{d}}$ is the redshift to the lens, $D_{\mathrm{d}}$,
$D_{\mathrm{s}}$ and $D_{\mathrm{ds}}$ are angular diameter distances to the
lens, to the source and between the lens and the source, respectively.
Incorporating further the effect of an external mass-sheet
($\kappa_{\mathrm{ext}}$) yields

$D_{\Delta t}=\frac{D_{{\Delta t},{\rm int}}}{1-\kappa_{\mathrm{ext}}},$ (7)

whereas $D_{\mathrm{d}}$ is fully invariant under the external MST555The
angular diameter distance to the lens $D_{\mathrm{d}}$ is also approximately
invariant under the internal MST with a single aperture averaged velocity
dispersion and for $\lambda_{\mathrm{int}}$ $\in[0.8,1.2]$(Chen et al., 2020).
(Jee et al., 2015)

$D_{\mathrm{d}}=D_{\mathrm{d},{\rm int}}.$ (8)

In concordance with B20, our adopted cored density distribution has the
special attribute of leaving the mean convergence within $\theta_{\rm E}$
($\sim$ 1.6″) virtually identical and making it thus inaccessible for strong
lensing studies alone. However, the above transformation introduces changes in
the radial 3D mass distribution, depending on the extent of the core radius,
to which spatially resolved stellar kinematics can be sensitive. Here, the
choice of a core radius is particularly important. Our core radius of 20″
might seem deliberate at first, but was adopted after intensive testing.
Larger core radii will mimic a scenario in which the original density profile
of the galaxy ($\kappa_{\rm composite}$) is effectively unchanged (see
Appendix C), i.e. both the lensing and kinematic data will be insensitive,
whereas core radii which are too small will have a more noticeable imprint on
the stellar kinematics, but might not be a valid mass-sheet transformation in
the sense that the lensing data will start to become sensitive to such a
transformation of the mass profile (see also Fig.3 in, Birrer et al., 2020).
The core radius of 20″ was chosen to probe a worst-case scenario of sorts,
where the MST is still intact from a lensing point of view while the imprint
on the kinematic data is big enough to allow for meaningful constraints on
$\lambda_{\mathrm{int}}$. Our distinction between an internal and external
mass-sheet is therefore a practical one; we assume that the internal MST does
not conspire in a way that it becomes undetectable for both lensing and
stellar kinematics within $\theta_{\rm E}$. We also note that the adopted core
radius of 20″ corresponds to a physical radius of $\approx 80$ kpc at our lens
distance. It is therefore sufficiently large to test physically motivated
mechanisms that could produce extended cored density distributions, such as
those put forward by ultralight dark matter (see e.g., Blum & Teodori, 2021).

###  3.2 Stellar dynamics under the internal MST

We briefly revisit the theoretical framework for mocking up and modelling the
stellar kinematics. As outlined in Y20, our models are based on the solution
of the axisymmetric Jeans equations (Jeans, 1922; Binney & Tremaine, 1987),
assuming alignment of the velocity ellipsoid with the cylindrical coordinate
($R,\phi,z$) system as well as a flattening in the meridional plane
(Cappellari et al., 2007; Cappellari, 2008). This yields the two equations,

$\frac{\beta_{z}\rho_{*}\overline{v_{z}^{2}}-\rho_{*}\overline{v_{\phi}^{2}}}{R}+\frac{\partial(\beta_{z}\rho_{*}\overline{v_{z}^{2}})}{\partial
R}=-\rho_{*}\frac{\partial\psi_{\mathrm{D}}}{\partial R}$ (9)
$\frac{\partial(\rho_{*}\overline{v_{z}^{2}})}{\partial
z}=-\rho_{*}\frac{\partial\psi_{\mathrm{D}}}{\partial z}$ (10)

which link the tracer density ($\rho_{*}$) and gravitational potential
($\psi_{D}$) to three intrinsic second-order velocity moments
($\overline{v_{z}^{2}},\overline{v_{\phi}^{2}},\overline{v_{R}^{2}}$). Once
the tracer density is measured (by means of the SB distribution) and a
gravitational potential is adopted, we can readily solve Eq. 10 for
$\overline{v_{z}^{2}}$ and Eq. 9 for $\overline{v_{\phi}^{2}}$ and obtain
$\overline{v_{R}^{2}}$ via
$\beta_{z}=1-\overline{v_{z}^{2}}/\overline{v_{R}^{2}}$. The intrinsic second-
order moments can then be projected along the LOS to predict a second-order
velocity moment

$\begin{split}\overline{v_{\mathrm{LOS}}^{2}}=&\frac{1}{\mu(x^{\prime},y^{\prime})}\;\int_{-\infty}^{\infty}\;\rho_{*}\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}$
(11)

where $x^{\prime}$ and $y^{\prime}$ are the cartesian coordinates on the plane
of the sky, $z^{\prime}$ is the coordinate along the LOS, $i$ is the
inclination angle (with 90° being edge-on), $\mu$ is the observed SB and $v$
and $\sigma$ the observed mean LOS velocity and velocity dispersion.

In our self-consistent modelling approach, both the intrinsic tracer and mass
density for the dynamics are constructed by deprojecting the SB and surface
mass density (SMD) distribution from the lens models, which are constrained by
the quasar image positions, time delays and the spatially extended Einstein
ring. As mentioned in Sec. 3.1, the projected SMD consists of a COMPOSITE and
PIEMD profile and is subject to the internal MST. More specifically,

$\Sigma_{\rm smd}(\theta)=\Sigma_{\rm crit}\times\ \kappa_{\rm
mst,int}(\theta),$ (12)

where $\Sigma_{\rm smd}$ denotes the (projected) physical surface mass density
and

$\Sigma_{\mathrm{crit}}=\frac{c^{2}}{4\pi G}\;\frac{D_{{\Delta t},{\rm
int}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm
int}}}\;\frac{1}{D_{\mathrm{d},{\rm int}}}$ (13)

is the critical surface mass density666We can express $\Sigma_{\mathrm{crit}}$
here in terms of $D_{{\Delta t},{\rm int}}$ and $D_{\mathrm{d},{\rm int}}$,
instead of $D_{\Delta t}$ and $D_{\mathrm{d}}$, because a subsequent external
MST (after the internal MST) does not affect $\Sigma_{\rm smd}$ as the terms
connected to $\kappa_{\mathrm{ext}}$ cancel out (Y20).. With Eq. 4, we rewrite
the projected SMD as

$\displaystyle\Sigma_{\rm smd}(\theta)$
$\displaystyle=\bigg{[}\frac{c^{2}}{4\pi G}\;\frac{{D_{{\Delta t},{\rm
int}}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm
int}}}\;\frac{1}{D_{\mathrm{d},{\rm int}}}\bigg{]}$ (14)
$\displaystyle\times\bigg{[}(1-\lambda_{\mathrm{int}})+\lambda_{\mathrm{int}}\;\kappa_{\rm
composite}(\theta)\bigg{]}.$

The SMD is expanded by multiple Gaussians (MGE; Emsellem et al., 1994) and
deprojected into a 3D mass density. Assuming an oblate axisymmetric system,
the intrinsic total density of each Gaussian in the fit (MGEfit; Cappellari,
2002) is obtained via

$\rho_{i}(R,z)=\frac{M_{i}}{(\sigma_{i}\sqrt{2\pi})^{3}\;q_{i}}\;\exp\bigg{[}-\frac{1}{2\sigma_{i}^{2}}\;\Bigg{(}R^{2}+\frac{z^{2}}{q_{i}^{2}}\Bigg{)}\bigg{]},$
(15)

and related to $\psi_{\rm D}$ through Poisson’s equation. We use index $i$ to
denote the MGE components of the mass density. It is of importance to note
that, based on our finite core size, our treatment and implementation of the
MST is an approximate one, with $\kappa_{\mathrm{int}}\approx\mathrm{const.}$
within $\theta_{\rm c}$. This, however, is a good enough approximation of the
perfect MST (i.e. $\theta_{\rm c}\to\infty$), where the deviation in e.g. the
deflection angles are small enough (i.e. $\delta_{\theta}(\theta_{\rm
c}\to\infty)\approx\ \delta_{\theta}(\theta_{\rm c}=20\arcsec)$), such that
the models fully capture the degeneracies associated with it (see Sec 4).

In Eq. 15,
$M_{i}=2\,\pi\,\sigma_{i}^{{}^{\prime}2}\,q^{\prime}_{i}\,\Sigma_{\mathrm{M},i}$
denotes the total mass of the individual Gaussians, $\sigma^{\prime}_{i}$ the
projected dispersion, $q^{\prime}_{i}$ the projected flattening and
$\Sigma_{\mathrm{M},i}$ the peak surface mass density. In our case of a MGE of
the SMD, the gravitational potential assigned to each Gaussian is

$\psi_{{\rm
D},i}(R,z)=-\frac{2G\,M_{i}}{\sqrt{2\pi}\sigma_{i}}\;\int_{0}^{1}\;\mathcal{F}_{i}(u)du$
(16)

with

$\mathcal{F}_{i}(u)=\exp\bigg{[}-\frac{u^{2}}{2\sigma^{2}_{i}}\Bigg{(}R^{2}+\frac{z^{2}}{\mathcal{Q}^{2}_{i}(u)}\Bigg{)}\bigg{]}\frac{1}{\mathcal{Q}_{i}(u)}$
(17)

and $\mathcal{Q}^{2}_{i}(u)=1-(1-q^{2}_{i})\;u^{2}$. Here, $\sigma_{i}$ is the
intrinsic dispersion of the Gaussians, $q_{i}$ their intrinsic long vs. short
axis ratio, and ($R$, $z$) the intrinsic cylindrical coordinate system.
Similarly, the SB distribution can also be expanded by multiple Gaussians.
Denoting the MGE components of the luminosity distribution by index $j$ (in
contrast to $i$ for the mass density), the intrinsic luminosity density of
each Gaussian is

$\rho_{*,j}(R,z)=\frac{L_{j}}{(\sigma_{j}\sqrt{2\pi})^{3}\;q_{j}}\;\exp\bigg{[}-\frac{1}{2\sigma_{j}^{2}}\;\Bigg{(}R^{2}+\frac{z^{2}}{q_{j}^{2}}\Bigg{)}\bigg{]},$
(18)

with $L_{j}$ now expressing the total luminosity of the individual Gaussians.
In general, the set of Gaussians describing the SB and SMD are not identical
(i.e. $i\neq j$ and hence $\sigma_{i}\neq\sigma_{j}$ and $q_{i}\neq q_{j}$)
unless mass-follows-light. The second-order moments of a single Gaussian
tracer distribution – embedded in the gravitational potential $\psi_{\rm
D}=\sum_{i}\psi_{{\rm D},i}$ – is then obtained via

$[\overline{v_{z}^{2}}]_{j}=4\pi\,G\sum_{i}\frac{\Sigma_{\mathrm{M},i}\,q^{\prime}_{i}}{\sqrt{2\pi}\sigma_{i}}\int_{0}^{1}q_{j}^{2}\sigma_{j}^{2}\frac{\mathcal{F}_{i}(u)u^{2}}{1-\mathcal{F}_{ji}u^{2}}du,$
(19)

with $\mathcal{F}_{ji}=1-q_{i}^{2}-q_{j}^{2}\sigma_{j}^{2}/\sigma_{i}$ (see
Cappellari, 2008; van de Ven et al., 2010, for analytical expressions of all
velocity moments involved in Eq. 9 and 10). Note that absolute lengths in the
plane of the sky scale linearly with the distance ($r\propto
D_{\mathrm{d},{\rm int}}$). Hence, when moving from angular to physical units,
this dependency translates into

$\overline{v_{\mathrm{LOS}}^{2}}\propto 1/D_{\mathrm{d},{\rm int}}.$ (20)

It is immediately evident from Eq. 14, that the deprojected (i.e. intrinsic)
3D mass distribution is subject to the mass-sheet transformation. The impact
on the radial 3D mass distribution is highlighted in Fig. 2, where we scale
$\lambda_{\mathrm{int}}$ (and hence $D_{{\Delta t},{\rm int}}$ and
$\kappa_{\rm mst,int}$) with respect to the mock input values. Differences in
the total enclosed mass within JWST’s FOV amount to roughly 5% for
$\lambda_{\mathrm{int}}=0.9$. Moreover, the mass-sheet transformation
translates not only to differences in the total enclosed mass, but also
changes the slope of the mass profile to which spatially resolved kinematics
are sensitive (Cappellari, 2020).

Figure 2: Top: total enclosed mass profile in 3D as a function of intrinsic
radius ($\sqrt{R^{2}+z^{2}}$), after deprojecting the SMD into an oblate
axisymmetric system. The black line shows the mock input mass profile (i.e.
the COMPOSITE model without a mass sheet) and the red line shows the mass-
sheet transformed profile with $\lambda_{\mathrm{int}}=0.9$. The profile
without a mass-sheet transformation, but with $D_{\mathrm{d},{\rm int}}$
scaled by +5%, is shown in blue. The vertical dashed lines indicate the extent
of literature, aperture averaged measurement of the velocity dispersion
($\sigma_{\rm ap}$) and the Einstein radius ($\theta_{\rm E}$) in projection.
The grey shaded area shows the region where the lensed imaging data constrain
the mass model and the striped area roughly corresponds to JWST’s FOV. Bottom:
the percent differences in the total enclosed mass as a function of radius
when compared to the mock input mass profile (black). The figure illustrates
that a single aperture averaged velocity dispersion is insufficient for
constraining the mass sheet, given the minuscule differences in the total
enclosed mass within that radius and errors on the kinematics of the order of
5%. On the other hand, spatially resolved kinematics are sensitive to both,
differences in the enclosed mass as well as the mass slope. As is evident from
the figure, $D_{\mathrm{d},{\rm int}}$ and $\lambda_{\mathrm{int}}$ are not
degenerate. That is, changes in $\lambda_{\mathrm{int}}$ (and hence
$D_{{\Delta t},{\rm int}}$) cannot be compensated by a corresponding scaling
of the lens distance $D_{\mathrm{d},{\rm int}}$. Numerical noise from the MGE
is present, but its contribution to the observed differences is subdominant
with respect to the final modelling uncertainties.

###  3.3 Models

We construct joint strong lensing and stellar dynamical models as introduced
in Y20. In brief, these models rely on a pixelated source fitting algorithm
(Suyu et al., 2006, 2010, 2013) and the solution of the Jeans equations in
axisymmetry (Cappellari, 2008, JAMPY). The model parameter space thus includes
those listed in Table 1 for the COMPOSITE mass model and SPEMD, including the
time-delay distance $D_{{\Delta t},{\rm int}}$, the lens distance
$D_{\mathrm{d},{\rm int}}$, the orbital anisotropy parameter $\beta_{z}$,
inclination $i$ and absolute strength/amplitude of the internal mass-sheet
parameter $\lambda_{\mathrm{int}}$. In line with Y20, we do not include the
satellite in our fits to both the lensing and kinematic data. Given that the
contribution of the satellite to the SMD is $\leq$ 1% at the lens centre, any
uncertainties due to the missing satellite are subdominant with respect to the
other sources of uncertainties in the models. Moreover, we do not constrain
the parameter range by physical arguments other than the fact that the mass-
sheet transformation is not supposed to yield negative convergences within a
FOV that is more than twice as large as the core radius, which penalises
models with $\lambda_{\mathrm{int}}\gtrapprox 1.05$777Our FOV within which we
carry out the MGE has to be large enough to avoid a premature exponential fall
off of the Gaussians, which can adversely impact the predicted
$\overline{v_{\mathrm{LOS}}^{2}}$. We find that a FOV $\geq 2\times\theta_{\rm
c}$ is sufficiently large to adequately describe the profile of the cored
density distribution. Within this FOV, however, the convergence profile is not
allowed to be negative, since this would yield negative SMDs. This translates
into an upper bound for our constraints on $\lambda_{\mathrm{int}}$ of
$\approx 1.05$. (see Sec. 4). Besides that, the priors are identical to those
adopted in Y20, with $\lambda_{\mathrm{int}}$ being flat within the
aforementioned prior range of $[0.5,1.5]$. Note that our kinematic data has
been mocked up without any contribution from a mass sheet (i.e.
$\lambda_{\mathrm{int}}=1.0$, as described in Sec. 2). Moreover, we do not
assume a ratio of $D_{\rm d}/D_{\rm ds}$. That is, our inferred values for
$D_{{\Delta t},{\rm int}}$ and $D_{\mathrm{d},{\rm int}}$ are based solely on
the data and are not constrained by assuming a cosmological world model.

Table 1: Model parameters and priors for our joint strong lensing & dynamical
models, including the cosmological distances, the SPEMD, the COMPOSITE mass
distribution, the PIEMD (parameterising the mass sheet) and the dynamical
variables, consisting of models with a constant ($\beta_{z}(r)=const.$) and
spatially varying ($\rm{log_{10}}\,(\beta_{z}(r))=\alpha+\beta\times\
\rm{log_{10}}\,(r)$) anisotropy profile. The mock IFU data set is based on the
best-fitting COMPOSITE lensing-only model with a source resolution of
$64\times 64$ pixels and arbitrary values for the dynamical parameters. The
mock cosmological distances are based on the best-fitting lensing-only model
for $D_{{\Delta t},{\rm int}}$ and assuming $H_{\mathrm{0}}$ $=82.5$ km s-1
Mpc-1, $\Omega_{\mathrm{m}}=0.27$, $\Omega_{\mathrm{\Lambda}}=0.73$,
$z_{\mathrm{d}}=0.295$ and $z_{\mathrm{s}}=0.654$ for determining
$D_{\mathrm{d},{\rm int}}$.

Description Parameters Mock input values Prior type Prior range Distances
Model time-delay distance [Mpc] $D_{{\Delta t},{\rm int}}$ 1823.42 Flat [1000,
3000] Model lens distance [Mpc] $D_{\mathrm{d},{\rm int}}$ 775.00 Flat [600,
1000] SPEMD Flattening $q$ Flat [0.2, 1.0] Einstein radius [arcsec]
$\theta_{\mathrm{E}}$ Flat [0.01, 2.0] Power law slope $\gamma_{\rm pl}$ Flat
[0.2, 0.8] External shear strength $\gamma_{\mathrm{{ext}}}$ Flat [0.0, 0.2]
External shear position angle [°] $\phi_{\mathrm{ext}}$ Flat [0.0, 360.0]
COMPOSITE Stellar M/L $\Upsilon_{\star}$ 2.09 Flat [0.5, 2.5] NFW flattening
$q$ 0.73 Flat [0.2, 1.0] NFW Einstein radius [arcsec] $\theta_{\mathrm{E}}$
0.20 Flat [0.01, 2.0] NFW scale radius [arcsec] $r_{\mathrm{s}}$ 22.53
Gaussian [18.6, 2.6] External shear strength $\gamma_{\,\mathrm{ext}}$ 0.08
Flat [0.0, 0.2] External shear position angle [°] $\phi_{\mathrm{ext}}$ 1.42
Flat [0.0, 360.0] PIEMD Einstein radius [arcsec] $\lambda_{\mathrm{int}}$ 1.00
Flat [0.5,1.5] Dynamics Anisotropy $\beta_{z}(r)=const.$ 0.15 Flat [$-0.3$,
0.3] Inclination [°] $i$ 84.26 Flat [80.0, 90.0] Dynamics Anisotropy
$\rm{log_{10}}(\beta_{z}(r))=\alpha+\beta\times\ \rm{log_{10}}(r)$ Pivot
$\alpha$ Flat [$-0.3$, 0.3] Slope $\beta$ Flat [$-0.5$, 0.5] Inclination [°]
$i$ Flat [80.0, 90.0]

Note: the stellar M/L is the constant multiplied to the observed light
intensity distribution in the mock imaging data to obtain the dimensionless
surface mass density $\kappa$. For a real system with real observations, this
value can then be converted to physical units such as M⊙ and L⊙.

We explore the posterior probability density function (PDF) by means of emcee
(Foreman-Mackey et al., 2013). We investigate different modelling runs by
adopting different source grid resolutions, which (for a given mass model
parameterisation and in the absence of a MSD) constitutes the biggest source
of uncertainty on $H_{\mathrm{0}}$ (Suyu et al., 2013, Y20). The grid
resolutions span a range of 56 $\times$ 56 to 68 $\times$ 68 pixels in the
source plane, with source pixel sizes roughly of the order of 0.05$\pm$0.01
″/pixel. The grid size is sufficiently large to allow for a proper
reconstruction of the extended source, without artificially constraining the
parameter space in $\lambda_{\mathrm{int}}$, as will be shown for the lensing-
only models in Sec. 4.1. In a last step, we then combine the constraints from
different source resolutions by assigning equal weights to the chains.

In the case of the FIDUCIAL data set, we sample over an additional velocity
fractional shift parameter $f_{\rm v}$ in the model to account for the
systematic uncertainties. For a given value of $f_{\rm v}$, the mock kinematic
data $\sqrt{\overline{v_{\mathrm{LOS,\,data,\,\textit{l}}}^{2}}}$ in Eq.1 is
scaled by $(1-f_{\rm
v})\sqrt{\overline{v_{\mathrm{LOS,\,data,\,\textit{l}}}^{2}}}$. We consider a
uniform distribution of $f_{\rm v}\in[-0.02,0.02]$, i.e., we assume the mock
kinematic data in all kinematic bins have equal likelihood in being
systematically shifted by an amount within 2%. Such systematic fractional
offset in the measured kinematic data could arise from, e.g., stellar template
mismatches. In practice, such as systematic shift in the kinematic data leads
to a shift in $D_{\mathrm{d},{\rm int}}$, leaving other model parameters
unchanged. To obtain the cosmological distance measurements in the FIDUCIAL
data set, we then marginalise over the parameter $f_{\rm v}$ from
$[-0.02,0.02]$.

##  4 Results and Discussions

###  4.1 Distance measurements without cosmological model assumptions

Figure 3: Measurements from our joint strong lensing and stellar dynamical
models, without any kinematic data (grey) and with mock IFU stellar kinematics
(orange and blue). The corner plot shows the most relevant parameters in the
fit, consisting of the lens distance $D_{\mathrm{d},{\rm int}}$, time-delay
distance $D_{{\Delta t},{\rm int}}$, internal mass-sheet parameter
$\lambda_{\mathrm{int}}$ and orbital anisotropy $\beta_{z}$. The contours
highlight the 1-, 2- and 3$\sigma$ confidence regions. The black points
(lines) depict the mock input values. Without any kinematics, $D_{{\Delta
t},{\rm int}}$ and $\lambda_{\mathrm{int}}$ are highly degenerate while
$D_{\mathrm{d},{\rm int}}$ and $\beta_{z}$ are unconstrained. The posterior
PDFs for $D_{\mathrm{d},{\rm int}}$ and $\beta_{z}$ are virtually flat in this
case, which is why we only show the 2D contours in the $D_{{\Delta t},{\rm
int}}$-$\lambda_{\mathrm{int}}$ space. IFU data will considerably tighten the
constraints. In the IDEAL case, the mock input values are recovered within
2$\sigma$. In the FIDUCIAL case, the mock input values are also recovered
within the 2$\sigma$ bands, but with larger uncertainties for the lens
distance ($\delta D_{\mathrm{d},{\rm int}}/D_{\mathrm{d},{\rm int}}\approx
4\%$). These measurements are independent of cosmological model assumptions.

The probability distributions for the various data sets are displayed in Fig.
3. As expected, the introduction of a constant internal mass sheet is
degenerate with the mass model. By scaling the mass model parameters
accordingly, the lensing data (consisting of high-resolution imaging and time
delays) can easily be recovered, while allowing for a large variation in the
contribution of $\lambda_{\mathrm{int}}$. Consequently, we observe a posterior
probability distribution function (PDF) for $\lambda_{\mathrm{int}}$ and
$D_{{\Delta t},{\rm int}}$ (grey) which is essentially flat within the prior
bounds. The PDF for the lens distance $D_{\mathrm{d},{\rm int}}$ and the
anisotropy parameter $\beta_{z}$ is also flat, since these are anchored by the
kinematic data only. It is worth noting here, that the effect depicted in Fig.
3 for the lensing-only data is neither new nor unexpected, but simply a
manifestation of the MSD. The degeneracy is well captured by our models and
this exercise also serves the purpose of illustrating the differences in the
constraining power of a mass sheet with and without high-spatially resolved
stellar kinematic data.

The corresponding probability distributions for our joint lensing and
dynamical models are also shown in Fig 3, where the IDEAL (orange) and
FIDUCIAL (blue) mock IFU observations were employed. In contrast to the highly
degenerate lensing-only run, an internal mass-sheet component is constrained
when IFU stellar kinematics are available. We observe tight constraints for
any contribution from an internal mass sheet with a core radius of 20″ of
$\lambda_{\mathrm{int}}=0.97^{+0.03}_{-0.05}$. The sensitivity of the 2D
kinematics to changes in the radial mass distribution translates to tight
measurements on the distances, too. We measure $D_{\mathrm{d},{\rm int}}$ and
$D_{{\Delta t},{\rm int}}$ with a precision of $\sim$ 1.6% and $\sim$ 3.7%
respectively (see Table 2) in the IDEAL case. Yet, the constraint for the lens
distance degrades significantly when the FIDUCIAL data is employed instead,
where $D_{\mathrm{d},{\rm int}}$ is only measured with a $\sim$ 4% precision.
In both cases, however, the mock input values are generally recovered within
the 1-2$\sigma$ uncertainty bands. The results can be explained as follows;
the statistical uncertainty in the mock kinematic data ($\delta
v_{\rm{stat}}$) sets the overall precision for measuring
$\lambda_{\mathrm{int}}$ and $D_{{\Delta t},{\rm int}}$, whereas the
systematic offset ($\delta v_{\rm{sys}}$) dictates the best precision to be
achieved for $D_{\mathrm{d},{\rm int}}$. Hence, with the same statistical
uncertainties of 2.5% in both the IDEAL and FIDUCIAL data, the final modelling
uncertainties for $\lambda_{\mathrm{int}}$ and $D_{{\Delta t},{\rm int}}$ are
identical (see Fig. 3), while the uncertainties on $D_{\mathrm{d},{\rm int}}$
in the FIDUCIAL case can readily be obtained from the IDEAL PDF, by utilising
the analytic relationship in Eq. 20888In detail, the final uncertainties for
$D_{\mathrm{d},{\rm int}}$ in the FIDUCIAL case are obtained by marginalising
over models with various realisations of the systematic offset $f_{\rm v}$. In
practice, we create 9 realisations corresponding to the values of $f_{\rm
v}=\\{-0.02,-0.015,-0.01,-0.005,0.,0.005,0.01,0.015,0.02\\}$, and combine the
resulting chains from these realisations with equal weights.. This yields the
final PDF in Fig. 3.

While we recover the key parameters of our mass models within the 1-2$\sigma$
regions as shown in Fig. 3, we note here some of the effects that affect the
shape of the posterior PDF of the model parameters. First of all, the amount
of noise in the mock data sets the size of the posterior PDF – the higher the
amount of noise, the broader the PDF (as further illustrated in Appendix A).
Secondly, the posterior distribution of $\lambda_{\mathrm{int}}$ is highly
asymmetric due to the imposed physical condition of positive mass density.
Finally, models with different source pixelisations lead to slightly different
optimisations, resulting in shifts of the PDF from individual source
pixelisations relative to each other (Suyu et al., 2013, Y20). The final PDF
is obtained by combining the individual PDFs for each source resolution,
giving equal weighting to all due to indistinguishable kinematic goodnes of
fit values ($\chi_{\rm D}$), in order to account for the model uncertainties
due to source pixelisation. This is in contrast to Y20, where the Bayesian
Information Criterion (BIC; Schwarz, 1978) was utilised effectively to
downweight individual source resolutions that are too far away from the mock
input values, but no longer applicable in the presence of a MST to further
constrain this source of uncertainty in the models.

For the orbital anisotropy, we picked a mock input value of $\beta_{z}=0.15$.
Our choice of mild radial anisotropy was motivated by the findings from
extensive modelling of massive elliptical galaxies in the local universe
(Cappellari et al., 2007, 2013). Regardless of this particular choice, Fig. 3
is a powerful illustration of how time delays and high-spatially resolved
kinematic data will enable us to lift the MSD and mass-anisotropy degeneracy
(van der Marel & Franx, 1993; Gerhard, 1993) simultaneously, providing tight
constraints also on the orbital anisotropy parameter, which is otherwise
difficult unless higher order moments of the line-of-sight velocity
distribution are available.

The anisotropy in our mock data and models was, so far, assumed to be constant
as a function radius. We also explored the scenario of allowing for spatially
varying anisotropy profiles to assess its impact on the inference of the
cosmological distances in the fit. We present the results in Fig. 4, where we
compare the constraints for the two different parameterisations of the orbital
anisotropy profile. Models which adopt a constant value for the anisotropy fit
the data exceptionally well and recover the mock input values within the
1$\sigma$ confidence levels (indigo). This is identical to Fig. 3, but shown
here for a single source resolution of $64\times\ 64$ pixels. Models with a
radially varying anisotropy profile have been implemented, where the
anisotropy profile is parameterised according to:

$\rm{log_{10}}\,(\beta_{z}(r))=\rm{log_{10}}\,(\alpha)+\beta\times\rm{log_{10}}\,(r),$
(21)

with $\alpha\in[-0.3,0.3]$ and $|\,\beta\,|\in[0,0.5]$. This parameterisation
provides enough freedom to allow for a wide range of anisotropy profiles,
while ensuring a smooth transition from the central to the outermost bins. It
is also tailored to account for the coarser sampling in the remote regions.
Models with a radially varying anisotropy profile are capable of recovering
the mock input value within 1$\sigma$ (teal), with $\alpha=0.13\pm 0.02$ and
$\beta=0.0\pm 0.1$. Despite the increased freedom in the anisotropy profile,
we do not observe inflated errors for the cosmological distances in the fit,
showing that the orbit distribution has limited power for compensating the
effects of the distances and mass-sheet, as illustrated in Fig. 2.

Figure 4: Measurements from our joint strong lensing and stellar dynamical
models with mock IFU stellar kinematics (IDEAL), but for two different
parameterisations of the orbital anisotropy profile $\beta_{z}$. The corner
plot shows the relevant parameters in the fit, consisting of the lens distance
$D_{\mathrm{d},{\rm int}}$, time-delay distance $D_{{\Delta t},{\rm int}}$ and
internal mass-sheet parameter $\lambda_{\mathrm{int}}$. The mock data is
mildly anisotropic (i.e., $\beta_{z}(r)=0.15$ is non-zero and independent of
radius) and fitted with both a radially constant (indigo) and radially varying
anisotropy profile (teal). The contours highlight the 1-, 2- and 3$\sigma$
confidence regions, with the black points (lines) depicting the mock input
values. The radially constant mock input anisotropy can be recovered from
models with a radially varying anisotropy profile. Table 2: Parameter
constraints from our lensing-only and joint strong lensing & stellar dynamical
models. The IFU stellar kinematics have been mocked with (FIDUCIAL) and
without (IDEAL) systematic errors included and a statistical uncertainty of
$\delta v_{\rm{stat}}$ = 2.5% per bin. The last column represents the inferred
$H_{\mathrm{0}}$ values for the various models, under the assumption of flat
$\Lambda$CDM. The mock input value is 82.5 km s-1 Mpc-1. We quote the 50th and
16th/84th percentiles of the distribution.

Model $D_{\mathrm{d},{\rm int}}$ [Mpc] $D_{{\Delta t},{\rm int}}$ [Mpc]
$\lambda_{\mathrm{int}}$ $H_{\mathrm{0}}$ [km s-1 Mpc-1] (flat $\Lambda$CDM)
Lensing - 1749${}^{+897}_{-891}$ 1.00${}^{+0.50}_{-0.50}$ 82.7${}^{+24}_{-22}$
Lensing & Dynamics 761${}^{+12}_{-12}$ 1816${}^{+81}_{-53}$
0.97${}^{+0.03}_{-0.05}$ 84.1${}^{+2.8}_{-2.6}$ IDEAL with $\delta v_{\rm
stat}=2.5\%$ Lensing & Dynamics 791${}^{+26}_{-26}$ 1816${}^{+81}_{-53}$
0.97${}^{+0.03}_{-0.05}$ 82.2${}^{+3.1}_{-3.1}$ FIDUCIAL with $\delta v_{\rm
stat}=2.5\%$

We further highlight the difference in the 2D kinematics between the IDEAL
mock data, the best-fitting model without a mass-sheet
($\lambda_{\mathrm{int}}\approx 1.0$) and the best-fitting model with a
(fixed) mass-sheet of $\lambda_{\mathrm{int}}\approx 0.85$ in Fig. 5. The mock
data is exceptionally well fitted without a mass sheet and generally recovered
within the 1-2$\sigma$ errors. A value of $\lambda_{\mathrm{int}}\leq 0.9$ is
deemed necessary to align the H0LiCOW $H_{\mathrm{0}}$ constraints of
RXJ1131$-$1231 with those from the cosmic microwave background (B20), but such
significant deviations from lens galaxies without a sheet can be ruled out at
the $\sim 3\sigma$ level, based on our IDEAL mock kinematics. The difference
between the overall best-fitting model without a mass sheet (i.e.
$\lambda_{\mathrm{int}}=1.0$) and the best-fitting model with a significant
mass-sheet contribution of $\lambda_{\mathrm{int}}\approx 0.85$ is
$\Delta\chi^{2}=6$. Given the same priors for the two models above, the ratio
of the posterior probabilities is determined by the Bayes factor, which we
approximate by means of the BIC difference (Y20), and which yields a posterior
odd of 0.05 for the model with a mass sheet against the model without a mass
sheet.

Figure 5: First: IDEAL mock IFU $\overline{v_{\mathrm{LOS}}^{2}}$ data of
RXJ1131$-$1231, at JWST NIRSpec resolution. Following the procedure in Y20,
the data has been mocked up by means of the best-fitting lensing-only model
and a set of random dynamical parameters. Second: normalised residuals (in
absolute values) of the best-fitting model without a mass-sheet component
($\lambda_{\mathrm{int}}\approx 1.0$), showing that the mock data is
remarkably well recovered mostly within the 1-$\sigma$ uncertainties. Third:
normalised residuals (in absolute values) of the best-fitting model with a
(fixed and significant) mass-sheet component ($\lambda_{\mathrm{int}}\approx
0.85$), which is at the edge of the 3-$\sigma$ uncertainty bands. Fourth: same
as the two panels before, but for the best-fitting model with a SPEMD. Various
mock data points are not recovered even within the 2-3$\sigma$ uncertainties.
The $\chi^{2}$ difference between this model and the best-fitting COMPOSITE
model without a sheet (second panel) is $\Delta\chi^{2}\approx 155$.

We also constructed models with a softened power-law elliptical mass
distribution (SPEMD; Barkana, 1998) and, again, added a PIEMD profile to mimic
the contribution of an internal mass-sheet. In agreement with our earlier
studies (Y20), the SPEMD provides a considerably worse fit to the mock
kinematics. The best-fitting SPEMD $\chi_{\rm D}^{2}\approx 235$, as opposed
to $\chi_{\rm D}^{2}\approx 78$ for the COMPOSITE mass model. Whereas the BIC
differences are insufficient to discriminate between the various source
resolutions of a given mass model, the BIC differences are sufficient to rule
out mass density parameterisations that deviate strongly from the input model,
even in light of the increased freedom provided by the MSD. Therefore,
including the SPEMD models in our final chains leaves the constraints for the
various parameters in the fit unchanged, given that these are assigned zero
weighting.

The final parameter constraints for $D_{\mathrm{d},{\rm int}}$, $D_{{\Delta
t},{\rm int}}$ and $\lambda_{\mathrm{int}}$ are summarised in Table 2. In
concordance with Y20, we observe tight measurements for $D_{\mathrm{d},{\rm
int}}$ and $D_{{\Delta t},{\rm int}}$ when IFU data is available. The
uncertainties on $D_{\mathrm{d},{\rm int}}$ are dominated by those of the
kinematic data, yet significantly better than expected from e.g. a single
aperture averaged velocity dispersion (Chen et al., 2020; Jee et al., 2019),
given the roughly 80 (almost) independent measurements of the LOSVD across the
entire FOV. The precision of $D_{{\Delta t},{\rm int}}$ also greatly benefits
from high-spatially resolved kinematics (Y20), with $|\delta D_{{\Delta
t},{\rm int}}/D_{{\Delta t},{\rm
int}}|\sim|\delta\lambda_{\mathrm{int}}/\lambda_{\mathrm{int}}|$, but is
highly sensitive to the statistical errors in the measurement of the LOSVD.
For illustration purposes, we also present an even more conservative scenario
in Appendix A, in which we increase $\delta v_{\rm{stat}}$ to 3.5% instead,
which further highlights the impact of $\delta v_{\rm{stat}}$ in constraining
$\lambda_{\mathrm{int}}$ and $D_{{\Delta t},{\rm int}}$.

###  4.2 Cosmographic forecast in flat $\Lambda$CDM

Unlike earlier findings without an internal MST, where much of the
constraining power in measuring $H_{\mathrm{0}}$ came from $D_{{\Delta t},{\rm
int}}$ only (with the small $\delta D_{{\Delta t},{\rm int}}$ propagating
almost linearly into the error budget on $H_{\mathrm{0}}$),
$D_{\mathrm{d},{\rm int}}$ is now expected to add substantial constraining
power to $D_{{\Delta t},{\rm int}}$ in measuring $H_{\mathrm{0}}$ (Jee et al.,
2016; Chen et al., 2020).

To this end, we utilise the joint distance measurements of $D_{\mathrm{d},{\rm
int}}$ and $D_{{\Delta t},{\rm int}}$ to obtain cosmological forecasts in flat
$\Lambda$CDM. We first fit the marginalised $D_{{\Delta t},{\rm
int}}$-$D_{\mathrm{d},{\rm int}}$ distribution (from all chains across various
source resolutions) with a kernel density estimator and account for mass
structures along the lens LOS with an external convergence
($\kappa_{\mathrm{ext}}$) distribution, following the number counts approach
from numerical simulations in Suyu et al. (2013, 2014). We draw random samples
of {$H_{\mathrm{0}}$, $\Omega_{\rm m},\kappa_{\mathrm{ext}}\\}$ with uniform
priors on $H_{\mathrm{0}}$ of [50,120] km s-1 and on $\Omega_{\rm{m}}$ of
[0.05, 0.5] and with the $P(\kappa_{\mathrm{ext}})$ distribution, compute the
corresponding $D_{{\Delta t},{\rm int}}$ and $D_{\mathrm{d},{\rm int}}$
following Eqs. 7 and 8, and importance sample given the $P(D_{\mathrm{d},{\rm
int}},D_{{\Delta t},{\rm int}})$ posterior. This yields a final value of
$H_{\mathrm{0}}$ $=82.2^{+3.1}_{-3.1}$ km s-1 Mpc-1 ($H_{\mathrm{0}}$
$=84.1^{+2.8}_{-2.6}$ km s-1 Mpc-1), i.e. a 3.2% precision measurement (3.7%)
for the IDEAL (FIDUCIAL) case, recovering the input value of 82.5 km s-1
Mpc-1.

Fig. 6 shows how the inclusion of 2D stellar kinematics of the lens not only
helps in tightly constraining the lens distance, but also leads to improved
constraints on cosmological parameters beyond $H_{\mathrm{0}}$, even for a
single lens. This conclusion has already been drawn in Y20 and is now
corroborated even in the presence of a MST, assuming that the stellar
kinematics can be measured with high accuracy and precision. Otherwise, more
systems will be needed to provide more meaningful constraints beyond
$H_{\mathrm{0}}$.

Figure 6: $H_{\mathrm{0}}$ and $\Omega_{\rm{m}}$ constraints from our models
in flat $\Lambda$CDM cosmology, for the IDEAL (left) and FIDUCIAL (data set)
with statistical uncertainties of $\delta v_{\rm{stat}}$ = 2.5%. We adopted
uniform priors on $H_{\mathrm{0}}$ of [50, 120] km s-1 Mpc-1and
$\Omega_{\rm{m}}$ of [0.05, 0.5]. The green shaded regions show the 1, 2 and
3$\sigma$ confidence regions based on the constraint on $D_{\Delta t}$ only,
i.e., $P(D_{\Delta t})$ after marginalising over other mass model parameters
including $D_{\mathrm{d}}$. The red shaded regions show the corresponding
constraints based on $P(D_{\mathrm{d}})$, after marginalising over $D_{\Delta
t}$. The grey histogram shows the constraints based on
$P(D_{\mathrm{d}},D_{\Delta t})$, where we have omitted showing the 2D
contours for ease of viewing. Different tilts in the
$H_{\mathrm{0}}$-$\Omega_{\rm{m}}$ space, when marginalised over
$D_{\mathrm{d}}$ and $D_{\Delta t}$, respectively, can break some of the
degeneracies, leading to improved constraints also on $\Omega_{\rm{m}}$.

###  4.3 View of $D_{\mathrm{d}}$ and core size

It is of great importance to recapitulate our findings within a more
comprehensive and broader context: despite employing a model that is maximally
degenerate with $D_{{\Delta t},{\rm int}}$, and hence with $H_{\mathrm{0}}$
through the MST, we have shown that high-quality, high-spatially resolved
stellar kinematics can break the internal MSD and constrain $D_{{\Delta
t},{\rm int}}$. Furthermore, time-delays and (spatially resolved) kinematics
also provide leverage for constraining $D_{\mathrm{d}}$ (Paraficz & Hjorth,
2009) regardless of the MST (Jee et al., 2015, 2019; Chen et al., 2020).
Previous analyses in TDCOSMO using two families of mass models with single-
aperture velocity dispersion measurements and without MST yield joint
constraints on $D_{\Delta t}$ and $D_{\mathrm{d}}$, but with substantially
tighter constraints on $D_{\Delta t}$ compared to $D_{\mathrm{d}}$. When
allowing for a mass-sheet transformation, the precision on $D_{\Delta t}$
degrades. On the other hand, $D_{\mathrm{d}}$ is nearly insensitive to the
internal MST (and completely insensitive to the external MST). With high-
quality spatially-resolved kinematics, the precision on $D_{\mathrm{d}}$ can
become better than that of $D_{\Delta t}$ (Sec. 4.1, 4.2 and Y20). Therefore,
the $D_{\Delta t}$-centric view in TDSL can be complemented by the new
$D_{\mathrm{d}}$-centric view; namely, the combination of the two distances
provide robust $H_{0}$ even in the presence of the MST.

The results presented here rest upon the assumption of an extended, yet finite
core size of $\theta_{\rm c}=20\arcsec$, but it may be useful to emphasise
again the role played by $D_{\mathrm{d}}$ in the two limiting cases of
$\theta_{\rm c}\ll\theta_{\rm E}$ and $\theta_{\rm c}\gg\theta_{\rm E}$.
First, when $\theta_{\rm c}\ll\theta_{\rm E}$, the approximation of a pure MST
is no longer valid (Birrer et al., 2020) and lensing data are sensitive to
transformations of the lens mass profile, particularly the
$\lambda_{\mathrm{int}}$ values. In this scenario, the precision of $D_{\Delta
t}$ is mainly limited by the precision of the lensing and time delay
observations, whereas the precision of $D_{\mathrm{d}}$ mainly depends on the
time-delays and kinematics. Both $D_{\Delta t}$ and $D_{\mathrm{d}}$ can be
measured almost equally well in this case (see also Y20), thus leading to
tight constraints on the inferred cosmological parameters. Second, when
$\theta_{\rm c}\gg\theta_{\rm E}$, the kinematic data no longer constrains the
internal mass sheet, as the corresponding density profile effectively vanishes
(Appendix C). Under these circumstances, no improvement in the measurement of
$D_{\Delta t}$ and $H_{\mathrm{0}}$ can be expected, leaving it prone to the
full degeneracy introduced by the MST. However, the kinematics still improve
the constraint on $D_{\mathrm{d}}$ significantly on a single lens basis, which
then becomes the dominant source of information in the final cosmological
inference. While this was expected from previous studies on the subject (Jee
et al., 2015, 2016, 2019; Chen et al., 2020), our results demonstrate this
explicitly and quantitatively using mock JWST observations.

###  4.4 Comparison to previous TDCOSMO studies

We compare our approach to previous analyses of TDCOSMO, particularly TDCOSMO
I (Millon et al., 2020b), IV (Birrer et al., 2020), V (Birrer & Treu, 2021)
and VI (Chen et al., 2020) that involve lens mass modelling.We start by
summarising previous studies, and explain the similarities and differences.

In TDCOSMO I, Millon et al. (2020b) investigated potential sources of
systematic uncertainties connected to lens mass modelling by employing two
families of physically motivated mass models, the power-law and composite mass
models, without internal MST. They find the two families of models yield the
same $H_{0}$ within 1% uncertainty from the sample of lens systems that were
previously modeled by the H0LiCOW, SHARP and STRIDES collaborations (Suyu et
al., 2013, 2014; Wong et al., 2017; Birrer et al., 2019; Chen et al., 2019;
Rusu et al., 2020; Shajib et al., 2020). The comparison of these two families
of model demonstrates that there is no evidence for significant residual
systematic uncertainties associated with the lens mass modelling.

To further explore the MST as a source of uncertainty for measuring $H_{0}$,
Birrer et al. (2020) in TDCOSMO IV relaxed the assumptions on the lens mass
profiles and allowed for an internal MST of the power-law model. The MST
parameter $\lambda_{\mathrm{int}}$ is then constrained through spherical Jeans
modelling (with a radially varying anisotropy profile) and the stellar
kinematic data of the lens galaxy, which are mostly single-aperture averaged
stellar velocity dispersions. A subset of the SLACS lens galaxies is further
combined with the 7 TDCOSMO lenses through a hierarchical analysis to
constrain $\lambda_{\mathrm{int}}$ and infer $H_{0}$ in flat $\Lambda$CDM. The
uncertainties on $H_{0}$ are 8% from TDCOSMO lenses only, or 5% from TDCOSMO
and SLACS lenses. Birrer & Treu (2021) in TDCOSMO V then estimated the
precisions of $H_{0}$ measurements from current and future samples of time-
delay lenses, with various levels of velocity dispersion measurements, ranging
from single aperture-averaged velocity dispersions to IFU kinematic maps with
30 radially averaged velocity bins. With 7 TDCOSMO lenses and future velocity
dispersion measurements in 10 radial bins for each lens, Birrer & Treu (2021)
anticipate a 3.3% precision on $H_{0}$ in flat $\Lambda$CDM. An expansion of
the sample to 40 time-delay lenses and 200 non-time-delay lenses would improve
the precision to 1.5% and 1.2%, respectively.

The analyses of TDCOSMO IV and V assumed the flat $\Lambda$CDM model, which
reduces the range of plausible $\lambda_{\mathrm{int}}$ in the MST given the
tight correlation between $D_{\mathrm{d}}$ and $D_{\Delta t}$. Chen et al.
(2020) in TDCOSMO VI relaxed the assumption of flat $\Lambda$CDM by
considering the measurements of $D_{\Delta t}$ and $D_{\mathrm{d}}$ in generic
cosmological models, while allowing for the MST. By combining the lenses
together with relative distance indicators such as Type Ia Supernovae and
Baryon Acoustic Oscillations, Chen et al. (2020) showed that
$\lambda_{\mathrm{int}}$ and $D_{\Delta t}$ can be constrained in generic
cosmological models. In this work, spherical Jeans modelling and aperture-
averaged kinematic measurements are used.

Our current paper is similar to TDCOSMO IV, V, VI in that we have also
considered a parametrised profile with an internal mass-sheet transformation
characterised by $\lambda_{\mathrm{int}}$. While TDCOSMOS IV, V, VI used only
the power-law profile with mass sheets, we have considered both families of
mass models (composite and power-law profiles that were used in TDCOSMO I)
with additional mass sheets. We find that spatially resolved kinematics not
only allow us to constrain the MST parameter $\lambda_{\mathrm{int}}$, but
also to distinguish between the families of mass models. Our set up is similar
to TDCOSMO VI in that we aim to measure the distances $D_{\Delta t}$ and
$D_{\mathrm{d}}$ without assuming a background cosmological model, in contrast
to TDCOSMO IV and V that assumed flat $\Lambda$CDM.

The key difference and novelty of our current work are the use of a more
flexible, self-consistent lensing and dynamical mass model. Instead of
spherical Jeans modelling for the kinematics that were used in previous
TDCOSMO studies, we have employed an axisymmetric mass distribution (in 3D)
for both the lensing and kinematic analysis. Furthermore, we have considered
2D spatially resolved kinematic data with angular structure, which help to
constrain both $\lambda_{\mathrm{int}}$ and $D_{\Delta t}$, even when allowing
for radially varying anisotropy profiles. With our forecast of JWST-like
spatially resolved kinematics of $\sim$80 bins for RXJ1131$-$1231, we find
that we can constrain $H_{0}$ with this single lens with better than 4%
precision. This is based on the assumption of a uniformly distributed
systematic uncertainty of 2% with $f_{\rm v}$ varying between [-0.02,0.02],
although we show in Appendix B that a Gaussian distributed systematic
uncertainty of 2% increases the $H_{0}$ uncertainty only slightly from 3.8% to
4.1% for our mock kinematic data with $\delta v_{\rm{stat}}$ = 2.5%. In
comparison to TDCOSMO V that adopted a 3% systematic uncertainty that is
larger than our 2%, our forecast of $\leq 4\%$ uncertainty on $H_{0}$
constrained from a single lens is more precise, as expected. Once we account
for the difference in the assumed kinematic data quality, the two studies
(this and TDCOSMO V) should provide consistent results. A direct comparison of
the two dynamical modeling approaches (axisymmetric Jeans in this paper and
spherical Jeans in TDCOSMO V) would require modeling the same kinematic,
lensing and time delay data set, which is not the focus of our current paper.

One other key difference to TDCOSMO IV, V and VI is that with the anticipated
high-quality JWST kinematic data, we can constrain the parameters for the mass
sheet ($\lambda_{\mathrm{int}}$) and for the anisotropy of the stellar orbits
from the JWST kinematic data, without having to rely on external samples of
non-time-delay lenses (in TDCOSMO IV and V) and/or other relative distance
indicators (in TDCOSMO VI). We therefore see the individual analyses of lenses
with high-quality data to be complementary to the hierarchical analyses with
external samples – the combination of both approaches in future studies would
give us the maximal amount of cosmological information especially in measuring
$H_{0}$.

##  5 Summary

In order to assess the constraining power on a mass-sheet component in TDSL
studies, we have performed a suite of simulations. For this purpose, we used
the prominent lens system RXJ1131$-$1231 as a test-bed for our studies. We
mocked up high-quality, high-spatially resolved IFU stellar kinematics (as is
expected from the next generation of ground and space-based telescopes)
together with mock HST imaging, and used long-baseline time-delay
measurements. Following the prescription of a cored density distribution
(B20), a variable internal mass-sheet component was included in our models.
Our findings can be summarised as follows:

  1. 1.

Without any kinematic data (i.e. for pure strong lens models), we observe the
well known degeneracy between the mass model parameters and any internal mass
sheet component $\lambda_{\mathrm{int}}$. The lensing data allow for a wide
range of $\lambda_{\mathrm{int}}$ values, which can adversely affect the time-
delay distance measurement $D_{{\Delta t},{\rm int}}$, according to Eq. 6.

  2. 2.

Future IFU stellar kinematics with the next generation of telescopes (such as
JWST and E-ELT) will be able to effectively break the internal MSD, by placing
tight constraints on the contribution of any internal mass-sheet component
with a core radius $\lessapprox 20\arcsec$ (This finding is specific to
RXJ1131$-$1231, However, with $\theta_{\rm E}/R_{e}\approx 0.9$,
RXJ1131$-$1231 is not a particular outlier in the H0LiCOW sample (Millon et
al., 2020b).).

  3. 3.

Using realistic mock IFU observations, provided by assuming reasonable on-
source integration times with JWST, we show that future measurements of
$\lambda_{\mathrm{int}}$ (and hence $D_{{\Delta t},{\rm int}}$) could be as
precise as 4% per lens. These uncertainties propagate almost linearly into the
error budget on $H_{\mathrm{0}}$, whereas the precise $D_{\mathrm{d},{\rm
int}}$ measurements ($\delta D_{\mathrm{d},{\rm int}}/D_{\mathrm{d},{\rm
int}}\leq 4\%$) allows further reduction in the $H_{\mathrm{0}}$
uncertainties. This yields a high-precision measurement on $H_{\mathrm{0}}$ of
$\leq 4\%$ from this single lens system, assuming flat $\Lambda$CDM cosmology.

  4. 4.

High-spatially resolved stellar kinematics of the lens can constrain
cosmological parameters beyond $H_{\mathrm{0}}$. However, more systems are
needed to provide more meaningful information on e.g. $\Omega_{m}$.

  5. 5.

The modelling approach presented here is fully self-consistent and more
flexible than literature studies where spherical symmetry is commonly employed
for modelling the stellar kinematic data. The 2D angular structure in the
kinematic data together with the lensing information allow us to recover the
anisotropy parameters of the stellar orbits, even with radial variations in
the stellar anisotropy parametrisation.

  6. 6.

We do not rely and make no use of informative priors at the galaxy population
level. This approach presents an alternative and complementary method to that
presented in Birrer et al. (2020), to gain back much of the precision lost by
accounting for the MSD. In addition, spatially-resolved kinematics can be
used, in turn, to verify if external data sets can be combined with the
TDCOSMO lenses from the galaxy population point of view.

  7. 7.

Our final inference of the angular diameter distances $D_{\mathrm{d},{\rm
int}}$ and $D_{{\Delta t},{\rm int}}$ are independent of cosmological model
assumptions.

It is beyond the scope of this paper to speculate on the physical plausibility
and origin of such a mass sheet999The internal mass-sheet does not have to be
a physical sheet of mass in the strictest terms, but could also reflect minor
deviations from the actual input mass model.. But, our analysis shows that
upper limits on its spatial extent can be provided by means of JWST-like
stellar kinematics, unless $r_{\rm c}\gtrapprox 100$ kpc. This will be helpful
in further constraining physical models that are put forward for producing
cored density distributions. Given this forecast, much of the precision lost
by incorporating a maximally degenerate model can be regained, such that only
a few systems with similar precision will be sufficient to address the $H_{0}$
tension between late and early time cosmological probes. As new IFU facilities
come online in the next years, we advocate for utilising this unique
combination of lensing and dynamics, particularly with spatially resolved
kinematics, as a powerful probe of cosmology.

## Acknowledgements

The authors thank S. Birrer, F. Courbin, D. Sluse and T. Treu for helpful
discussions and comments on the manuscript. AY and SHS thank the Max Planck
Society for support through the Max Planck Research Group for SHS. G C.-F. C
acknowledges support from the National Science Foundation through grants NSF-
AST-1906976 and NSF-AST-1836016. G. C.-F. C. acknowledges support from the
Moore Foundation through grant 8548. SHS and EK are supported in part by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311. The Kavli IPMU is
supported by World Premier International Research Center Initiative (WPI),
MEXT, Japan.

## References

  * Abbott et al. (2018) Abbott, T. M. C., Abdalla, F. B., Annis, J., et al. 2018, MNRAS, 480, 3879
  * Agrawal et al. (2019) Agrawal, P., Cyr-Racine, F.-Y., Pinner, D., & Randall, L. 2019, arXiv e-prints, arXiv:1904.01016
  * Barkana (1998) Barkana, R. 1998, ApJ, 502, 531
  * Barnabè et al. (2012) Barnabè, M., Dutton, A. A., Marshall, P. J., et al. 2012, MNRAS, 423, 1073
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S. 1987, Galactic Dynamics
  * Birrer et al. (2020) Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, A165
  * Birrer & Treu (2021) Birrer, S. & Treu, T. 2021, A&A, 649, A61
  * Birrer et al. (2019) Birrer, S., Treu, T., Rusu, C. E., et al. 2019, MNRAS, 484, 4726
  * Blum et al. (2020) Blum, K., Castorina, E., & Simonović, M. 2020, ApJ, 892, L27
  * Blum & Teodori (2021) Blum, K. & Teodori, L. 2021, arXiv e-prints, arXiv:2105.10873
  * Cappellari (2002) Cappellari, M. 2002, MNRAS, 333, 400
  * Cappellari (2008) Cappellari, M. 2008, MNRAS, 390, 71
  * Cappellari (2020) Cappellari, M. 2020, MNRAS, 494, 4819
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345
  * Cappellari et al. (2007) Cappellari, M., Emsellem, E., Bacon, R., et al. 2007, MNRAS, 379, 418
  * Cappellari et al. (2013) Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709
  * Chen et al. (2019) Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743
  * Chen et al. (2020) Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2020, arXiv e-prints, arXiv:2011.06002
  * Emsellem et al. (1994) Emsellem, E., Monnet, G., & Bacon, R. 1994, A&A, 285, 723
  * Falco et al. (1985) Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Freedman et al. (2019) Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34
  * Gerhard (1993) Gerhard, O. E. 1993, MNRAS, 265, 213
  * Jeans (1922) Jeans, J. H. 1922, MNRAS, 82, 122
  * Jee et al. (2015) Jee, I., Komatsu, E., & Suyu, S. H. 2015, J. Cosmology Astropart. Phys., 11, 033
  * Jee et al. (2016) Jee, I., Komatsu, E., Suyu, S. H., & Huterer, D. 2016, J. Cosmology Astropart. Phys., 4, 031
  * Jee et al. (2019) Jee, I., Suyu, S. H., Komatsu, E., et al. 2019, Science, 365, 1134
  * Kassiola & Kovner (1993) Kassiola, A. & Kovner, I. 1993, ApJ, 417, 450
  * Knox & Millea (2020) Knox, L. & Millea, M. 2020, Phys. Rev. D, 101, 043533
  * Kochanek (2020) Kochanek, C. S. 2020, MNRAS, 493, 1725
  * Millon et al. (2020a) Millon, M., Courbin, F., Bonvin, V., et al. 2020a, A&A, 640, A105
  * Millon et al. (2020b) Millon, M., Galan, A., Courbin, F., et al. 2020b, A&A, 639, A101
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Paraficz & Hjorth (2009) Paraficz, D. & Hjorth, J. 2009, A&A, 507, L49
  * Pesce et al. (2020) Pesce, D. W., Braatz, J. A., Reid, M. J., et al. 2020, ApJ, 891, L1
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Poulin et al. (2019) Poulin, V., Smith, T. L., Karwal, T., & Kamionkowski, M. 2019, Phys. Rev. Lett., 122, 221301
  * Refsdal (1964) Refsdal, S. 1964, MNRAS, 128, 307
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85
  * Rusu et al. (2020) Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440
  * Schneider & Sluse (2013) Schneider, P. & Sluse, D. 2013, A&A, 559, A37
  * Schneider & Sluse (2014) Schneider, P. & Sluse, D. 2014, A&A, 564, A103
  * Schwarz (1978) Schwarz, G. 1978, Ann. Statist., 6, 461
  * Shajib et al. (2020) Shajib, A. J., Birrer, S., Treu, T., et al. 2020, MNRAS, 494, 6072
  * Shanks et al. (2019) Shanks, T., Hogarth, L. M., & Metcalfe, N. 2019, MNRAS, 484, L64
  * Sluse et al. (2003) Sluse, D., Surdej, J., Claeskens, J.-F., et al. 2003, A&A, 406, L43
  * Sonnenfeld (2018) Sonnenfeld, A. 2018, MNRAS, 474, 4648
  * Suyu et al. (2013) Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70
  * Suyu et al. (2017) Suyu, S. H., Bonvin, V., Courbin, F., et al. 2017, MNRAS, 468, 2590
  * Suyu et al. (2018) Suyu, S. H., Chang, T.-C., Courbin, F., & Okumura, T. 2018, Space Sci. Rev., 214, 91
  * Suyu et al. (2010) Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201
  * Suyu et al. (2006) Suyu, S. H., Marshall, P. J., Hobson, M. P., & Blandford, R. D. 2006, MNRAS, 371, 983
  * Suyu et al. (2014) Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35
  * Tewes et al. (2013a) Tewes, M., Courbin, F., & Meylan, G. 2013a, A&A, 553, A120
  * Tewes et al. (2013b) Tewes, M., Courbin, F., Meylan, G., et al. 2013b, A&A, 556, A22
  * Treu & Marshall (2016) Treu, T. & Marshall, P. J. 2016, A&A Rev., 24, 11
  * Unruh et al. (2017) Unruh, S., Schneider, P., & Sluse, D. 2017, A&A, 601, A77
  * van de Ven et al. (2010) van de Ven, G., Falcón-Barroso, J., McDermid, R. M., et al. 2010, ApJ, 719, 1481
  * van der Marel & Franx (1993) van der Marel, R. P. & Franx, M. 1993, ApJ, 407, 525
  * Verde et al. (2019) Verde, L., Treu, T., & Riess, A. G. 2019, Nature Astronomy, 3, 891
  * Wertz et al. (2018) Wertz, O., Orthen, B., & Schneider, P. 2018, A&A, 617, A140
  * Wong et al. (2017) Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420
  * Yıldırım et al. (2020) Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, MNRAS, 493, 4783

## Appendix A

We present the distance constraints and cosmological forecast from our joint
strong lensing and stellar dynamical models, based on the IDEAL and FIDUCIAL
data. Unlike the results shown in Fig. 3 and Fig. 6, the statistical errror
($\delta v_{\rm stat}$) in the mock kinematics (Eq. 1) has been increased from
2.5% to 3.5%, while $\delta v_{\rm sys}$ remains $-2$% for the FIDUCIAL data.
This roughly corresponds to a target S/N of 30 per bin. As before, we assume a
uniform distribution of $f_{\rm v}\in[-0.02,0.02]$. Fig. 7 effectively
illustrates how the increased statistical uncertainty in the kinematics
translates into larger modelling uncertainties on $D_{{\Delta t},{\rm int}}$
and $\lambda_{\mathrm{int}}$, whereas the uncertainty on $D_{\mathrm{d},{\rm
int}}$ depends mainly on the fractional offset $f_{\rm v}$ distribution that
accounts for the systematic uncertainties in the kinematic measurements. Given
the degraded mock kinematics, the uncertainties on $D_{{\Delta t},{\rm int}}$
(further amplified by the degeneracy between $D_{{\Delta t},{\rm int}}$ and
$\lambda_{\mathrm{int}}$) become significant, such that most of the
constraining power in the cosmological inference hinges on the accuracy and
precision of $D_{\mathrm{d},{\rm int}}$ (Fig.8). This is most obvious when
drawing a direct comparison between the IDEAL and FIDUCIAL cosmological
forecast here. We obtain a final $H_{\mathrm{0}}$ measurement of
$82.2^{+3.0}_{-2.8}$ km s-1 Mpc-1 in the IDEAL case, whereas the FIDUCIAL data
yields $H_{0}=79.3^{+3.6}_{-3.6}$ km s-1 Mpc-1. The shift in the mean value
for the FIDUCIAL data is mainly caused by the significant drop in the
constraining power of $D_{{\Delta t},{\rm int}}$. With the systematic offset
in the velocity moments for the FIDUCIAL case, the mean $H_{\mathrm{0}}$ is
also biased. This is in contrast to our findings without a mass sheet in Y20
(see Table 3) where, despite increased statistical errors in the kinematic
data, both $D_{{\Delta t},{\rm int}}$ and $D_{\mathrm{d},{\rm int}}$
contribute almost equally to an accurate and precise $H_{\mathrm{0}}$
measurement. We summarise the results in Table 3.

Figure 7: Measurements from our joint strong lensing and stellar dynamical
models, without any kinematic data (grey) and with mock IFU stellar kinematics
(orange and blue) and an increased statistical error of $\delta v_{\rm
stat}=3.5\%$ for the kinematics. In concordance with Fig. 3, the plot shows
the constraints for the most relevant parameters in the fit, consisting of the
lens distance $D_{\mathrm{d},{\rm int}}$, time-delay distance $D_{{\Delta
t},{\rm int}}$, internal mass-sheet parameter $\lambda_{\mathrm{int}}$ and
orbital anisotropy $\beta_{z}$. The contours highlight the 1-, 2- and
3$\sigma$ confidence regions. The black points (lines) depict the mock input
values. The increased error of 3.5% for each velocity bin translates into
larger uncertainties on $D_{{\Delta t},{\rm int}}$ and
$\lambda_{\mathrm{int}}$. The uncertainties on $D_{\mathrm{d},{\rm int}}$ are
dictated by the systematic offset ($\delta v_{\rm sys}$) and are unchanged
with respect to the models in Sec. 4.1, given the same bias of $-2$% in
$\overline{v_{\mathrm{LOS,\,data,\,\textit{l}}}^{2}}$.

Figure 8: $H_{\mathrm{0}}$ and $\Omega_{\rm{m}}$ constraints from our models
in flat $\Lambda$CDM cosmology, for the IDEAL (left) and FIDUCIAL (data sets),
with increased statistical errors of $\delta v_{\rm stat}=3.5\%$. We adopted
uniform priors on $H_{\mathrm{0}}$ of [50, 120] km s-1 Mpc-1and on
$\Omega_{\rm m}$ of [0.05, 0.5]. The green shaded regions show the 1, 2 and
3$\sigma$ confidence regions based on the constraint on $D_{\Delta t}$ only,
i.e., $P(D_{\Delta t})$ after marginalising over other mass model parameters
including $D_{\mathrm{d}}$. The red shaded regions show the corresponding
constraints based on $P(D_{\mathrm{d}})$, after marginalising over $D_{\Delta
t}$. The grey histogram shows the constraints based on
$P(D_{\mathrm{d}},D_{\Delta t})$, where we have omitted showing the 2D
contours for ease of viewing. In principle, the different tilts in the
$H_{\mathrm{0}}$-$\Omega_{\rm m}$ space, when marginalised over
$D_{\mathrm{d}}$ and $D_{\Delta t}$ respectively, can break some of the
degeneracies. However, this requires highly accurate and precise measurements
of the LOSVD. Errors of $\delta v_{\rm stat}=3.5\%$ are already insufficient
to provide powerful constraints on e.g. $\Omega_{\rm m}$ on a single lens
basis when accounting for the MST. Table 3: Parameter constraints and
cosmological forecast from our lensing-only and joint strong lensing & stellar
dynamical models. The IFU stellar kinematics have been mocked with (FIDUCIAL)
and without (IDEAL) systematic errors included, but an increased statistical
error of $\delta v_{\rm stat}=3.5\%$. The last column represents the inferred
$H_{\mathrm{0}}$ values for the various models, under the assumption of flat
$\Lambda$CDM. The mock input value is 82.5 km s-1 Mpc-1. We quote the 50th and
16th/84th percentiles of the distribution.

Model $D_{\mathrm{d},{\rm int}}$ [Mpc] $D_{{\Delta t},{\rm int}}$ [Mpc]
$\lambda_{\mathrm{int}}$ $H_{\mathrm{0}}$ [km s-1 Mpc-1] (flat $\Lambda$CDM)
Lensing - 1749${}^{+897}_{-891}$ 1.00${}^{+0.50}_{-0.50}$ 82.7${}^{+24}_{-22}$
Lensing & Dynamics 782${}^{+15}_{-15}$ 1965${}^{+217}_{-153}$
0.89${}^{+0.09}_{-0.10}$ 82.2${}^{+3.0}_{-2.8}$ IDEAL with $\delta v_{\rm
stat}=3.5\%$ Lensing & Dynamics 813${}^{+28}_{-28}$ 1965${}^{+217}_{-153}$
0.89${}^{+0.08}_{-0.10}$ 79.3${}^{+3.6}_{-3.6}$ FIDUCIAL with $\delta v_{\rm
stat}=3.5\%$

## Appendix B

To facilitate a direct comparison with the cosmological forecast in Birrer &
Treu (2021), we perform the joint modelling of the strong lensing and stellar
kinematic data with $\delta v_{\rm stat}=2.5\%$ and $\delta v_{\rm sys}=2\%$
(orange) and $\delta v_{\rm stat}=3.5\%$ and $\delta v_{\rm sys}=2\%$ (blue),
and show the corresponding parameter constraints in Fig. 9. In contrast to
Section 4.1 where $f_{\rm v}$ was uniformly distributed between
$[-0.02,0.02]$, $f_{\rm v}$ is now assumed to be a Gaussian distribution with
a standard deviation of 0.02, for which we create 25 realisations of $f_{\rm
v}$, covering a range of [-0.06, 0.06] in steps of 0.005, and apply a
weighting of each chain given by $w=\rm exp\big{(}-0.5\,(f_{\rm v}/\sigma_{\rm
v})^{2}\big{)}$, with $\sigma_{\rm v}=0.02$. The corresponding cosmological
forecast for flat $\Lambda$CDM is presented in Fig. 10, with the results
summarised in Table 4.

Figure 9: Measurements from our joint strong lensing and stellar dynamical
models, without any kinematic data (grey) and with mock IFU stellar kinematics
(orange and blue). In both cases, we adopt $\delta v_{\rm sys}=-2\%$ and a
statistical error of $\delta v_{\rm stat}=2.5\%$ (orange) and $\delta v_{\rm
stat}=3.5\%$ (blue) respectively. In concordance with Fig. 3, the plot shows
the constraints for the most relevant parameters in the fit, consisting of the
lens distance $D_{\mathrm{d},{\rm int}}$, time-delay distance $D_{{\Delta
t},{\rm int}}$, internal mass-sheet parameter $\lambda_{\mathrm{int}}$ and
orbital anisotropy $\beta_{z}$. The contours highlight the 1-, 2- and
3$\sigma$ confidence regions. The black points (lines) depict the mock input
values. The increased error of 3.5% for each velocity bin translates into
larger uncertainties on $D_{{\Delta t},{\rm int}}$ and
$\lambda_{\mathrm{int}}$. The uncertainties on $D_{\mathrm{d},{\rm int}}$ are
dictated by the systematic offset ($\delta v_{\rm sys}$) and are almost
identical for both mock kinematic data sets, given the same bias of $-2$% in
$\overline{v_{\mathrm{LOS,\,data,\,\textit{l}}}^{2}}$.

Figure 10: $H_{\mathrm{0}}$ and $\Omega_{\rm{m}}$ constraints from our models
in flat $\Lambda$CDM cosmology, for the FIDUCIAL data sets with $\delta v_{\rm
stat}=2.5\%$ (left) and $\delta v_{\rm stat}=3.5\%$ (right). The systematic
uncertainty is $\delta v_{\rm sys}=-2.0\%$, but assuming a Gaussian
distribution for $f_{\rm v}$ within the range of [-0.06, 0.06] for both. We
adopted uniform priors on $H_{\mathrm{0}}$ of [50, 120] km s-1 Mpc-1and on
$\Omega_{\rm m}$ of [0.05, 0.5]. The green shaded regions show the 1, 2 and
3$\sigma$ confidence regions based on the constraint on $D_{\Delta t}$ only,
i.e., $P(D_{\Delta t})$ after marginalising over other mass model parameters
including $D_{\mathrm{d}}$. The red shaded regions show the corresponding
constraints based on $P(D_{\mathrm{d}})$, after marginalising over $D_{\Delta
t}$. The grey histogram shows the constraints based on
$P(D_{\mathrm{d}},D_{\Delta t})$, where we have omitted showing the 2D
contours for ease of viewing. In principle, the different tilts in the
$H_{\mathrm{0}}$-$\Omega_{\rm m}$ space, when marginalised over
$D_{\mathrm{d}}$ and $D_{\Delta t}$ respectively, can break some of the
degeneracies. However, this requires highly accurate and precise measurements
of the LOSVD. Errors of $\delta v_{\rm stat}=3.5\%$ are already insufficient
to provide powerful constraints on e.g. $\Omega_{\rm m}$ on a single lens
basis when accounting for the MST. Table 4: Parameter constraints and
cosmological forecast from our lensing-only and joint strong lensing & stellar
dynamical models. The IFU stellar kinematics have been mocked with systematic
errors included (FIDUCIAL), and two realisations for the statistical error of
$\delta v_{\rm stat}=2.5\%$ and $\delta v_{\rm stat}=3.5\%$, respectively. The
last column represents the inferred $H_{\mathrm{0}}$ values for the various
models, under the assumption of flat $\Lambda$CDM. The mock input value is
82.5 km s-1 Mpc-1. We quote the 50th and 16th/84th percentiles of the
distribution.

Model $D_{\mathrm{d},{\rm int}}$ [Mpc] $D_{{\Delta t},{\rm int}}$ [Mpc]
$\lambda_{\mathrm{int}}$ $H_{\mathrm{0}}$ [km s-1 Mpc-1] (flat $\Lambda$CDM)
Lensing - 1749${}^{+897}_{-891}$ 1.00${}^{+0.50}_{-0.50}$ 82.7${}^{+24}_{-22}$
Lensing & Dynamics 791${}^{+34}_{-34}$ 1816${}^{+81}_{-53}$
0.97${}^{+0.03}_{-0.05}$ 82.5${}^{+3.2}_{-3.5}$ FIDUCIAL with $\delta v_{\rm
stat}=2.5\%$ Lensing & Dynamics 813${}^{+36}_{-36}$ 1965${}^{+217}_{-153}$
0.89${}^{+0.08}_{-0.10}$ 79.3${}^{+4.0}_{-4.1}$ FIDUCIAL with $\delta v_{\rm
stat}=3.5\%$

## Appendix C

In this Appendix, we show in more detail that a MST with a large core radius
for the cored surface mass-density profile (that approximates a mass sheet)
effectively yields a vanishing density profile and does not influence the
stellar kinematics of the lens galaxy.

We start from equation (14) for the projected SMD,

$\displaystyle\Sigma_{\rm smd}(\theta)$
$\displaystyle=\bigg{[}\frac{c^{2}}{4\pi G}\;\frac{{D_{{\Delta t},{\rm
int}}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm int}}^{2}}\bigg{]}$ (22)
$\displaystyle\times\bigg{[}(1-\lambda_{\mathrm{int}})+\lambda_{\mathrm{int}}\;\kappa_{\rm
composite}(\theta)\bigg{]}.$

The term $(1-\lambda_{\mathrm{int}})$ represents the mass sheet, and following
B20, we approximate it as a cored dimensionless surface-mass-density
distribution101010Our cored density profile has a different form compared to
the one used in B20.,

$\kappa_{\rm core}(\theta)=\frac{a_{0}\theta_{\rm
c}}{\sqrt{\theta^{2}+\theta_{\rm c}^{2}}}\simeq(1-\lambda_{\mathrm{int}}),$
(23)

where $a_{0}$ is a normalisation parameter that is directly proportional to
the surface mass density of the mass sheet. For $\theta_{\rm c}\gg\theta$, the
cored mass-density profile approximates well a mass sheet of constant
convergence. Assuming spherical symmetry, the mass density corresponding to
the (projected) dimensionless surface-mass density above is

$\rho_{\rm core}(r)=\frac{a_{0}r_{\rm
c}\Sigma_{\mathrm{crit}}}{\pi(r^{2}+r_{\rm c}^{2})},$ (24)

where $\Sigma_{\mathrm{crit}}$ is given by Equation (13). In other words,
$\rho_{\rm core}$ is the deprojection of $\kappa_{\rm core}$.

Substituting equation (23) into equation (22), we obtain

$\displaystyle\Sigma_{\rm smd}(\theta)$
$\displaystyle=\bigg{[}\frac{c^{2}}{4\pi G}\;\frac{{D_{{\Delta t},{\rm
int}}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm
int}}^{2}}\bigg{]}\bigg{[}\kappa_{\rm
core}+\lambda_{\mathrm{int}}\;\kappa_{\rm composite}(\theta)\bigg{]}.$ (25)

In order to predict the $\overline{v_{\mathrm{LOS}}^{2}}$ through the Jeans
equation, we would first deproject the SMD to obtain the mass density

$\displaystyle\rho(r)$ $\displaystyle=\bigg{[}\frac{c^{2}}{4\pi
G}\;\frac{{D_{{\Delta t},{\rm int}}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm
int}}^{2}}\bigg{]}$ (26) $\displaystyle\times\bigg{[}\rho_{\rm
core}(r)+\lambda_{\mathrm{int}}\;{\rm Deproj}\left[\kappa_{\rm
composite}(\theta)\right]\bigg{]},$

where the deprojection for ${\rm Deproj}\left[\kappa_{\rm
composite}(\theta)\right]$ can be done through the MGEs described in Section
3.2. For large $r_{\rm c}$ ($\gtrsim 400\,\rm{kpc}$ for our case), $\rho_{\rm
core}(r)\rightarrow$ 0 and we obtain

$\displaystyle\rho(r)$ $\displaystyle=\frac{c^{2}}{4\pi G}\;\frac{{D_{{\Delta
t},{\rm int}}}\;\lambda_{\mathrm{int}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm
int}}^{2}}\;{\rm Deproj}\left[\kappa_{\rm composite}(\theta)\right].$ (27)

Substituting equation (6) into the above equation for $D_{{\Delta t},{\rm
int}}$, the $\lambda_{\mathrm{int}}$ in the numerator and denominator cancels
and we find that

$\displaystyle\rho(r)$ $\displaystyle=\frac{c^{2}}{4\pi G}\;\frac{{D_{\Delta
t,{\rm composite}}}}{(1+z_{\mathrm{d}})\;D_{\mathrm{d},{\rm int}}^{2}}\;{\rm
Deproj}\left[\kappa_{\rm composite}(\theta)\right].$ (28)

Therefore, when $r_{\rm c}$ is large, this mass density is independent of
$\lambda_{\mathrm{int}}$ and does not affect the stellar kinematics111111We
used $\kappa_{\rm composite}$ for consistency with Sec. 3.1, but the results
presented here are valid irrespective of the parameterisation of the lens mass
profile..