11institutetext:  Department of Astronomy & Astrophysics, University of
Chicago, Chicago, IL 60637, USA; 11email<EMAIL_ADDRESS>22institutetext: Kavli Institute for Cosmological Physics, University of
Chicago, Chicago, IL 60637, USA 33institutetext: Department of Physics and
Astronomy, University of California, Davis, CA 95616, USA 44institutetext:
Department of Physics and Astronomy, University of California, Los Angeles, CA
90095, USA 55institutetext: Sub-Department of Astrophysics, Department of
Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford,
OX1 3RH, UK 66institutetext: Technical University of Munich, TUM School of
Natural Sciences, Department of Physics, James-Franck-Str. 1, Garching, 85748,
Germany 77institutetext: Max Planck Institute for Astrophysics, Karl-
Schwarzschild-Str. 1, Garching, 85748, Germany 88institutetext: Institute of
Astronomy and Astrophysics, Academia Sinica, 11F of ASMAB, No.1, Section 4,
Roosevelt Road, Taipei, 10617, Taiwan 99institutetext: Physics Department,
California Polytechnic State University, San Luis Obispo, CA 93407, USA
1010institutetext: Fermi National Accelerator Laboratory, P.O. Box 500,
Batavia, IL 60510, USA 1111institutetext: STAR Institute, Quartier Agora,
Allée du Six Août, 19c, 4000 Liége, Belgium 1212institutetext: Kavli Institute
for Particle Astrophysics and Cosmology and Department of Physics, Stanford
University, Stanford, CA 94305, USA 1313institutetext: SLAC National
Accelerator Laboratory, Menlo Park, CA, 94025 1414institutetext: Department of
Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA
1515institutetext: Institute of Physics, Laboratory of Astrophysics, Ecole
Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290
Versoix, Switzerland

# TDCOSMO

XII. Improved Hubble constant measurement from lensing time delays using
spatially resolved stellar kinematics of the lens galaxy††thanks: Reduced Keck
Cosmic Web Imager data analyzed in this paper are also available at the CDS
via anonymous ftp to cdsarc.cds.unistra.fr (ftp://130.79.128.5) or via
https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/673/A9. jupyter notebooks and
python scripts used in this analysis are available at
https://github.com/TDCOSMO/RXJ1131_KCWI/.

Anowar J. Shajib  , Pritom Mozumdar NFHP Einstein Fellow11 2 2 Geoff C.-F.
Chen 33 Tommaso Treu 44 Michele Cappellari 44 Shawn Knabel 55 Sherry H. Suyu
44 Vardha N. Bennert 66 7 7 8 8 Joshua A. Frieman 99 Dominique Sluse 11 2 2 10
10 Simon Birrer 1111 Frederic Courbin 1212 13 13 14 14 Christopher D.
Fassnacht 1515 Lizvette Villafaña 33 Peter R. Williams 4444

(Received xxx, xxxx; accepted xxx, xxxx)

Strong-lensing time delays enable measurement of the Hubble constant ($H_{0}$)
independently of other traditional methods. The main limitation to the
precision of time-delay cosmography is mass-sheet degeneracy (MSD). Some of
the previous TDCOSMO analyses broke the MSD by making assumptions about the
mass density profile of the lens galaxy, reaching 2% $H_{0}$ precision from
seven lenses. However, this approach could potentially bias the $H_{0}$
measurement or underestimate the errors. In this work, for the first time, we
break the MSD using spatially resolved kinematics of the lens galaxy in
RXJ1131$-$1231 obtained from the Keck Cosmic Web Imager spectroscopy, in
combination with previously published time delay and lens models derived from
Hubble Space Telescope imaging. This approach allows us to robustly estimate
$H_{0}$, effectively implementing a maximally flexible mass model. Following a
blind analysis, we estimate the angular diameter distance to the lens galaxy
$D_{\rm d}=865_{-81}^{+85}$ Mpc and the time-delay distance $D_{\Delta
t}=2180_{-271}^{+472}$ Mpc, giving $H_{0}=77.1_{-7.1}^{+7.3}$ km s-1 Mpc-1 –
for a flat $\Lambda$ cold dark matter cosmology. The error budget accounts for
all uncertainties, including the MSD inherent to the lens mass profile and the
line-of-sight effects, and those related to the mass–anisotropy degeneracy and
projection effects. Going from single-aperture to spatially resolved stellar
kinematics improves the $H_{0}$ constraint from 13% to 9% for this single
lens, using maximally flexible models. Our new measurement is in excellent
agreement with those obtained in the past by H0LiCOW using simple parametric
models for this single system ($H_{0}$ = $78.3^{+3.4}_{-3.3}$ km s-1 Mpc-1),
and for seven lenses by TDCOSMO with single-aperture kinematics using the same
maximally flexible models used by us ($H_{0}$ = $73.3^{+5.8}_{-5.8}$ km s-1
Mpc-1), corroborating the methodology.

###### Key Words.:

cosmology: distance scale – gravitational lensing: strong – Galaxy: kinematics
and dynamics – Galaxies: elliptical and lenticular, cD – Galaxies: individual:
RXJ1131$-$1231

##  1 Introduction

The Hubble constant, $H_{0}$, the current value of the Universe’s expansion
rate, is a crucial cosmological parameter that also sets the extragalactic
distance scale. Recently, tension has emerged between early- and late-Universe
estimates of $H_{0}$ (e.g., Freedman 2021; Abdalla et al. 2022). The
temperature and polarisation fluctuations in the cosmic microwave background
(CMB) provide an estimate of the Hubble parameter at the last scattering
surface $H(z\approx 1100)$, which can be extrapolated to the current epoch
using the $\Lambda$ cold dark matter ($\Lambda$CDM) cosmology. The CMB
measurements from Planck give $H_{0}=67.4\pm 0.5$ km s-1 Mpc-1 (Planck
Collaboration 2020) and $H_{0}=67.6\pm 1.1$ km s-1 Mpc-1(Aiola et al. 2020).
In the local Universe, $H_{0}$ can be estimated using the cosmic distance
ladder, which uses luminosity distances of type Ia supernovae (SNe Ia) with
their absolute brightness calibrated using different classes of stars. The
Supernova $H_{0}$ for the Equation of State of the dark energy (SH0ES) team
uses Cepheids and parallax distances for this calibration, and they find
$H_{0}=73.04\pm 1.04$ km s-1 Mpc-1 (Riess et al. 2022). This value is in
5$\sigma$ tension with the Planck CMB-based measurements. If this difference
is not due to systematic errors in either of these measurements (e.g.,
Efstathiou 2021), then this tension could point to new physics beyond the
$\Lambda$CDM cosmological model (e.g., Knox & Millea 2020).

To determine whether this “Hubble tension” is due to systematics or new
physics, multiple independent methods to measure $H_{0}$ are needed (e.g.,
Verde et al. 2019; Di Valentino et al. 2021; Freedman 2021). The
Carnegie–Chicago Hubble Project uses the tip of the red giant branch (TRGB) to
calibrate the SNe Ia absolute brightness and measures $H_{0}=69.6\pm 1.9$ km
s-1 Mpc-1 (Freedman et al. 2019, 2020). This TRGB-calibrated measurement is
statistically consistent with both the SH0ES measurement and the CMB-based
measurements. However, several independent local probes strengthen the
“$H_{0}$ tension” by measuring values consistent with the SH0ES value. For
example, the Megamaser Cosmology Project (MCP) estimates $H_{0}=73.9\pm 3.0$
km s-1 Mpc-1(Pesce et al. 2020), the surface brightness fluctuation (SBF)
method measures $H_{0}=73.7\pm 0.7\pm 2.4$ km s-1 Mpc-1 (Blakeslee et al.
2021), and the Tully–Fisher-relation-based method calibrated with Cepheids
measures $H_{0}=75.1\pm 0.2\pm 0.3$ km s-1 Mpc-1 (Kourkchi et al. 2020).

Strong-lensing time delays provide an independent measurement of $H_{0}$
(Refsdal 1964; for an up-to-date review, see Birrer et al. 2022b; Treu et al.
2022, for a historical perspective, see Treu & Marshall 2016). In strong
lensing, a background source appears as multiple images due to the
gravitational deflection of photons by a massive foreground galaxy or galaxy
cluster. The photons that were emitted at the same time from the background
source arrive in different images with a relative time delay. This time delay
carries cosmological information through a combination of angular diameter
distances involved in the lensing system. This combination is referred to as
the “time-delay distance”, which is inversely proportional to $H_{0}$ (Refsdal
1964; Schneider et al. 1992; Suyu et al. 2010). The Time-Delay COSMOgraphy
(TDCOSMO) collaboration has analyzed seven time-delay lenses to measure
$H_{0}$ with 2% error, $H_{0}=74.2\pm 1.6$ km s-1 Mpc-1 assuming a power-law
or composite (i.e., stars and Navarro–Frenk–White (NFW) halo; Navarro et al.
1996, 1997) mass profile for the lensing galaxies (Millon et al. 2020b). The
TDCOSMO collaboration encompasses the COSmological MOnitoring of GRAvItational
Lenses (COSMOGRAIL; Courbin et al. 2005; Millon et al. 2020a), the $H_{0}$
Lenses in COSMOGRAIL’s Wellspring (H0LiCOW; Suyu et al. 2010, 2013; Bonvin et
al. 2017; Birrer et al. 2019; Rusu et al. 2020; Wong et al. 2020), the Strong-
lensing High Angular Resolution Programme (SHARP; Chen et al. 2019), and the
STRong-lensing Insights into the Dark Energy Survey (STRIDES; Treu et al.
2018; Shajib et al. 2020) collaborations.

The simple parametric lens models, e.g., the power law, adopted in the TDCOSMO
analyses are “industry standard” consistent with non-lensing measurements. The
TDCOSMO collaboration has performed various systematic checks on the adopted
lens modeling procedure. These checks find potential systematic biases to be
lower than the acceptable limit ($\sim$1%) from the choice of mass model
parametrization (i.e., power law or composite, Millon et al. 2020b), from
ignoring dark substructures in the lens galaxy’s halo (Gilman et al. 2020),
from ignoring disky or boxy-ness in the baryonic distribution (Van de Vyvere
et al. 2022a), from using different lens modeling software (Shajib et al.
2022a), and from ignoring potential isodensity twists and ellipticity
gradients in the lens galaxy (Van de Vyvere et al. 2022b). However, a
significant source of potential systematics could arise due to the relatively
simple parametrization of the lens mass profile (Kochanek 2020). The well-
known mass-sheet degeneracy (MSD) does not allow one to constrain the mass
profile shape of the deflector galaxy from lens imaging observables alone
(Falco et al. 1985; Schneider & Sluse 2013, 2014). Non-lensing observables,
such as the deflector galaxy’s velocity dispersion or the source’s unlensed
intrinsic brightness, are required to break the mass-sheet degeneracy and
simultaneously constrain $H_{0}$ and the mass profile shape (Treu & Koopmans
2002; Shajib et al. 2018; Yıldırım et al. 2020, 2021; Birrer et al. 2020,
2022a).

The TDCOSMO collaboration has redesigned the experiment to mitigate this
systematic by relaxing the simple parametric assumptions in the mass profile
and constraining the profile shape solely from stellar velocity dispersion
measurements of the lensing galaxies (Birrer et al. 2020). Relaxing the
assumption on the mass profile leads to an increase in the $H_{0}$ uncertainty
from 2 to 8% – which is dominated by the uncertainty of the measured velocity
dispersion – giving $H_{0}=74.5_{-6.1}^{+5.6}$ km s-1 Mpc-1. One approach to
improving the precision is to incorporate prior information on the mass
profile shape from the measured velocity dispersions of a larger sample of
external lenses without measured time delays. Assuming that the Sloan Lens ACS
(SLACS) survey’s lens galaxies are drawn from the same population as the
TDCOSMO lens galaxies and using their velocity dispersions to constrain the
mass profile shape, the uncertainty on $H_{0}$ improves to 5%, giving
$H_{0}=67.4^{+4.1}_{-3.2}$ km s-1 Mpc-1(Birrer et al. 2020). Note that this
estimate is statistically consistent within $1\sigma$ with the larger 8%
$H_{0}$ measurement above. However, the shift could also arise from systematic
differences, e.g., a difference between the parent populations of time-delay
and non-time-delay lenses (Gomer et al. 2022). Such differences could arise,
for example, from evolutionary effects, as the SLACS sample is at lower
redshift than the TDCOSMO lenses (see, e.g., Sonnenfeld et al. 2015, for a
discussion of the evolution of mass density profiles of massive elliptical
galaxies).

Spatially resolved velocity dispersion measurements of lens galaxies for
systems with measured time delays are critical to drastically improving the
$H_{0}$ precision, given the limited sample size of time-delay lenses (Shajib
et al. 2018; Yıldırım et al. 2021). The spatially resolved nature of the
measured velocity dispersion is especially powerful in simultaneously breaking
the MSD and the mass-anisotropy degeneracy (Cappellari 2008; Barnabè et al.
2009, 2012; Collett et al. 2018; Shajib et al. 2018). Spatially resolved
velocity dispersion measurements for $\sim$40 time-delay lens galaxies will
yield an independent $\lesssim$2% $H_{0}$ measurement without any mass profile
assumption (Birrer & Treu 2021). Additional constraints from velocity
dispersion measurements of non-time-delay lens galaxies or magnification
information for standardizable lensed type Ia supernovae can further improve
the uncertainty to $\lesssim 1$% (Birrer & Treu 2021; Birrer et al. 2022a).

In this paper, we measure the spatially resolved velocity dispersion for the
lens galaxy in the strongly lensed quasar system RXJ1131$-$1231using the Keck
Cosmic Web Imager (KCWI) integral field spectrograph on the W. M. Keck
Observatory (Morrissey et al. 2012, 2018). and constrain $H_{0}$ without any
mass profile assumption from this single time-delay lens system. This is the
first application of spatially resolved velocity dispersion from a time-delay
lens to measure $H_{0}$. This lens system was previously used to measure
$H_{0}$ by combining the observed imaging data, single-aperture velocity
dispersion, time delays, and analysis of the line-of-sight environment (Suyu
et al. 2013, 2014). However, these previous studies assumed simple
parametrizations for the mass profile, such as a power law or a combination of
the NFW profile and the stellar profile with constant mass-to-light, which is
the industry standard in modeling of galaxy-scale lenses (Shajib et al.
2022b). Birrer et al. (2016) marginalized over the MSD effect for the system
RXJ1131$-$1231 to constrain $H_{0}$ using a single-aperture velocity
dispersion measurement. Here, we allow the maximal freedom in the MSD by
introducing one free parameter on top of the simply parametrized mass profile
constrained by lens modeling, which is completely degenerate with $H_{0}$.

This paper is organized as follows. In Section 2, we describe the
observational strategy and data reduction. In Section 3, we describe the
procedures to extract the spatially resolved kinematics map from the KCWI
data. In Section 4, we briefly review the lensing and dynamical formalisms and
how we combine the two to mitigate the MSD in our analysis. Then in Section 5,
we describe our dynamical models and present results. We infer the
cosmological parameters from our analysis in the Section 6. We discuss our
results in Section 7 and conclude the paper in Section 8.

We performed the cosmological inference blindly in this paper. The measurement
of velocity dispersion was not blinded. However, we blinded the cosmological
and other model parameters directly related to cosmological parameters in the
dynamical modeling. Before unblinding, this analysis went through an internal
collaboration-wide review and code review. After all the coauthors had agreed
that the necessary systematic checks were satisfactorily performed, we froze
the analysis and unblinded on 5 January 2023. All the sections in this paper
except for the final discussion in Section 7 and summary in Section 8 were
written before unblinding. After unblinding, we only made minor edits for
clarity and grammatical corrections in the previous sections and added the
unblinded numbers where relevant in the abstract, main text, and plots.

##  2 Observations and data reduction

In this section, we provide a brief description of the lens system
RXJ1131$-$1231 (Section 2.1), the spectroscopic observation with KCWI (Section
2.2), and the data reduction procedure (Section 2.3).

Figure 1:  HST/ACS image of RXJ1131$-$1231 in the F814W band. The four quasar
images are labeled with A, B, C, and D. The central deflector is marked with
G, of which we are measuring the spatially resolved velocity dispersion. An
arrow points to the nearby satellite S, which we mask out in the velocity
dispersion measurement.. The North and East directions and 1s̈cale are also
illustrated.

###  2.1 Description of lens system

The quadruply imaged quasar lens system RXJ1131$-$1231 was discovered by Sluse
et al. (2003). The deflector in this system is an elliptical galaxy with
redshift $z_{\rm d}=0.295$, and the source redshift is $z_{\rm s}=0.657$
(Sluse et al. 2003). Due to its low redshifts, the system is relatively bright
and large in angular size. The Einstein ring in this system contains intricate
features, providing a wealth of information to constrain the lens mass model
(see Figure 1). Due to its early discovery and information-rich features, this
system is one of the most studied lensed quasar systems. The time delays for
this system were measured by Tewes et al. (2013). Suyu et al. (2013, 2014)
performed cosmographic analyses of this system. These authors combined simply
parametrized lens models based on the high-resolution imaging from the HST’s
Advanced Camera for Surveys (ACS) instrument (HST-GO 9744; PI: Kochanek), the
measured time delays, single-aperture velocity dispersion, and external
convergence estimate to infer $H_{0}=80.0_{-4.7}^{+4.5}$ km s-1 Mpc-1.
However, such simply parameterized lens models implicitly break the MSD.
Birrer et al. (2016) performed an independent mass modeling of this system
while marginalizing the MSD with a prior on the source size. These authors
found the prior choice on the anisotropy in the dynamical modeling to be the
dominant systematic in inferring $H_{0}$.

###  2.2 KCWI Spectroscopy

We obtained integral field unit (IFU) spectroscopy of RXJ1131$-$1231 on 16 May
and 7 June 2021 with the KCWI instrument on the Keck Observatory (Morrissey et
al. 2012, 2018). We chose KCWI with the small IFU slicer and the low-
resolution blue grating (BL) with a field-of-view (FoV) of
$8\aas@@fstack{\prime\prime}4\times 20\aas@@fstack{\prime\prime}4$. The
spectral resolution is $R\approx 3600$, corresponding to an instrumental
dispersion $\sigma_{\rm inst}\sim 35$ km s-1. The reciprocal dispersion is
$0.5\ \AA$ per pixel. The observed wavelength range 3600–5600 Å covers the Ca
H&K lines with wavelengths $\lambda\lambda 3933$, $3968$ Å at the redshift of
the lens galaxy ($z_{\rm d}=0.295$). We primarily use these lines to determine
the stellar velocity dispersion. The redshifted 4304 Å G-band is beyond the
observed range, so it is not accessible with the KCWI for the RXJ1131$-$1231
system.

We aligned the FoV’s longer side with the North direction (i.e.,
$\textrm{PA}=0\degr$) and dithered the individual exposures by $9\arcsec$
along the North-South direction. As the extent of the RXJ1131$-$1231 system is
smaller than the FoV, each exposure contained the entire lens system within
the FoV. In different exposures, the lens system occupied the upper or lower
portion of the FoV. Thus, the sky in an exposure with the system occupying the
upper portion can be subtracted using another exposure with the system
occupying the lower portion, and vice versa. We obtained six exposures with a
total integration time of 10,560 s on 16 May and three with a total
integration time of 5,400 s on 7 June. Therefore, the total exposure time is
$t_{\textrm{exp}}=15,960$ s. The airmass ranged from 1.2 to 1.48 over the
integrating period.

###  2.3 Data Reduction

We use the official Python-based data reduction pipeline111developed by Luca
Rizzi, Don Neill, Max Brodheim; https://kcwi-drp.readthedocs.io/ (DRP) to
reduce our data. The DRP converts the 2D raw data captured on the detector
into a 3D datacube. It performs geometry correction, differential atmospheric
refraction correction, and wavelength calibration and produces a final
standard-star-calibrated 3D datacube for each exposure. The calibration with
the standard star corrects for instrumental response and scales the data to
flux units (Morrissey et al. 2018). We use the final output file with the
suffix “_icubes” for further analysis.

We stack the dithered datacubes through drizzling (Fruchter & Hook 2002).
Since the exposures are obtained on different dates, the world coordinate
system information is not accurate enough to determine the relative positions
of the dithered exposures. We follow Chen et al. (2021b) to determine the
relative positions by simultaneously fitting the point spread function (PSF)
to the four quasar image positions. To perform the drizzling on the datacubes,
we repurpose the drizzling routine of the DRP for OSIRIS, another IFU
spectrograph on the Keck Observatory222https://github.com/Keck-
DataReductionPipelines/OsirisDRP. For the drizzling process, we set
$\texttt{pixfrac}=0.7$ as recommended to reduce correlated uncertainties
between the drizzled pixels (Avila et al. 2015). We calculate the drizzled
weight image and ensure that the ratio of RMS/median $<0.2$ in the region of
interest so that the trade-off is balanced between improving the image
resolution and increasing the background noise (Gonzaga et al. 2012). The
rectangular pixel size $0\aas@@fstack{\prime\prime}1457\arcsec\times
0\aas@@fstack{\prime\prime}3395$ of the KCWI is kept the same in the drizzled
output. We transform the datacube to have square pixels of size
$0\aas@@fstack{\prime\prime}1457\times 0\aas@@fstack{\prime\prime}1457$
through resampling while conserving the total flux. We converted the pixels
into square sizes for the convenience of Voronoi binning the spectra using the
software vorbin as described in Section 3.2.

We directly estimate the PSF from the observed data. We produce a 2D image
from the datacube by summing along the wavelength axis (see Figure 2). We
create a model for this KCWI image using a high-resolution template from the
HST imaging (Figure 1) that has a pixel size $0\aas@@fstack{\prime\prime}05$
and PSF full width at half maximum (FWHM) $0\aas@@fstack{\prime\prime}10$. In
the model, the template is convolved with a Gaussian PSF with a free FWHM
parameter, and the positioning of the template on the KCWI image grid is
fitted with two additional free parameters. By optimizing the model, we
estimate that the PSF FWHM is 0$\aas@@fstack{\prime\prime}$96.

Figure 2:  Left: 2D representation (median-collapsed) of the 3D KCWI datacube
for RXJ1131$-$1231\. The yellow contour traces the region with
$1\aas@@fstack{\prime\prime}5$ radial extent from the center selected for
stellar kinematic measurement. A circular region with
$0\aas@@fstack{\prime\prime}5$ radius around image D and the spaxel containing
the satellite S are excluded from this selected region. All the individual
spaxels within this region have continuum $S/N>1.4\ \AA^{-1}$ for the lens
galaxy’s light within 3985–4085 $\AA$ (the purple shaded range in the right
panel). Right: The spectra (grey) from an example pixel (grey box in the left
panel) and the estimate of the signal from the lens galaxy’s spectra (orange)
after removing the contribution from the quasar light (blue). The full model
of the spectra is presented with the red line, and the model’s residual is
plotted in emerald color. The vertical purple shaded region marks where we
compute the continuum $S/N$. Figure 3:  Distribution of the stellar spectral
types in the XSL according to the Simbad database. Unspecified stars are
grouped in the ‘$\sim$’ class. The dark grey color represents the full library
of 628 stars. Set 1 (orange) refers to the 39 stars selected by pPXF out of
the full library to construct an optimal template 1. Set 2 (blue) refers to 32
stars selected from a random half of the full library and set 3 (emerald)
refers to 33 stars selected from the other half. Sets 2 and 3 have 15 and 17
stars, respectively, in common with Set 1. The alternating light grey and
white vertical regions divide the spectral classes for easier visualization.

##  3 Kinematics maps

This section describes our procedure to obtain the final kinematics map. We
use the pPXF package333https://pypi.org/project/ppxf/ to fit the spectra with
a library of stellar templates and extract the velocity dispersion (Cappellari
2017, 2022). In Section 3.1, we describe the stellar templates used for the
analysis. In Section 3.2, we present the measurement of the spatially-resolved
kinematics map of the lens galaxy. In Section 3.3, we test the systematics of
the velocity dispersion measurement.

###  3.1 Library of Stellar Templates

The popularly used template libraries Medium-resolution Isaac Newton Telescope
library of empirical spectra (MILES; Sánchez-Blázquez et al. 2006) and INDO-US
templates (Valdes et al. 2004) are both too low resolution to fit our
datasets. The KCWI’s instrumental resolution of $R\approx 3600$ leads to
$\sigma_{\rm inst}\sim 35$ km s-1 for a Gaussian line spread function
(LSF)444We quantitatively verified that the shape of the instrumental LSF is
Gaussian (cf. Figure 28 of Morrissey et al. 2018). Thus, the treatment of the
instrumental LSF in pPXF is self-consistent and avoids any systematic bias due
to inconsistent definitions of the LSF’s FWHM (Robertson 2013).. MILES has a
resolution of $\sigma_{\rm template}\sim 64$ km s-1 (i.e., $R\sim 2000$), and
the INDO-US templates have an approximately constant-wavelength resolution of
1.2 Å, which corresponds to $\sigma_{\rm template}=39$ km s-1 over the Ca H&K
wavelength range. Therefore, we choose the X-shooter Spectral Library (XSL),
which contains 628 stars covering three segments, including UVB, Vis, and NIR
bands (Gonneau et al. 2020). As our data cover the rest-frame blue/UV range,
we only use the UVB segment to fit the data, where its resolution is $R\sim
9700$ and $\sigma_{\rm template}\sim 13$ km s-1.

Figure 4:  Left: Voronoi binning of the selected spaxels within
$1\aas@@fstack{\prime\prime}5$ from the galaxy center that avoid lensed arcs,
quasar images, and the satellite galaxy S. The different colors illustrate the
regions for each Voronoi bin in a cartographic manner for easier
visualization, with the bin number specified within each bin. We perform the
binning with a target $S/N\approx 23\ \AA^{-1}$ for each bin, which results in
41 bins in total. Right: Resultant $S/N$ for each Voronoi bin (red points).
Figure 5: pPXF fitting to the spectra from four examples of Voronoi bins. The
bin number and the measured velocity dispersion for the corresponding bin are
specified in each panel. The grey line presents the full spectra, the red line
traces the best-fit model, and the blue line shows the quasar component in the
best-fit model.

###  3.2 Measuring the velocity dispersion

We choose a cutout centred on the lens system with
$6\aas@@fstack{\prime\prime}235\times 6\aas@@fstack{\prime\prime}235$ (43
pixels $\times$ 43 pixels) to initiate the analysis (see Figure 2). We
estimate the lens galaxy light’s signal-to-noise ratio ($S/N$) in each spatial
pixel (hereafter, spaxel) within this initial cutout. We then select a region
with sufficient $S/N$ from the lens galaxy and relatively low quasar
contamination for measuring the velocity dispersion (the yellow contour in
Figure 2’s left panel). We perform Voronoi binning within this selected region
to preserve the maximal spatial resolution and reduce the bias in the
lower-$S/N$ region (Cappellari & Copin 2003). We elaborate on these steps
below.

To estimate the lens galaxy’s $S/N$ in each spaxel, we first simultaneously
fit the quasar and the lens galaxy in each spaxel to calculate the signal from
each of them. We perform this fitting within the wavelength range 3400–4300 Å.
As the four quasar images surround the lens galaxy, each spaxel receives a
different contribution from the quasar light. We take spectra at the central
spaxel of image A as the quasar template, ignoring the lens galaxy’s small
contribution. Later in Section 3.3, we also choose the quasar template from
images B and C to account for the associated systematic uncertainty, i.e., the
potential impact of chromatic microlensing that may change the contrast
between the line and the continuum (e.g., Sluse et al. 2007).

We determine a single optimal template spectrum for the lens galaxy template.
For this purpose, we binned the spectra from spaxels within a circular region
of radius $0\aas@@fstack{\prime\prime}5$ centered on the lens galaxy and fit
it with pPXF using the 628 stellar templates from the XSL and the quasar
template. We also include a Legendre polynomial of degree 3 as a component in
the fitting to account for any residual gradient in the continuum. pPXF
chooses 39 of the stellar templates and builds the optimal template by taking
a weighted linear combination of them. See Figure 3 for the weighted
distribution of spectral types of the full template library and that of the 39
stars selected by pPXF. Among those stars in the XSL with stellar classes
specified by the Simbad database (Wenger et al. 2000), G-type stars are
selected with the highest total weight, consistent with the fact that massive
elliptical galaxy spectra are dominated by G and K-type stars. In the pPXF
fitting procedure, the stellar templates are broadened, corresponding to a
freely varying velocity dispersion, but the velocity dispersion does not
broaden the quasar template.

Once the optimal galaxy template is constructed, we use this template and the
quasar template to fit the spectrum of each spaxel individually. We use this
optimal template to fit the galaxy spectra in individual spaxels instead of
the full template library to avoid large spurious fluctuations in the measured
velocity dispersion from spaxel to spaxel. We show the decomposition of the
spectra from one example spaxel into different components after fitting with
pPXF in Figure 2. We calculate the signal of the lens galaxy’s spectrum in
each spaxel by subtracting the modeled quasar component from the observed
spectra. The noise is estimated by adding in quadrature the Poisson noise of
the total signal and the background noise estimated from an empty patch of the
sky. The noise values are multiplied by $\sqrt{2}$ to account for the fact
that the square pixels are created from the rectangular pixels about double
the size through resampling. We estimate the $S/N$ using the restframe
wavelength range 3985–4085 Å, slightly above the Ca H&K absorption lines in
wavelength (see the purple shaded region in Figure 2).

To perform Voronoi binning before the velocity dispersion measurement, we
select the spaxels within a radius of $1\aas@@fstack{\prime\prime}5$555For
reference, $1\aas@@fstack{\prime\prime}5$ corresponds to 6.6 kpc at $z_{\rm
d}=0.295$ for a fiducial flat $\Lambda$CDM cosmology with $H_{0}=70$ km s-1
Mpc-1 and $\Omega_{\rm m}=0.3$. from the lens galaxy center that avoid the
brightest spaxels containing images A, B, and C and the lensed arcs. We also
exclude a circular region around image D with radius
$0\aas@@fstack{\prime\prime}5$. To avoid any potential bias due to
contamination from the satellite galaxy S, we exclude the spaxel at its
position ($\Delta{\rm RA}=0\aas@@fstack{\prime\prime}09$, $\Delta{\rm
Dec}=0\aas@@fstack{\prime\prime}54$ from the galaxy center, Suyu et al. 2013).
We also exclude pixels with $S/N<1\ \AA^{-1}$. In the end, the spaxels within
the selected region have $S/N>1.4\ \AA^{-1}$ (see Figure 2 for the selected
region). We perform Voronoi binning using
vorbin666https://pypi.org/project/vorbin/ given the estimated $S/N$ values for
each spaxel. In Figure 4, we show the 41 Voronoi bins obtained by setting the
target $S/N\approx 23\ \AA^{-1}$ for each bin. This target $S/N$ was chosen so
that the resultant $S/N\gtrsim 20\ \AA^{-1}$ for each bin, which is standard
practice (Figure 4, only bin 16 has $S/N\approx 18\ \AA^{-1}$).

For each Voronoi bin, we measure the velocity dispersion by fitting the binned
spectra using pPXF using the optimal galaxy template described above, the
quasar template, and the additive Legendre polynomial to model any slight
gradient in the population. A few examples of pPXF fit of the binned spectra
are shown in Figure 5.

###  3.3 Estimation of systematic uncertainty

To estimate the systematic uncertainties in the velocity dispersion
measurement, we consider a range of plausible choices in the extraction
procedure: the degrees of the additive Legendre polynomial used to correct the
template continuum shape between 2 to 4; the quasar template obtained from
images A, B, and C; the fitted wavelength range chosen from 3300–4200 Å,
3350–4250 Å, and 3400–4300 Å; and three sets of template spectra used in the
fitting. The first set of template spectra contains the complete XSL of 628
stars. The second set contains half of the entire sample that is randomly
selected, and the third set contains the other half. The numbers of stars
selected by pPXF in the three sets are 39, 32, and 33, respectively. Sets 2
and 3 have 15 and 17 stars, respectively, in common with Set 1. Figure 3 shows
the distribution of spectral types in all three sets and the entire library.
We do not take the quasar template from image D as it is much fainter than the
other images, and thus the galaxy contribution in the brightest spaxel on
image D is non-negligible. Taking a combination of all of these choices yields
81 different setups. We illustrate the shift in the extracted velocity
dispersion maps for one change of setting at a time in Figures 6 and 7.

We estimate the variance-covariance matrix of the binned velocity dispersions
from these 81 setups. To do this, we generate 1,000 random realizations of the
measured velocity dispersion map for each of the 81 setups using the
corresponding statistical uncertainty. We create the variance-covariance
matrix from the 81,000 realizations combined from all the setups. In this way,
the diagonal terms of the variance-covariance matrix encode the total variance
from systematic and statistical uncertainties, and the off-diagonal terms
encode the systematic covariances. For example, if all 81 setups
hypothetically provided the same velocity dispersion map and uncertainty, then
the off-diagonal terms would be zero, and the diagonal terms would reflect
only the statistical uncertainties. We show the systematic variance-covariance
relative to the statistical variance in Figure 8. The systematic variance is
subdominant relative to the statistical variance (with a median of 0.47 of the
ratio between systematic and statistical covariances along the diagonal)
except for bins 29 and 31. These two bins are closest to quasar images A and
C, and thus largely susceptible to the choice of quasar template (see Figures
6 and 7.)

We show the velocity dispersion and mean velocity maps averaged over the 81
setups in Figure 9. We estimate a systematic velocity of 182 km s-1 using the
pafit777https://pypi.org/project/pafit/ software program (Krajnović et al.
2006) and subtract it from the mean velocity map. The systematic velocity is
the result of a slight deviation in the true redshift from the fiducial value.
The mean velocity map does not show any significant evidence of ordered
rotation above the systematic and statistical noise levels. Thus it is
consistent with the lens galaxy being a slow rotator. We use this systematic-
averaged velocity dispersion map and the variance-covariance matrix estimated
above when computing the likelihood function for dynamical modeling in Section
5.

To test the impact of our choice for the Voronoi binning scheme, we adopt an
alternative target $S/N\approx 28\ \AA^{-1}$ for each bin, which results in 27
bins. We similarly produce another set of 81 model setups in this binning
scheme and produce the variance-covariance matrix for these binned velocity
dispersions. We test the systematic impact of this different binning scheme on
the cosmological measurement later in Section 5.2. We show the difference in
the extracted kinematics between the two binning schemes in Figure 10.

Figure 6:  Absolute difference in km s-1 between the extracted velocity
dispersion from two setups that differ by one setting. The baseline setup has
the range: 3400–4300 Å, polynomial degree: 3, stellar template set 1, and
quasar template from image A. The different setting for each case is specified
at the top of each panel.  Figure 7:  Same as Figure 6, but the difference is
normalized by the statistical uncertainty of the baseline setup.  Figure 8:
Illustration of the systematic covariance relative to the statistical
covariance. $\Sigma$ is the variance-covariance matrix of the Voronoi-binned
velocity dispersions (with target $S/N\approx 23\ \AA^{-1}$ for each bin),
$\sigma_{\rm stat,\it x}$ is the statistical uncertainty in bin number $x$
from our fiducial setup, and ${\rm diag}(\sigma_{\rm stat})$ is a diagonal
matrix. Note that we assume no covariance in the statistical uncertainty from
each setup for kinematic measurement. Thus the off-diagonal terms in the
variance-covariance matrix purely represent the systematic covariance. Most
diagonal terms are $<1$ (with a median of 0.47), showing that the systematic
variances are subdominant to the statistical variances except for bins 29 and
31. These bins are close to images A and C. Thus, they are largely susceptible
to the choice of the quasar template, as seen in Figures 6 and 7.  Figure 9:
Maps of extracted velocity dispersion (top row) and mean velocity (bottom row)
in Voronoi bins along with the corresponding uncertainties (right column). The
Voronoi binning was tuned to achieve $S/N\approx 23\ \AA^{-1}$ for each bin.
The illustrated maps (left column) correspond to the average values after
combining 81 model setups, and the uncertainty maps correspond to the square
root of the diagonal of the variance-covariance matrices. A systematic
velocity of 182 km s-1 was subtracted from the mean velocity map. Figure 10:
Absolute (left) and uncertainty-normalized (right) difference in the extracted
velocity dispersion between two Voronoi binning schemes. The two binning
schemes are obtained by setting the target $S/N$ to 23 $\AA^{-1}$ and 28
$\AA^{-1}$ for each bin. We take the case with target $S/N\approx 23\
\AA^{-1}$ for each bin as the baseline in our analysis.

##  4 Overview of lens and dynamical modeling

This section reviews the theoretical formalism of lens and dynamical modeling.

###  4.1 Lensing observables and modeling

We briefly review the strong lensing formalism in the context of time-delay
cosmography in Section 4.1.1, describe the mass-sheet transform (MST) in
Section 4.1.2, and explain the internal and external components of the MST in
Section 4.1.3.

####  4.1.1 Strong lensing formalism

In the thin lens approximation applicable in this case, lensing observables
are described using the surface mass density $\Sigma(R)$ projected from the 3D
mass density distribution $\rho(r)$ in the lens galaxy. Formally, the lensing
observables depend on the dimensionless convergence defined as

$\kappa(\bm{\theta})\equiv\frac{\Sigma(\bm{\theta})}{\Sigma_{\rm cr}},$ (1)

which is the surface mass density normalized by the critical density

$\Sigma_{\rm cr}\equiv\frac{c^{2}D_{\rm s}}{4\uppi GD_{\rm d}D_{\rm ds}}.$ (2)

Here, $c$ is the speed of light, $G$ is the gravitational constant, $D_{\rm
s}$ is the angular diameter distance between the observer and the source,
$D_{\rm d}$ is the angular diameter distance between the observer and the lens
galaxy, and $D_{\rm ds}$ is the angular diameter distance between the lens
galaxy and the source. The on-sky deflection angle $\bm{\alpha}(\bm{\theta})$
relates to the convergence as

$\kappa(\bm{\theta})=\frac{1}{2}\nabla\cdot\bm{\alpha}(\bm{\theta}).$ (3)

The time delay between two quasar images labeled A and B is given by

$\begin{split}\Delta t_{\rm AB}&=\frac{D_{\Delta
t}}{c}\left[\frac{(\bm{\theta}_{\rm
A}-\bm{\varsigma})^{2}}{2}-\frac{(\bm{\theta}_{\rm
B}-\bm{\varsigma})^{2}}{2}-\psi(\bm{\theta}_{\rm A})+\psi(\bm{\theta}_{\rm
B})\right],\end{split}$ (4)

where $\bm{\theta}_{\rm A}$ is the angular position of image A,
$\bm{\varsigma}$ is the source’s angular position, $\psi(\bm{\theta})$ is the
lensing potential, and the time-delay distance $D_{\Delta t}$ is defined as

$D_{\Delta t}\equiv(1+z_{\rm d})\frac{D_{\rm d}D_{\rm s}}{D_{\rm ds}}.$ (5)

####  4.1.2 Description of the MST

The MST is a mathematical transform of the convergence profile that leaves
invariant all the imaging observables, such as the image positions and the
flux ratios (Falco et al. 1985; Schneider & Sluse 2014). This transform scales
the convergence and the unknown source position as

$\displaystyle\kappa\to\kappa^{\prime}=\lambda_{\rm MST}\kappa+(1-\lambda_{\rm
MST}),$ (6) $\displaystyle\varsigma\to\varsigma^{\prime}=\lambda_{\rm
MST}\varsigma.$

where $\lambda_{\rm MST}$ is the transformation parameter. The predicted time
delay $\Delta t$ scales under the transform as

$\displaystyle\Delta t\to\Delta t^{\prime}=\lambda_{\rm MST}\Delta t.$ (7)

Then, the inferred time-delay distance $D_{\Delta t}$ and the Hubble constant
$H_{0}$ based on the observed time delays will change as

$\displaystyle D_{\Delta t}^{\prime}=\frac{D_{\Delta t}}{\lambda_{\rm MST}},$
(8) $\displaystyle H_{0}^{\prime}=\lambda_{\rm MST}H_{0}.$

However, the MST changes the predicted velocity dispersion, thus measuring it
breaks the MSD. Notably, the MST also rescales the lensing magnifications.
Thus, standardizable candles can also be used to break the MSD (Bertin &
Lombardi 2006; Birrer et al. 2022a) provided that microlensing and
millilensing can be mitigated (e.g., Yahalomi et al. 2017; More et al. 2017;
Foxley-Marrable et al. 2018).

####  4.1.3 Internal and external MST

We can express the “true” (i.e., physically present) lensing mass distribution
as

$\kappa_{\rm true}=\kappa_{\rm gal}+\kappa_{\rm ext},$ (9)

where $\kappa_{\rm gal}$ is the mass distribution of the central lens galaxy
(or galaxies) that is (are) considered in the lens modeling, and $\kappa_{\rm
ext}$ is called the external convergence, which approximates the projected
mass distribution of line-of-sight structures as a mass sheet. Since
$\lim_{\theta\to\infty}\kappa_{\rm gal}=0$ has to be satisfied, we find that
$\lim_{\theta\to\infty}\kappa_{\rm true}=\kappa_{\rm ext}$, hence the
interpretation of $\kappa_{\rm ext}$ as the lensing mass far from (or,
“external” to) the central deflector(s).

All the lensing observables including imaging observables result from
$\kappa_{\rm true}$. However, since only the central galaxies are usually
considered in lens modeling with imaging observables, the lens model provides
$\kappa_{\rm model}^{\prime}$ with $\lim_{\theta\to\infty}\kappa_{\rm
model}^{\prime}=0$. This $\kappa_{\rm model}^{\prime}$ is an MST of
$\kappa_{\rm true}$ for $\lambda_{\rm MST}=1/(1-\kappa_{\rm ext})$ as

$\kappa_{\rm model}^{\prime}=\frac{\kappa_{\rm gal}+\kappa_{\rm
ext}}{1-\kappa_{\rm ext}}+1-\frac{1}{1+\kappa_{\rm ext}}=\frac{\kappa_{\rm
gal}}{1-\kappa_{\rm ext}}.$ (10)

Lens mass models are usually described with simply parametrized models, such
as the power law or a combination of the NFW profile and the observed stellar
distribution. In that case, the assumption of a simple parametric form
implicitly breaks the MSD. Therefore, the simply parametrized model
$\kappa_{\rm model}$ can be expressed as another approximate MST of the
$\kappa_{\rm model}^{\prime}$ as

$\kappa_{\rm model}^{\prime}\approx\lambda_{\rm int}\kappa_{\rm
model}+(1-\lambda_{\rm int})\kappa_{\rm s}(\theta),$ (11)

where $\lambda_{\rm int}$ is called the internal MST parameter, and
$\kappa_{\rm s}$ is a “variable” mass sheet with
$\lim_{\theta\to\infty}\kappa_{\rm s}(\theta)=0$ to ensure that both
$\lim_{\theta\to\infty}\kappa_{\rm model}^{\prime}=0$ and
$\lim_{\theta\to\infty}\kappa_{\rm model}=0$ are satisfied. However, for
Equation (11) to be an approximate MST, the variable mass-sheet needs to
satisfy $\kappa_{\rm s}(\theta)\simeq 1$ within the central region that
lensing observables are sensitive to ($\theta\lesssim 2\theta_{\rm E}$,
Schneider & Sluse 2013). This can be achieved with the formulation (Blum et
al. 2020)

$\kappa_{\rm s}(\theta)=\frac{\theta_{\rm s}^{2}}{\theta^{2}+\theta_{\rm
s}^{2}},$ (12)

where $\theta_{\rm s}\gg\theta_{\rm E}$ is a scale radius where the variable
mass-sheet smoothly transitions from $1-\lambda_{\rm int}$ to 0. This
approximate MST converges to the pure MST in the limit $\theta_{\rm
s}\to\infty$. Thus, the actual mass distribution of the central deflector(s)
relates to the modeled mass distribution as

$\kappa_{\rm gal}\approx(1-\kappa_{\rm ext})\left[\lambda_{\rm int}\kappa_{\rm
model}+(1-\lambda_{\rm int})\kappa_{\rm s}(\theta)\right].$ (13)

The external convergence $\kappa_{\rm ext}$ can be estimated by using relative
number counts of line-of-sight galaxies near the central deflector(s) (e.g.,
Suyu et al. 2010; Greene et al. 2013; Rusu et al. 2017; Buckley-Geer et al.
2020), or by using weak lensing of distant galaxies by the line-of-sight mass
distribution (e.g., Tihhonova et al. 2018). The measured velocity dispersion
then constrains the internal MST parameter $\lambda_{\rm int}$ (Birrer et al.
2020; Yıldırım et al. 2021).

###  4.2 Dynamical modeling

In this section, we describe the Jeans anisotropic multi-Gaussian-expansion
(JAM) framework to model our dynamical observable, which is the spatially
resolved stellar velocity dispersion measured in Section 3. The orbital
motions of the stars, i.e., the distribution function $f(\bm{x},\bm{v})$ of
position $\bm{x}$ and velocity $\bm{v}$, in the galactic potential $\Phi$ is
described by the steady-state collisionless Boltzmann equation (Binney &
Tremaine 1987, Eq. 4-13b)

$\sum_{i=1}^{3}\left(v_{i}\frac{\partial f}{\partial
x_{i}}-\frac{\partial\Phi}{\partial x_{i}}\frac{\partial f}{\partial
v_{i}}\right)=0.$ (14)

We assume an axisymmetric case (i.e., $\partial\Phi/\partial\phi=\partial
f/\partial\phi=0$ with $\phi$ being the polar angle in the spherical
coordinate system), a spherically aligned velocity ellipsoid, and the
anisotropy for each Gaussian component in the multi-Gaussian expansion (MGE;
Emsellem et al. 1994; Cappellari 2002) to be spatially constant. Slow rotators
such as the deflector galaxy in RXJ1131$-$1231 are in general expected to be
weakly triaxial or oblate but never flat and instead quite close to spherical
in their central parts (e.g., Cappellari 2016). For this reason, we expect the
spherical alignment of the velocity ellipsoid of jamsph (Cappellari 2020) to
provide a better approximation to the galaxy dynamics than the cylindrical
alignment jamcyl solution (Cappellari 2008). Then, the above equation gives
two Jeans equations in spherical coordinates (Jeans 1922; Bacon et al. 1983;
de Zeeuw et al. 1996; Cappellari 2020)

$\displaystyle\frac{\partial\left(\zeta\langle{v_{r}^{2}}\rangle\right)}{\partial
r}+\frac{(1+\beta)\zeta\langle{v_{r}^{2}}\rangle-\zeta\langle{v_{\phi}^{2}}\rangle}{r}$
$\displaystyle=-\zeta\frac{\partial\Phi}{\partial r},$ (15)
$\displaystyle(1-\beta)\frac{\partial\left(\zeta\langle{v_{r}^{2}}\rangle\right)}{\partial\theta}+\frac{(1-\beta)\zeta\langle{v_{r}^{2}}\rangle-\zeta\langle{v_{\phi}^{2}}\rangle}{\tan\theta}$
$\displaystyle=-\zeta\frac{\partial\Phi}{\partial\theta},$

where the following notations are used

$\displaystyle\zeta\langle{v_{p}v_{q}}\rangle$ $\displaystyle\equiv\int
v_{p}v_{q}f{\rm d}^{3}\bm{v},$ (16) $\displaystyle\beta$ $\displaystyle\equiv
1-\frac{\langle{v_{\theta}^{2}}\rangle}{\langle{v_{r}^{2}}\rangle}.$

Here, $\beta$ is the anisotropy parameter, and the velocity dispersion
ellipsoid is assumed to be spherically aligned, giving
$\langle{v_{r}v_{\theta}}\rangle=0$.

The line-of-sight second moment $\langle{v_{\rm los}^{2}}\rangle$ is the
integral given by

$S\langle{v_{\rm los}^{2}}\rangle(x,y)=\int_{-\infty}^{\infty}{\rm d}z\
\zeta\langle{v_{z}^{2}}\rangle.$ (17)

where $S(x,y)$ is the surface density of the dynamical tracer. Given that
there is no evidence of significant ordered rotation and the only
significantly nonzero velocities are likely due to systematic errors (see
Figure 9), we assume $\langle{v_{\rm los}}\rangle=0$ and define
$\langle{v_{\rm los}^{2}}\rangle=\sigma_{\rm los}^{2}$. The observed line-of-
sight velocity dispersion is given a luminosity-weighted integral as

$\left[\sigma^{2}_{\rm los}\right]_{\rm obs}=\left[\langle{v_{\rm
los}^{2}}\rangle\right]_{\rm obs}=\frac{\int_{\rm ap}{\rm d}x{\rm d}y\
I\langle{v_{\rm los}^{2}}\rangle\otimes{\rm PSF}}{\int_{\rm ap}{\rm d}x{\rm
d}y\ I\otimes{\rm PSF}},$ (18)

where the symbol “$\otimes\ {\rm PSF}$” denotes a convolution with the PSF. In
the equation above, we have chosen the surface brightness profile $I(x,y)$ as
a substitute for the surface density $S(x,y)$ of the dynamical tracer since
the constant factor between surface brightness and surface number density
cancels out in this expression.

We use the dynamical modeling software jampy888https://pypi.org/project/jampy/
to compute the observed velocity dispersion by solving the Jeans equation from
Equation (15) for a given 3D potential $\Phi(r)$ and anisotropy profile
$\beta(r)$. Specifically, we use the jam_axi_proj() routine with the keyword
align=‘sph’. See Cappellari (2008, 2020) for a detailed formalism in computing
Equation (18) by jampy.

###  4.3 Cosmological inference from combining dynamical and lensing
observables

We parametrize the 3D potential $\Phi(r)$ using the lens model parameters
$\xi_{\rm mass}$ and the internal MST parameter $\lambda_{\rm int}$ to
conveniently use the lens model posterior from Suyu et al. (2013) as a mass
model prior in the dynamical modeling. Thus from Equation (13), the surface
mass density for our dynamical model is given by

$\Sigma(\bm{\theta})=\Sigma_{\rm cr}(1-\kappa_{\rm ext})\left[\lambda_{\rm
int}\ \kappa_{\rm model}(\bm{\theta})+(1-\lambda_{\rm int})\kappa_{\rm
s}(\theta)\right].$ (19)

We include $D_{\Delta t}$ and $D_{\rm d}$ as free parameters in our model,
which give the critical density $\Sigma_{\rm cr}$ as

$\Sigma_{\rm cr}=\frac{c^{2}}{4\uppi G}\frac{D_{\Delta t}}{(1+z_{\rm d})D_{\rm
d}^{2}}=\frac{c^{2}}{4\uppi G}\frac{D_{\Delta t}^{\rm model}}{(1+z_{\rm
d})(1-\kappa_{\rm ext})\lambda_{\rm int}D_{\rm d}^{2}},$ (20)

where $D_{\Delta t}^{\rm model}$ is the time-delay distance predicted by the
lens mass model $\kappa_{\rm model}(\bm{\theta})$ for the time delays observed
by Tewes et al. (2013).

We approximate the surface mass density $\Sigma(\bm{\theta})$ with an MGE
(Emsellem et al. 1994; Cappellari 2002; Shajib 2019) using the software
program mgefit999https://pypi.org/project/mgefit/. jampy deprojects the MGE
components into an oblate or prolate spheroid with an inclination angle $i$
(Cappellari 2002). The deprojected 3D mass density provides the 3D potential
$\Phi$ for the kinematic computation. We also take the MGE of the surface
brightness $I(x,y)$ for deprojection to 3D with the inclination angle $i$ for
the kinematic computation by jampy.

The combination of lens imaging observables and the stellar kinematics is
sensitive to $\lambda_{\rm int}(1-\kappa_{\rm ext})D_{\rm s}/D_{\rm ds}$
(Birrer et al. 2016; Chen et al. 2021a). We apply a prior on $\kappa_{\rm
ext}$ using the estimated $\kappa_{\rm ext}$ distribution from Suyu et al.
(2014) to help break the degeneracy in distributing the total MSD into
external and internal components.

###  4.4 Bayesian framework

According to Bayes’ theorem, the posterior of the model parameters
$\Xi=\\{\xi_{\rm mass},\xi_{\rm light},D_{\Delta t}^{\rm model},i,\kappa_{\rm
ext},\lambda_{\rm int},D_{\rm d},\beta\\}$ as

$p(\Xi\mid\mathcal{D})\propto p(\mathcal{D}\mid\Xi)\ p(\Xi),$ (21)

where $p(\mathcal{D}\mid\Xi)$ is the likelihood given data $\mathcal{D}$ and
$p(\Xi)$ is the prior. In this study, the data $\mathcal{D}$ is the measured
velocity dispersions in Voronoi bins (Figure 9). The observational information
from the published time delays, lens models using HST imaging, and the line-
of-sight effects (Tewes et al. 2013; Suyu et al. 2013, 2014) is incorporated
by adopting those previous posteriors as the prior on our model parameters.
The likelihood of the observed velocity dispersion vector $\bm{\sigma}_{\rm
los}\equiv[\sigma_{1},\dots,\sigma_{N_{\rm bin}}]$, with $N_{\rm bin}$ being
the number of Voronoi bins, is given by

$\mathcal{L}(\bm{\sigma}_{\rm
los}\mid\Xi)\propto\exp\left[-\frac{1}{2}\bm{\sigma}_{\rm los}^{\rm
T}\Sigma^{-1}\bm{\sigma}_{\rm los}\right],$ (22)

where $\Sigma$ is the variance-covariance matrix. Specific priors used in this
Bayesian framework are given in Section 5. We obtain the posterior probability
distribution function (PDF) of the model parameters using the Markov-chain
Monte Carlo (MCMC) method using the affine-invariant ensemble sampler emcee
(Goodman & Weare 2010; Foreman-Mackey et al. 2013). We ensure the MCMC chains’
convergence by running the chains for $\gtrsim$20 times the autocorrelation
length after the chains have stabilized (Foreman-Mackey et al. 2013).

##  5 Dynamical models

We first describe our baseline dynamical model in Section 5.1 and then perform
various checks on systematics in Section 5.2.

###  5.1 Baseline dynamical model

This subsection describes the baseline settings in our dynamical model, namely
the specific parametrization of the mass model (Section 5.1.1), the dynamical
tracer profile (Section 5.1.2), the probability of oblate or prolate
axisymmetry (Section 5.1.3), the inclination angle (Section 5.1.4), and the
choice of anisotropy profile (Section 5.1.5).

####  5.1.1 Parametrization of the mass model

We adopt the power-law mass model as our baseline model. In this model, the
mass profile is defined with Einstein radius $\theta_{\rm E}$, logarithmic
slope $\gamma$, projected axis ratio $q_{\rm m}$, and position angle
$\varphi_{\rm mass}$. The convergence profile $\kappa_{\rm model}$ in Equation
(19) for the power-law model is given by

$\kappa_{\rm model}^{\rm
pl}(\theta_{1},\theta_{2})=\frac{3-\gamma}{2}\left(\frac{\theta_{\rm
E}}{\sqrt{q_{\rm m}\theta_{1}^{2}+\theta_{2}^{2}/q_{\rm
m}}}\right)^{\gamma-1}$ (23)

Here, the coordinates $(\theta_{1},\ \theta_{2})$ are rotated by $\varphi_{\rm
mass}$ from the (RA, Dec) coordinate system. We adopt the lens model posterior
from Suyu et al. (2013) as a prior in our dynamical model. For simplicity, we
set the position angle $\varphi_{\rm mass}$ the same as the observed position
angle of light $\varphi_{\rm light}$. We use the estimated $\kappa_{\rm ext}$
distribution for the power-law model as the prior (see Figure 3 of Suyu et al.
2014).

We set $\theta_{\rm s}=12\arcsec$ ($\simeq 7.5\theta_{\rm E}$) in the
approximate mass-sheet $\kappa_{\rm s}$ (Equation 12) so that the imaging
constraints alone cannot differentiate the power-law mass profile and its
approximate MST from Equation (11). We obtain this lower limit by running the
jupyter notebook that produces Figure 3 of (Birrer et al.
2020).101010https://github.com/TDCOSMO/hierarchy_analysis_2020_public/blob/6c293af582c398a5c9de60a51cb0c44432a3c598/MST_impact/MST_pl_cored.ipynb
However, we adjusted the fiducial lens model parameters in the notebook to
match with those for RXJ1131$-$1231\. We take a uniform prior for the internal
MST parameter $\lambda_{\rm int}\sim\mathcal{U}(0.5,1.13)$. The upper limit of
1.13 is set by the requirement that the transformed mass profile under the
approximate MST must be monotonic so that the MGE can approximate the
transformed profile sufficiently well (Shajib et al. 2019). Previous studies
also found similar or more restrictive upper limits for $\lambda_{\rm int}$ to
satisfy the physical requirement of non-negative density (Birrer et al. 2020;
Yıldırım et al. 2021).

The appropriate number of MGE components for the mass or light profile is
automatically chosen by jampy with a maximum of 20 components. We check that
the MGE approximates the input mass or light profile very well (with a maximum
1% deviation at $<10\arcsec$ and maximum 10% deviation between
10$\arcsec$–50$\arcsec$). These deviations from the density profile have an
oscillatory pattern due to the MGE approximation’s nature, except near the end
of the fitted ranges. Thus the deviation in the integrated mass profile often
averages out in the line-of-sight integration up to a very large radius. We
perform the MGE fitting up to 100$\arcsec$. Thus, the large mismatch between
the MGE approximation and the original profile occurs largely outside the
integration limit $\sim 70\arcsec$. The chosen number of maximum Gaussian
components is not a dominant source of numerical error. Setting this maximum
number to a very high value, such as 100, shifts the computed velocity
dispersion by only $<0.5$% within the observed region, which is insignificant
compared to the 1% numerical stability targeted by jampy.

####  5.1.2 Dynamical tracer profile

We update the light profile fitting for the lens galaxy from Suyu et al.
(2013) using a larger HST image cutout than that therein, which did not
contain the full extent of the lens galaxy’s light profile (see Figure 11).
The lensed arcs and quasar images are first subtracted from the cutout using
the prediction of the best-fit lens model from Suyu et al. (2013). We use the
software package lenstronomy111111https://github.com/lenstronomy/lenstronomy
to fit the residual light distribution attributed to the lens galaxy (Birrer &
Amara 2018; Birrer et al. 2021). Following Suyu et al. (2013), we use the
double Sérsic model to fit the light profile, which is a superposition of two
concentric Sérsic profiles. The Sérsic profile is defined as

$I(\theta_{1},\theta_{2})=I_{\rm 0}\exp\left[-b_{n}\left(\frac{\sqrt{q_{\rm
l}\theta_{1}^{2}+\theta_{2}^{2}/q_{\rm l}}}{\theta_{\rm eff}}\right)^{1/n_{\rm
s}}+b_{n}\right],$ (24)

where $I_{0}$ the amplitude, $q_{\rm l}$ is the axis ratio, $\theta_{\rm eff}$
is the effective radius, $n_{s}$ is the Sérsic index, and $b_{n}=1.999n-0.327$
is a normalizing factor so that $\theta_{\rm eff}$ becomes the half-light
radius (Sérsic 1968). The coordinates $(\theta_{1},\ \theta_{2})$ are rotated
by $\varphi_{\rm light}$ from the (Ra, Dec) coordinate system.

We first mask circular regions at the quasar image positions due to slightly
saturated pixels producing significant residuals in the subtracted cutout (see
Figure 11). We then iteratively mask the other pixels with significant
residuals above statistical expectations to effectively perform an outlier
rejection while preserving the shape of a Gaussian tail. For each iteration of
this process, we take a discrepancy threshold, which we decrease from
5$\sigma$ to 2$\sigma$ with step size 0.5$\sigma$ across these iterations. We
then randomly mask a subset of the pixels with residuals more than the
discrepancy level at the given iteration such that the number of remaining
pixels with such high residuals is statistically expected. The final masked
area after the iterations is illustrated in Figure 11. We tabulate the best-
fit light model parameters in Table 1 and compare them with those from Suyu et
al. (2013). The circularized half-light radius for our best-fit model is
$\theta_{\rm eff}=1\aas@@fstack{\prime\prime}91$, which is slightly larger
than the value $\theta_{\rm eff}=1\aas@@fstack{\prime\prime}85$ from Suyu et
al. (2013) based on the same imaging data but from a smaller cutout
(illustrated in 11). We then take the MGE of the fitted double Sérsic profile
as the light distribution $I(x,y)$ in our dynamical modeling. We propagate the
uncertainties and covariances from the light profile fitting into the
dynamical modeling. To do that, we sample from the multivariate normal
distribution corresponding to all the light model parameters for each call of
the likelihood function within the MCMC process and then take the MGE of the
light profile given the sampled parameters.

Table 1:  Values of the light model parameters for the double Sérsic model in our fitting of a large cutout and those from Suyu et al. (2013). The position angle $\varphi_{\rm light}$ is defined as East of North. Parameter | This analysis | Suyu et al. (2013)
---|---|---
Sérsic profile 1
$I_{0}$ (e-1 s-1 pixel-1) | 32.8 $\pm$0.1 | 36.4 $\pm$0.4
$\theta_{\rm eff}$ (″) | 2.437 $\pm$0.005 | 2.49 $\pm$0.01
$n_{\rm s}$ | 1.10 $\pm$0.01 | 0.93 $\pm$0.03
$q_{\rm l}$ | 0.865 $\pm$0.001 | 0.878 $\pm$0.004
Sérsic profile 2
$I_{0}$ (e-1 s-1 pixel-1) | 441 $\pm$7 | 356 $\pm$12
$\theta_{\rm eff}$ (″) | 0.300 $\pm$0.003 | 0.362 $\pm$0.009
$n_{\rm s}$ | 1.60 $\pm$0.02 | 1.59 $\pm$0.03
$q_{\rm l}$ | 0.847 $\pm$0.002 | 0.849 $\pm$0.004
$\varphi_{\rm light}$ (∘) | 120.5 $\pm$0.3 | 121.6 $\pm$0.5
Figure 11:  Fit of the lens galaxy’s surface brightness profile. Left: the
HST/ACS imaging in the F814W filter of the lens system RXJ1131$-$1231 with the
quasar images and the lensed arcs subtracted using the prediction from the
best-fit lens model from Suyu et al. (2013), thus leaving only the lens
galaxy’s light to be fitted. The orange circle shows the large circular region
considered for fitting in our analysis, and the yellow square shows the
smaller cutout used for lens modeling by Suyu et al. (2013). The cyan annulus
contains the region where pixels were fitted to reconstruct the source by Suyu
et al. (2013). Thus the lensed arcs from the quasar host galaxy were
subtracted only within this annulus. The red contours mark quasar image
positions with significant residuals due to saturated pixels, which we mask.
Middle: The fitted light profile with a double Sérsic model. The black pixels
correspond to masked pixels. The additional masked pixels within the orange
circle not described above are randomly selected through an iterative process
that performs outlier rejection while preserving the Gaussian tail (see
Section 5.1.2 for details). Right: Normalized residual of the best-fit model.

####  5.1.3 Oblate or prolate shape of the axisymmetry

The oblateness, prolateness, or triaxiality of a slow rotator galaxy can, in
principle, be constrained from the kinematic misalignment angle
$\Delta\varphi_{\rm kin}\equiv\left|\varphi_{\rm kin}-\varphi_{\rm
light}\right|$. However, we do not detect any significant rotational pattern
in the $v_{\rm mean}$ map (Figure 9). Thus, the uncertainty for the
constrained kinematic major axes is too large to be meaningful, and we cannot
directly constrain this galaxy’s oblateness from the data. Instead, we obtain
the probability of oblateness from a population prior based on 189 slow
rotator elliptical galaxies that are in the Sloan Digital Sky Survey’s
(SDSS’s) Mapping Nearby Galaxies at APO (MaNGA) sample (Abolfathi et al. 2018;
Graham et al. 2018). We take the distribution of $\Delta\varphi_{\rm kin}$ for
this sample of slow rotators (Li et al. 2018), where $\Delta\varphi_{\rm
kin}=0^{\circ}$ corresponds to a purely oblate shape, and $\Delta\varphi_{\rm
kin}=90^{\circ}$ corresponds to a purely prolate shape. Li et al. (2018) find
two distinct peaks in the distribution at $\Delta\varphi_{\rm kin}=0^{\circ}$
and $\Delta\varphi_{\rm kin}=90^{\circ}$ (see Figure 12). We, therefore, fit
the data points with a double Gaussian profile with the means set at
$\Delta\varphi_{\rm kin}=0^{\circ}$ and $\Delta\varphi_{\rm kin}=90^{\circ}$
(see the fit in Figure 12). Although the slow rotators with
$0^{\circ}<\Delta\varphi_{\rm kin}<90^{\circ}$ have triaxial shapes, we choose
only oblate or prolate axisymmetric shapes in our dynamical modeling for
computational simplicity. Therefore, we take $\Delta\varphi_{\rm
kin}<45^{\circ}$ as the oblate case and $\Delta\varphi_{\rm kin}>45^{\circ}$
as the prolate case. We obtain the prior probability $p({\rm oblate})_{\rm
pop}$ of the galaxy being oblate as

${{p({\rm oblate})_{\rm pop}=\int_{0^{\circ}}^{45^{\circ}}{{\rm
d}}(\Delta\varphi_{\rm kin})\ p(\Delta\varphi_{\rm kin})_{\rm pop}\simeq
0.65,}}$ (25)

and thus $p({\rm prolate})_{\rm pop}=1-p({\rm oblate})_{\rm pop}\simeq 0.35$.

The jampy software package, by default, adopts the oblate case for
deprojection. We implement the prolate case in jampy by setting $q_{\rm
prolate}=1/q>1$ and switching the $x$ and $y$ axes in the input coordinate
system. Due to the switching of $x$ and $y$ axes, $\sigma$ parameters of the
MGEs for mass and light models need to be scaled as $\sigma_{\rm
prolate}=q\sigma$.

Figure 12:  Population prior on kinematic misalignment angle
$\Delta\varphi_{\rm kin}\equiv\left|\varphi_{\rm kin}-\varphi_{\rm
light}\right|$ for a sample of slow rotator elliptical galaxies from the
SDSS’s MaNGA dataset (Li et al. 2018). Here, $\Delta\varphi_{\rm
kin}=0^{\circ}$ corresponds to a purely oblate shape, and $\Delta\varphi_{\rm
kin}=90^{\circ}$ corresponds to a purely prolate shape. The vertical dashed
grey lines mark $\Delta\varphi_{\rm kin}=0^{\circ}$, $45^{\circ}$, and
$90^{\circ}$. The red points with error bars show the measurements from Li et
al. (2018). We fit this distribution with a double Gaussian model (blue line)
with the means fixed to $\Delta\varphi_{\rm kin}=0^{\circ}$ and
$\Delta\varphi_{\rm kin}=90^{\circ}$. We take $\Delta\varphi_{\rm
kin}<45^{\circ}$ as the oblate case and $\Delta\varphi_{\rm kin}>45^{\circ}$
as the prolate case. Integrating the double Gaussian model from $0^{\circ}$ to
$45^{\circ}$ gives the prior probability of oblateness $p({\rm oblate})_{\rm
pop}\simeq 0.65$.

####  5.1.4 Inclination

The observed axis ratio of light $q_{\rm l,obs}=0.850\pm 0.002$ relates to
$q_{\rm l,int}$ through the inclination angle $i$ as

$q_{\rm l,obs}^{2}=q_{\rm l,int}^{2}\sin^{2}i+\cos^{2}i.$ (26)

We impose a prior on the intrinsic axis ratio $q_{\rm l,int}$ from a sample of
massive elliptical galaxies in the SDSS with stellar mass
$10.8<\log_{10}(M_{\star}/M_{\odot})<11.5$ at $0.04<z<0.08$ (Chang et al.
2013). The distribution of $q_{\rm l,int}$ by Chang et al. (2013) is different
for oblate and prolate assumptions. Therefore, we adopt the specific prior
corresponding to the oblate or the prolate case (see Figure 13).

Figure 13:  Prior on the intrinsic axis ratio $q_{\rm l,int}$ of light for
oblate (solid line) and prolate (dashed line) cases from Chang et al. (2013).
The priors correspond to massive elliptical galaxies from the SDSS survey at
$0.04<z<0.08$ with $10.8<\log_{10}(M_{\star}/M_{\odot})<11.5$.

####  5.1.5 Anisotropy profile

We investigate two choices to parametrize the anisotropy profile. The first
choice is a single spatially constant
$\beta=1-\sigma_{\theta}^{2}/\sigma_{r}^{2}$ value for all the light MGE
components. Numerically, we sample $\sigma_{\theta}/\sigma_{r}$ with a uniform
prior $(\sigma_{\theta}/\sigma_{r})\sim\mathcal{U}(0.78,1.14)$. This range of
$\sigma_{\theta}/\sigma_{r}$ allows $-0.31<\beta<0.38$. We adopt this range
using the $\beta$ values of eight slow rotator galaxies measured by Cappellari
et al. (2007, see Figure 2). These measurements of $\beta$ by Cappellari et
al. (2007) are from Schwarzschild modeling of data with one of the highest
$S/N$ values in the literature, allowing to constrain the Gauss–Hermite
moments up to order six. Applying the student’s $t$-distribution on the sample
mean of this small sample, we find the 95% confidence interval of the
population mean for $\beta$ to be [-0.10, 0.17] and the standard deviation to
be 0.16. These values infer that 95% of the population is contained within
$\beta\in[-0.31,0.38]$, which we take as the boundaries of our prior range.
The second choice of the anisotropy profile has two free parameters: the inner
light MGE components with $\sigma<r_{\rm break}=\theta_{\rm
eff}=1\aas@@fstack{\prime\prime}91$ are assigned one value for
$(\sigma_{\theta}/\sigma_{r})_{\rm inner}$ and the outer light MGE components
with $\sigma\geq r_{\rm break}$ are assigned another independent value of
$(\sigma_{\theta}/\sigma_{r})_{\rm outer}$. Thus, this parametrization with
two free parameters allows radial variability in the anisotropy profile. Both
the inner and outer ratios have uncorrelated uniform priors
$(\sigma_{\theta}/\sigma_{r})\sim\mathcal{U}(0.78,1.14)$. For these two
choices of parametrization, we compute the Bayesian information criterion
(BIC) given by

${\rm BIC}\equiv k\log(N_{\rm bin})-2\log\hat{\mathcal{L}},$ (27)

where $k$ is the number of free model parameters, $N_{\rm bin}$ is the number
of data points, and $\hat{\mathcal{L}}$ is the maximum likelihood. We
approximate $\hat{\mathcal{L}}$ from the highest likelihood value sampled in
the MCMC chain. The single-parameter $\beta$ model provides the lowest BIC
value excluding the two-parameter $\beta$ model with $\Delta{\rm BIC}\approx
3.7$ (i.e., positively excluded; Raftery 1995). We check that the difference
between the highest and the second highest likelihood values among the MCMC
samples is $\ll\Delta$BIC, thus this $\Delta$BIC value is robust against our
approximation of $\hat{\mathcal{L}}$ from the highest likelihood value in the
sampled chain. The non-detection of varying anisotropy in our data is
consistent with that observed in nearby elliptical galaxies, as even
high-$S/N$ SAURON data for a large sample of galaxies are accurately described
by JAM models with constant anisotropy, as used here, within the noise of the
kinematics (e.g., Cappellari et al. 2013). We compare the posterior
distributions of the model parameters for the two anisotropy models in Figure
14. An example of a best-fit kinematic model and the corresponding residual
with the single-parameter $\beta$ model and oblate axisymmetry is illustrated
in Figure 15. The reduced $\chi_{\nu}^{2}$ value is 0.83 with $\nu=41$ degrees
of freedom. The distribution of residuals is similar to a normal distribution
expected from a perfect model for data with Gaussian noise, illustrating that
our model is appropriate for the data. We show the range of velocity
dispersion radial profiles sampled by our model in Figure 16 and compare it
with the radially averaged measurements of the velocity dispersion. This
illustration shows that our model reproduces the uncertainty range of the
measurement.

Figure 14:  Constraints from axisymmetric JAM modeling on the power-law mass
model parameters ($\theta_{\rm E}$, $\gamma$, and $q_{\rm m}$), internal MST
parameter $\lambda_{\rm int}$, external convergence $\kappa_{\rm ext}$,
anisotropy profile parameter(s), and the cosmological distances $D_{\Delta t}$
and $D_{\rm d}$. assuming two anisotropy parametrizations: (i) one single
constant $\beta\equiv 1-(\sigma_{\theta}/\sigma_{r})^{2}$ for all light MGE
components (orange contours), and (ii) one free
$(\sigma_{\theta}/\sigma_{r})_{\rm inner}\equiv(\sigma_{\theta}/\sigma_{r})$
for light MGE components with $\sigma<r_{\rm break}=\theta_{\rm
eff}=1\aas@@fstack{\prime\prime}91$ and another free
$(\sigma_{\theta}/\sigma_{r})_{\rm outer}$ for light MGE components with
$\sigma>r_{\rm break}$(blue contours). The blinded parameters are blinded as
$p_{\rm blinded}\equiv p/\langle p\rangle-1$ so that the distributions only
reveal fractional uncertainties. The darker and lighter shaded regions in the
2D plots trace 68% and 95% credible regions, respectively. The mass model
parameters Einstein radius $\theta_{\rm E}$, power-law slope $\gamma$, axis
ratio $q$, and position angle PA are additionally constrained through a prior
from the imaging data from Suyu et al. (2013). The two anisotropy
parametrizations provide equally good fits to the kinematics data. However,
the BIC selects the constant-$\beta$ anisotropy model over the other one with
one additional free parameter ($\Delta$BIC value is 3.5).  Figure 15:
Observed velocity dispersion map in Voronoi bins (first panel), the best-fit
dynamical model with a power-law mass model, constant $\beta$ anisotropy
profile, and oblate shape (second panel), the normalized residual for the
best-fit dynamical model (third panel), and the distribution of the normalized
residual (orange, fourth panel). The reduced $\chi^{2}$ quantity is
${\chi^{2}_{\nu}}=0.83$ with degrees of freedom $\nu=41$. The grey dashed line
in the fourth panel shows a normal distribution expected for residuals from a
perfect model to the data with Gaussian noise. The residual distribution for
41 points is similar to this Gaussian distribution.  Figure 16: Radial profile
of the line-of-sight velocity dispersion. The red points are radially binned
values from the 2D maps, with the horizontal error bars illustrating the
widths of the annuli. The lines show the radial profiles for random samples
from the dynamical model posterior. The radial profile of the model is
averaged over the major, minor, and intermediate axes. The solid purple lines
correspond to 65 random samples for the oblate case, and the dashed green
lines correspond to 35 random samples for the prolate case. Note that the
model was fit to the 2D kinematics data. However, we illustrate the 1D radial
profile only for visualization.

###  5.2 Checking potential systematics due to modeling choices

In this section, we perform several checks on potential systematics for
different choices in the dynamical model setup.

####  5.2.1 Comparison between power-law and composite mass models

In addition to the power-law mass model, Suyu et al. (2014) also adopted a
composite mass model individually describing the lens galaxy’s dark matter and
baryonic components. The dark matter distribution was modeled with an
elliptical NFW profile in the potential. The parameters in this profile are
the normalization of the NFW component $\kappa_{\rm s}$, the NFW scale radius
$r_{\rm scale}$, and the mass axis ratio $q_{\rm m}$. The baryonic component
was modeled with a mass-follow-light profile with a free mass-to-light ratio
($M/L$) parameter. Thus, this mass model parametrization has one more free
parameter than the power-law model. See Suyu et al. (2014) for parametric
definitions of these profiles. We implement this composite mass profile as
$\kappa_{\rm model}^{\rm comp}$ in Equation (19) and adopt the model posterior
from Suyu et al. (2014) as a prior in our model. We appropriately convert the
ellipticity defined in the potential by Suyu et al. (2014) to an ellipticity
defined in the convergence in our model. We take the MGE of this composite
surface density model as done for the power-law surface density model.
However, since the dark matter and baryonic components have different
ellipticities, we take the MGE of each component separately to preserve the
ellipticity information in deprojection. Specifically, We take the MGE of the
approximate MST with $\lambda_{\rm int}$ of the dark matter profile and the
MGE of an accordingly rescaled baryonic profile, which effectively results in
the total mass profile being transformed as the approximate MST with
$\lambda_{\rm int}$.

This mass model with one more free parameter than the power-law model has a
higher BIC score with $\Delta$BIC $=3.8$. Thus, the BIC excludes the composite
model with positive evidence (Raftery 1995). The median values of $D_{\rm d}$
from the power-law and composite mass models differ by 0.9% ($0.07\sigma$,
Figure 17), and the median $D_{\Delta t}$ values differ by 1.26%
($0.06\sigma$). Therefore, we conclude that our power-law mass model with an
additional degree of freedom to scale with the MST robustly describes the
observed data.

Figure 17:  Comparison of the constrained $D_{\rm d}$ from power-law (blue
contours) and composite (orange contours) mass models. The blinded parameters
are blinded as $p_{\rm blinded}\equiv p/\langle p\rangle-1$ so that the
distributions only reveal fractional uncertainties. The darker and lighter
shaded regions in the 2D plots trace 68% and 95% credible regions,
respectively.

####  5.2.2 Comparison between prolate and oblate axisymmetry

We compare the inferred $D_{\rm d}$ between the purely oblate and purely
prolate cases in the deprojected 3D spheroidal shape of the mass and light
models (Figure 18). The median $D_{\rm d}$ values from these two cases differ
by 3.6% ($0.3\sigma$), and the median $D_{\Delta t}$ values differ by 0.94%
($0.04\sigma$). Our final distance posterior is the combination of oblate and
prolate cases, with weights $p({\rm oblate})_{\rm pop}=0.65$ and $1-p({\rm
oblate})_{\rm pop}=0.35$, respectively. Thus, this difference between the
oblate and prolate cases is marginalized in our final cosmological distance
posterior.

We also compare the predictions from axisymmetric and spherical mass models in
Figure 18. The median $D_{\rm d}$ from the spherical model matches very well
with the axisymmetric prolate model, but the median $D_{\Delta t}$ differs by
2.0% ($0.08\sigma$). The galaxy is only mildly elliptical in projection
($q_{\rm l}\sim 0.85$), and the resulting axisymmetric models are not very
flat. For this reason, the relatively small difference between the
axisymmetric and spherical models is not surprising.

Figure 18:  Comparison of the constrained $D_{\rm d}$ between oblate (blue)
and prolate (blue) cases of the deprojected spheroidal shape in the dynamical
model. The blinded parameters are blinded as $p_{\rm blinded}\equiv p/\langle
p\rangle-1$ so that the distributions only reveal fractional uncertainties.
The darker and lighter shaded regions in the 2D plots trace 68% and 95%
credible regions, respectively.

####  5.2.3 Comparison between Voronoi binning schemes

Here, we compare the Voronoi binning schemes with two choices for the target
$S/N$ in each bin: $\approx 23$ Å-1 and $\approx 28$ Å-1. The two cases match
very well with only a 0.21% difference ($0.02\sigma$) in the median values of
$D_{\rm d}$ (Figure 19) and 0.28 % difference (0.01$\sigma$) in the median
$D_{\Delta t}$ values. As a result, we conclude that our choice of the Voronoi
binning scheme is not a significant source of systematic error in our
analysis.

Figure 19:  Comparison of the constrained $D_{\rm d}$ between two choices of
the target $S/N$ for each bin in the Voronoi binning scheme. The blinded
parameters are blinded as $p_{\rm blinded}\equiv p/\langle p\rangle-1$ so that
the distributions only reveal fractional uncertainties. The darker and lighter
shaded regions in the 2D plots trace 68% and 95% credible regions,
respectively.

Based on the systematics tests performed above, we adopt a robust final
distance posterior from the model with the power-law parametrization for the
mass profile that the approximate internal MST is applied to. We marginalize
the oblate and prolate axisymmetrical cases by combining the posteriors from
these two choices with weights of 0.65 and 0.35, respectively. In the next
section, we present the unblinded values from the distance posterior and infer
the value of $H_{0}$ from it.

##  6 Cosmological inference

In this section, we infer cosmological parameters from the joint distribution
of $D_{\rm d}$ and $D_{\Delta t}$, accounting for their covariance. The
unblinded point estimates of these distances are $D_{\rm d}=865_{-81}^{+85}$
Mpc (a 9.6% measurement) at $z_{\rm d}=0.295$, and $D_{\Delta
t}=2180_{-271}^{+472}$ Mpc (a 17% measurement) for $z_{\rm s}=0.657$.

We infer $H_{0}$ and $\Omega_{\rm m}$ from our distance posterior for a flat
$\Lambda$CDM cosmology (see Figure 20, left panel). We leave the exploration
of more exotic cosmologies based on our distance posterior for future studies.
We approximate the likelihood function $\mathcal{L}(H_{0},\Omega_{\rm m}\mid
D_{\rm d},D_{\Delta t})$ of the cosmological parameters using a 2D Gaussian
kernel density estimate (KDE) from the 2D distance posterior. We adopt two
choices of prior for $\Omega_{\rm m}$: one is a uniform prior $\Omega_{\rm
m}\sim\mathcal{U}(0.05,\ 0.5)$, and the other is a Gaussian prior $\Omega_{\rm
m}\sim\mathcal{N}(0.334,\ 0.018)$ from the Pantheon$+$ analysis of type Ia
supernovae relative distances (Brout et al. 2022). We infer the posterior
joint PDF of $H_{0}$ and $\Omega_{\rm m}$ by performing MCMC sampling using
emcee, given the likelihood function and prior choice.

We infer $H_{0}$ = $77.1_{-7.1}^{+7.3}$ km s-1 Mpc-1(a 9.4% measurement) with
the uniform $\Omega_{\rm m}$-prior, and $H_{0}$ = $76.0_{-6.6}^{+7.3}$ km s-1
Mpc-1 (a 9.1% measurement) with the Pantheon$+$ $\Omega_{\rm m}$-prior (solid
contours in the right panel of Figure 20). We show the $D_{\Delta t}$–$D_{\rm
d}$ region allowed by our priors in the left panel of Figure 20, which also
shows the region allowed by our distance posterior that provides information
for the cosmological inference. Other cosmological models beyond flat
$\Lambda$CDM (e.g., Bonvin et al. 2017; Wong et al. 2020) or combining other
cosmological probes in a cosmology-independent manner (e.g., Taubenberger et
al. 2019) can utilize the additional cosmological information contained by our
full 2D posterior outside the regions probed by our cosmological priors.

For comparison, we also perform cosmological inference using only the 1D
posterior of $D_{\rm d}$ (dashed contours in right panel of Figure 20). This
gives $H_{0}=75.5_{-7.2}^{+8.3}$ km s-1 Mpc-1 (a 10.3% measurement) for the
uniform $\Omega_{\rm m}$-prior, and $H_{0}=74.4_{-6.2}^{+8.1}$ km s-1 Mpc-1 (a
9.6% measurement) for the Pantheon$+$ $\Omega_{\rm m}$-prior. The $D_{\rm
d}$-only constraints are lower by $\sim$1.4% ($0.15\sigma$) than that from the
full 2D distance posterior (for the uniform $\Omega_{\rm m}$-prior). This
slight difference arises from the projection difference of the 2D posterior
along the $D_{\rm d}$ direction and along the narrow track allowed by our
choice of cosmological priors.

Figure 20:  Left: Final 2D posterior of the time-delay distance $D_{\Delta t}$
and the angular diameter distance $D_{\rm d}$ (emerald contour). The darker
and lighter shaded regions in the 2D plots trace 64% and 95% credible regions,
respectively. We infer $H_{0}$ and $\Omega$ from this distance posterior
accounting for the covariance in a flat $\Lambda$CDM cosmology. We take a wide
uniform prior on $H_{0}\sim\mathcal{U}(0,150)$ km s-1 Mpc-1. The blue-shaded
region corresponds to a uniform prior $\Omega_{\rm m}\sim\mathcal{U}(0.05,\
0.5)$ and the orange-shaded region corresponds to a Gaussian prior
$\Omega_{\rm m}\sim\mathcal{N}(0.334,\ 0.018)$ from the Pantheon$+$ analysis
of type Ia supernovae relative distances (Brout et al. 2022). Right: Posterior
PDF of $H_{0}$ and $\Omega_{\rm m}$ in flat $\Lambda$CDM cosmology. We
constrain $H_{0}$ to 9.4% and 9.1% precision for the uniform and Pantheon$+$
$\Omega_{\rm m}$-priors, respectively. We show the cosmological parameter
posterior from only the 1D $D_{\rm d}$ posterior with dashed contours with
colors matching the associated $\Omega$ prior. In this case, the $H_{0}$
precision is 10.3% and 9.6% for the uniform and Gaussian priors, respectively.
The $D_{\rm d}$-only constraint on the $H_{0}$ is lower by $\sim$1.4%
(0.15$\sigma$) than the constraint from the full 2D posterior, for the uniform
$\Omega_{\rm m}$-prior.

##  7 Discussion

We now compare our results with previous works (Section 7.1), discuss the
improvement of the constraint in this paper over single-aperture stellar
kinematics (Section 7.2), and describe the limitations of this work (Section
7.3).

###  7.1 Comparison with previous time-delay $H_{0}$ measurements

Our measured value $H_{0}=77.1_{-7.1}^{+7.3}$ km s-1 Mpc-1 is consistent with
previous measurements from lensing time delays with different treatments of
the MSD. These previous studies can be divided into two approaches: the first
breaks the MSD by assuming simple parametric mass profiles such as the power
law or composite (i.e., NFW halo and stars with constant mass-to-light ratio),
and the second breaks the MSD based solely on stellar kinematics. Our study
belongs to the second approach by allowing the freedom in the model to be
maximally degenerate with $H_{0}$ and constraining it solely from the
spatially resolved stellar kinematics. However, it is illustrative to compare
our result with the first approach to discuss the validity of their mass model
assumptions.

Following the first approach, Suyu et al. (2013, 2014) measured $H_{0}$ =
$80.0_{-4.7}^{+4.5}$ km s-1 Mpc-1 from this same system RXJ1131$-$1231 with
simple parametric mass profiles using HST imaging. Chen et al. (2019) combined
the HST imaging and adaptive-optics-assisted imaging from the Keck Telescope
to measure $H_{0}$ = $78.3^{+3.4}_{-3.3}$ km s-1 Mpc-1Ȧlthough these studies
used single-aperture stellar kinematics, the MSD was already broken by the
assumption of parametric mass profiles, and the single-aperture velocity
dispersion helped tighten the constraint and made the inferred $H_{0}$ values
from the power-law and composite models more consistent(Suyu et al. 2014). Our
measured value – albeit with a larger uncertainty due to the maximal freedom
allowed in the mass model – has a median value very close to these previous
measurements. Such a good agreement in the medians suggests that these
previous studies’ simple parametric mass models are close to the ground truth,
and no bias is detected within the precision afforded by the data. Future
spatially resolved velocity dispersion measurements for more time-delay lens
systems or better quality data for this system (e.g., from the James Webb
Space Telescope) will allow us to make a more definitive statement on the
validity of the parametric mass model assumptions.

Following the second approach, Birrer et al. (2016) analyzed this same system
RXJ1131$-$1231 using HST imaging and single-aperture velocity dispersion.
These authors marginalized the effect of MSD by incorporating a source on the
prior but found that the $H_{0}$ posterior strongly depends on the shape of
the anisotropy prior. These authors use two different choices for this prior
to find $H_{0}$ = $74.5_{-7.8}^{+8.0}$ km s-1 Mpc-1 and $H_{0}$ =
$86.6_{-6.9}^{+6.8}$ km s-1 Mpc-1. This large difference illustrates that
single aperture velocity dispersion imposes only a weak constraint on the
anisotropy profile and, thus, on the MSD. This result highlights the need for
spatially resolved velocity dispersion, such as the one presented in this
study. Our measured $H_{0}$ has a precision of 9% while allowing the data to
constrain the MSD effect that is maximally degenerate with $H_{0}$,
illustrating the power of spatially resolved kinematics in constraining the
anisotropy profile and the MSD, despite the seeing-limited nature of our data.
In the future, exquisite data from the James Webb Space Telescope (JWST) will
provide an even more dramatic improvement (4% $H_{0}$ precision forecasted,
Yıldırım et al. 2021).

We also compare our result with the measured values of $H_{0}$ from the
current TDCOSMO sample of seven time-delay lenses. With the power-law mass
model assumptions, the combination of seven time-delay lenses gives a 2%
measurement with $H_{0}$ = $74.2_{-1.6}^{+1.6}$ km s-1 Mpc-1 (Wong et al.
2020; Millon et al. 2020b). However, relaxing this mass profile assumption and
constraining the MSD solely from the single-aperture stellar kinematics of the
TDCOSMO sample leads to a 9% uncertainty on the resultant $H_{0}$ =
$74.5^{+5.6}_{-6.1}$ km s-1 Mpc-1. In this study, we achieve the same 9%
precision from a single system, highlighting the superb constraining power of
spatially resolved kinematics over single-aperture ones.

It is also worth comparing with the result obtained by Birrer et al. (2020)
when combining the seven TDCOSMO lenses with information obtained from the
external SLACS sample of non-time-delay lenses, $H_{0}$ =
67.4${}_{-3.2}^{+4.1}$ km s-1 Mpc-1. Given the uncertainties, our new
measurement is not statistically inconsistent with that result, although the
difference is clearly important from a cosmological standpoint. With the data
in hand, we cannot conclude whether (a) the difference is real and the SLACS
sample cannot, therefore, be combined with the TDCOSMO sample, or whether (b)
it is due to a statistical fluctuation. This study demonstrates that, as we
gather more and better data for spatially resolved kinematics and external
samples of non-lenses, we will soon be able to conclude whether the difference
is real or not.

In the context of the “Hubble tension”, our new measurement strengthens the
tension by reaffirming the previously obtained time-delay $H_{0}$ measurements
that agreed with other local measurement values, e.g., from SH0ES (Riess et
al. 2022). Although the 9% uncertainty in $H_{0}$ from our measurement alone
is not sufficient to resolve the tension, it demonstrates that time-delay
cosmography can provide a powerful independent perspective with the help of
future data from telescopes such as Keck, JWST, and the extremely large
telescopes (e.g., see forecasts from Shajib et al. 2018; Yıldırım et al. 2021;
Birrer & Treu 2021). We cannot help noticing that the median of our
measurement is somewhat higher than the mean of the local values ($\sim$73 km
s-1 Mpc-1). However, the difference is not significant, given the
uncertainties. Therefore our likely explanation is that the difference
originates from statistical fluctuation pertaining to this system, as the
initial $H_{0}$ measurements using simple parametric assumption all provided
such higher values (Suyu et al. 2013, 2014; Birrer et al. 2016; Chen et al.
2019). We conclude by stressing that some dispersion around the mean is, of
course, expected, and indeed Millon et al. (2020b) shows that the seven
TDCOSMO lenses scatter around the mean by an amount consistent with the
estimated errors.

Figure 21:  Comparison of our 9.4% $H_{0}$ measurement (red,
$77.1_{-7.1}^{+7.3}$ km s-1 Mpc-1) from the single system RXJ1131$-$1231 with
previous measurements from Chen et al. (2019, blue), Millon et al. (2020b,
grey), and Birrer et al. (2020, emerald). The distributions show the $H_{0}$
posteriors as described in the figure legend, and the points with error bars
mark the mean and 64% credible intervals of the corresponding posterior with
matching color. For the same flexible mass models, our analysis on a single
system provides a similar precision on $H_{0}$ with that from seven lenses
with only single-aperture stellar kinematics (emerald, $73.3^{+5.8}_{-5.8}$km
s-1 Mpc-1). Moreover, the median value of our measurement falls very close to
those from previous analyses on the same system but with simple parametric
assumption on the mass model breaking the MSD (blue, $78.3^{+3.4}_{-3.3}$km
s-1 Mpc-1, cf. also $80.0_{-4.7}^{+4.5}$ km s-1 Mpc-1by Suyu et al. 2014).

###  7.2 Improvement from the spatial resolution of the stellar kinematics

We investigate the improvement in constraints provided by the spatially
resolved nature of the stellar kinematics presented in this paper over the
unresolved or single-aperture case. Suyu et al. (2013) presents a single-
aperture measurement of the line-of-sight velocity dispersion $\sigma_{\rm
los}=323\pm 20$ km s-1 obtained within a $0\aas@@fstack{\prime\prime}81\times
0\aas@@fstack{\prime\prime}7$ aperture with a $0\aas@@fstack{\prime\prime}7$
seeing. This measurement was from the Low-Resolution Imaging Spectrometer
(LRIS; Oke et al. 1995) on the Keck Observatory. The probed wavelength range
was $\sim$3900–4700 Å, which probes mostly the redward range of the Ca H&K
lines with a little overlap with the range probed by our data (i.e., 3300–4200
Å). If we take a luminosity-weighted-sum of the spatially resolved velocity
dispersion map within the same $0\aas@@fstack{\prime\prime}81\times
0\aas@@fstack{\prime\prime}7$ aperture, we get $288\pm 5$ km s-1, which is
$1.7\sigma$ (11%) lower than the previous single-aperture measurement.
Although the $1.7\sigma$ difference is not statistically significant, some
parts of it can be due to potential systematics in the kinematic extraction
procedure or due to different wavelength ranges probed. It is generally
considered that the minimum error, considering systematics, on velocity
dispersion measurements is 5%, even for very high-$S/N$ data.

However, to illustrate the improvement in precision from the spatially
resolved nature of the velocity dispersion presented in this study, we take a
fiducial single-aperture measurement value of $288\pm 18$ km s-1. This mean
value is from the luminosity-weighted sum within the single aperture mentioned
above, and the 18 km s-1 uncertainty comes from applying the 6% uncertainty of
the $323\pm 20$ km s-1 measurement on the fiducial mean. We take the galaxy’s
major axis to align with the rectangular aperture’s longer side. Rotating the
aperture by $90\degr$ only changes the predicted velocity dispersion
integrated within the aperture by $\lesssim 0.1$%, which is unsurprising given
the mild ellipticity ($q_{\rm l}\sim 0.85$) of the galaxy and the
$0\aas@@fstack{\prime\prime}96$ seeing. We compare the key dynamical model
parameters between the spatially resolved and single-aperture cases in Figure
22. As expected, the internal MST parameter $\lambda_{\rm int}$ and the
anisotropy profile parameter $\sigma_{\theta}/\sigma_{r}$ are almost
completely unconstrained in the case of the single-aperture stellar kinematics
due to the mass-anisotropy degeneracy (Treu & Koopmans 2002; Courteau et al.
2014). However, the angular diameter distance $D_{\rm d}$ can be constrained
to 15.7% precision, largely by the anisotropy prior (cf. the 9.6% constraint
on $D_{\rm d}$ from the spatially resolved data). This single-aperture
precision level on $D_{\rm d}$ agrees very well with the 17.9% precision on
$D_{\rm d}$ ($=810_{-130}^{+160}$ Mpc) obtained by Jee et al. (2019) from the
same system RXJ1131$-$1231 based on the previously available single-aperture
stellar kinematics mentioned above. The Hubble constant $H_{0}$ can be
inferred to 12.5% precision with the uniform $\Omega_{\rm m}$ prior from the
full 2D posterior of the fiducial single-aperture case. Although the
improvement in $H_{0}$ precision (by $\sim$3%) from the spatially resolved
kinematics does not appear to be dramatic, this is due to the fact that the
projection of $D_{\Delta t}$–$D_{\rm d}$ posterior along the narrow track
allowed by our chosen prior happens to give a small difference between the two
cases. The improvement could have appeared more drastic if the full 2D
posterior had a different orientation from the prior region. In reality, the
full cosmological information (illustrated by the area enclosed within the 95%
contour) contained by the single-aperture data is much diluted than that from
the spatially resolved data presented in this study (see the $D_{\Delta
t}$–$D_{\rm d}$ contours in Figure 22).

Figure 22:  Comparison of the distance constraints between spatially resolved
velocity dispersion and single-aperture velocity dispersion. Here, the
integrated velocity dispersion is taken as the fiducial value of $287\pm 18$
km s-1 to match the mean of our spatially resolved measurement, but the
uncertainty of a single-aperture velocity dispersion measurement (Suyu et al.
2013). The darker and lighter shaded regions in the 2D plots trace 68% and 95%
credible regions, respectively. The single-aperture velocity dispersion cannot
constrain the anisotropy profile parameter $\sigma_{\theta}/\sigma_{r}$ and
the internal MST parameter $\lambda_{\rm int}$, with both limited by the
prior. As a result, the $D_{\Delta t}$–$D_{\rm d}$ posterior is constrained
much more weakly.

###  7.3 Limitations of this study

One limitation of our study is the data quality. Although our data are the
first of their kind from a cutting-edge ground-based facility such as the Keck
Observatory, there are opportunities to obtain better-quality data. The KCWI
instrument is seeing-limited. Thus the $S/N$ on the lensing galaxy is degraded
by contamination from the nearby quasars, and the spatial resolution of the
velocity dispersion map is limited by the seeing. Adaptive-optics-assisted IFU
spectroscopy from the ground or observations from space, e.g., with the JWST,
can deliver exquisite spatially resolved data for improved $H_{0}$ precision
in the future (Yıldırım et al. 2020, 2021).

Future data with higher spatial resolution will be particularly powerful in
constraining the anisotropy profile better. Our measurement has only weak
constraints on the anisotropy profile, which is largely bounded by the adopted
uniform prior (see Figure 14). This prior is obtained from a sample of eight
local massive ellipticals with one of the highest quality spatially resolved
kinematics. However, this is a small sample size. A tighter anisotropy prior
from larger samples of massive ellipticals, even better if they are from a
redshift range that matches with the one for our system, will be helpful to
mitigate further the degeneracy induced by the anisotropy profile, i.e., the
mass-anisotropy degeneracy (Treu & Koopmans 2002; Courteau et al. 2014).

##  8 Conclusion

We measured the spatially resolved stellar velocity dispersion of the lens
galaxy in RXJ1131$-$1231 using the KCWI IFU spectrograph on the Keck
Observatory. We combined the new spatially resolved stellar kinematics with
previously obtained lens models derived from HST imaging data, observed time
delays, and estimated line-of-sight lensing effects (i.e., the external
convergence) to infer $H_{0}$. Combining the spatially resolved velocity
dispersion with lens imaging and time delays simultaneously alleviates the MSD
in the measured $D_{\Delta t}$ and additionally measures the angular diameter
distance $D_{\rm d}$.

In order to prevent conscious or unconscious experimenter bias, we blindly
performed the dynamical modeling and the cosmographic inference. We unblinded
the $H_{0}$ value after all the co-authors had agreed on the modeling choices
after various checks on systematics, and the analysis was frozen. The main
conclusions from our study are as follows:

  * •

The 2D distance posterior of $D_{\rm d}$ and $D_{\Delta t}$ gives
$H_{0}=77.1_{-7.1}^{+7.3}$ km s-1 Mpc-1 for a uniform prior on $\Omega_{\rm
m}\sim\mathcal{U}(0.05,0.5)$, and $H_{0}=76.0_{-6.6}^{+7.3}$ km s-1 Mpc-1for a
Gaussian prior on $\Omega_{\rm m}$ from the Pantheon$+$ analysis (Brout et al.
2022).

  * •

Our 9.4% measurement from a single system with spatially resolved kinematics
provides a similar precision as, and is in excellent agreement with, the
current TDCOSMO sample of seven time-delay lenses based only on single-
aperture stellar kinematics ($H_{0}$ = $74.5^{+5.6}_{-6.1}$ km s-1 Mpc-1,
Birrer et al. 2020). Note that the system RXJ1131$-$1231 analyzed here is part
of that sample of seven.

  * •

The median value of $H_{0}$ from our analysis is very close to the previously
inferred values assuming simple parametric mass models (e.g., $H_{0}$ =
$78.3^{+3.4}_{-3.3}$ km s-1 Mpc-1, Chen et al. 2019). Thus we do not detect
any potential bias in those mass profile assumptions within the precision
afforded by our data.

In conclusion, our study provides an important validation of previous work by
our collaboration on the determination of $H_{0}$ from time-delay cosmography.
This analysis also showcases the power of spatially resolved kinematics in
breaking the degeneracies that limit the $H_{0}$ precision when mass profile
assumptions on the galaxy density profile are relaxed. As the first
application of such methodology performed on real data, this study stands as
an important proof of concept to pioneer future studies on many more time-
delay lens systems. A future sample of $\sim$40 lensed quasars can
independently provide $\sim$1.2% precision on $H_{0}$ that is necessary to
resolve or confirm the “Hubble tension” at 5$\sigma$ confidence level, thanks
to spatially resolved stellar kinematics (Birrer et al. 2021).

###### Acknowledgements.

We thank Elizabeth Buckley-Geer, Thomas E. Collett, Philip J. Marshall, and
Chiara Spiniello for useful discussions and comments that improved this study
and the manuscript. Support for this work was provided by NASA through the
NASA Hubble Fellowship grant HST-HF2-51492 awarded to AJS by the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. TT and GCFC acknowledge support by NSF through grants NSF-
AST-1906976 and NSF-AST-1836016, and from the Moore Foundation through grant
8548. SHS thanks the Max Planck Society for support through the Max Planck
Research Group and the Max Planck Fellowship. SHS is supported in part by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311. This project has
received funding from SNSF and from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme
(COSMICLENS : grant agreement No 787886) VNB gratefully acknowledges
assistance from National Science Foundation (NSF) Research at Undergraduate
Institutions (RUI) grant AST-1909297. Note that findings and conclusions do
not necessarily represent views of the NSF.
This work used computational and storage services associated with the Hoffman2
Shared Cluster provided by UCLA Institute for Digital Research and Education’s
Research Technology Group.
The data presented herein were obtained at the W. M. Keck Observatory, which
is operated as a scientific partnership among the California Institute of
Technology, the University of California and the National Aeronautics and
Space Administration. The Observatory was made possible by the generous
financial support of the W. M. Keck Foundation. The authors wish to recognize
and acknowledge the very significant cultural role and reverence that the
summit of Maunakea has always had within the indigenous Hawaiian community. We
are most fortunate to have the opportunity to conduct observations from this
mountain.
This research made use of jampy (Cappellari 2008, 2020), pPXF (Cappellari
2017, 2022), pafit (Krajnović et al. 2006), vorbin (Cappellari & Copin 2003),
mgefit (Cappellari 2002) , lenstronomy (Birrer & Amara 2018; Birrer et al.
2021), numpy (Oliphant 2015), scipy (Jones et al. 2001), astropy (Astropy
Collaboration 2013, 2018), jupyter (Kluyver et al. 2016), matplotlib (Hunter
2007), seaborn (Waskom et al. 2014), emcee (Foreman-Mackey et al. 2013), and
getdist (https://github.com/cmbant/getdist).

## References

  * Abdalla et al. (2022) Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, Journal of High Energy Astrophysics, 34, 49
  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, The Astrophysical Journal Supplement Series, 235, 42
  * Aiola et al. (2020) Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmology Astropart. Phys., 2020, 047–047
  * Astropy Collaboration (2013) Astropy Collaboration. 2013, A&A, 558, A33
  * Astropy Collaboration (2018) Astropy Collaboration. 2018, AJ, 156, 123
  * Avila et al. (2015) Avila, R., Koekemoer, A., Mack, J., & Fruchter, A. 2015, Optimizing pixfrac in Astrodrizzle: An example from the Hubble Frontier Fields, Tech. rep.
  * Bacon et al. (1983) Bacon, R., Simien, F., & Monnet, G. 1983, Astronomy and Astrophysics, Vol. 128, p. 405-410 (1983), 128, 405
  * Barnabè et al. (2009) Barnabè, M., Czoske, O., Koopmans, L. V. E., et al. 2009, Monthly Notices of the Royal Astronomical Society, 399, 21
  * Barnabè et al. (2012) Barnabè, M., Dutton, A. A., Marshall, P. J., et al. 2012, Monthly Notices of the Royal Astronomical Society, 423, 1073
  * Bertin & Lombardi (2006) Bertin, G. & Lombardi, M. 2006, ApJ, 648, L17
  * Binney & Tremaine (1987) Binney, J. & Tremaine, S. 1987, Galactic dynamics
  * Birrer & Amara (2018) Birrer, S. & Amara, A. 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2016) Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmology Astropart. Phys., 8, 020
  * Birrer et al. (2022a) Birrer, S., Dhawan, S., & Shajib, A. J. 2022a, ApJ, 924, 2
  * Birrer et al. (2022b) Birrer, S., Millon, M., Sluse, D., et al. 2022b, Time-Delay Cosmography: Measuring the Hubble Constant and other cosmological parameters with strong gravitational lensing
  * Birrer et al. (2020) Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, A165
  * Birrer et al. (2021) Birrer, S., Shajib, A. J., Gilman, D., et al. 2021, JOSS, 6, 3283
  * 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
  * Blakeslee et al. (2021) Blakeslee, J. P., Jensen, J. B., Ma, C.-P., Milne, P. A., & Greene, J. E. 2021, ApJ, 911, 65
  * Blum et al. (2020) Blum, K., Castorina, E., & Simonović, M. 2020, arXiv e-prints, arXiv:2001.07182
  * Bonvin et al. (2017) Bonvin, V., Courbin, F., Suyu, S. H., et al. 2017, MNRAS, 465, 4914
  * Brout et al. (2022) Brout, D., Scolnic, D., Popovic, B., et al. 2022, The Astrophysical Journal, 938, 110
  * Buckley-Geer et al. (2020) Buckley-Geer, E. J., Lin, H., Rusu, C. E., et al. 2020, MNRAS, 498, 3241
  * Cappellari (2002) Cappellari, M. 2002, MNRAS, 333, 400
  * Cappellari (2008) Cappellari, M. 2008, MNRAS, 390, 71
  * Cappellari (2016) Cappellari, M. 2016, ARA&A, 54, 597
  * Cappellari (2017) Cappellari, M. 2017, MNRAS, 466, 798
  * Cappellari (2020) Cappellari, M. 2020, Monthly Notices of the Royal Astronomical Society, 494, 4819
  * Cappellari (2022) Cappellari, M. 2022, Full spectrum fitting with photometry in ppxf: non-parametric star formation history, metallicity and the quenching boundary from 3200 LEGA-C galaxies at redshift z 0.8
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, Monthly Notices of the Royal Astronomical Society, 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
  * Chang et al. (2013) Chang, Y.-Y., van der Wel, A., Rix, H.-W., et al. 2013, The Astrophysical Journal, 773, 149
  * Chen et al. (2019) Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743
  * Chen et al. (2021a) Chen, G. C.-F., Fassnacht, C. D., Suyu, S. H., et al. 2021a, A&A, 652, A7
  * Chen et al. (2021b) Chen, G. C. F., Treu, T., Fassnacht, C. D., et al. 2021b, Monthly Notices of the Royal Astronomical Society, 508, 755
  * Collett et al. (2018) Collett, T. E., Oldham, L. J., Smith, R. J., et al. 2018, Science, 360, 1342
  * Courbin et al. (2005) Courbin, F., Eigenbrod, A., Vuissoz, C., Meylan, G., & Magain, P. 2005, 225, 297
  * Courteau et al. (2014) Courteau, S., Cappellari, M., de Jong, R. S., et al. 2014, Reviews of Modern Physics, 86, 47
  * de Zeeuw et al. (1996) de Zeeuw, P. T., Evans, N. W., & Schwarzschild, M. 1996, Monthly Notices of the Royal Astronomical Society, 280, 903
  * Di Valentino et al. (2021) Di Valentino, E., Mena, O., Pan, S., et al. 2021, Classical and Quantum Gravity, 38, 153001
  * Efstathiou (2021) Efstathiou, G. 2021, MNRAS, 505, 3866
  * Emsellem et al. (1994) Emsellem, E., Monnet, G., Bacon, R., & Nieto, J.-L. 1994, A&A, 285, 739
  * 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
  * Foxley-Marrable et al. (2018) Foxley-Marrable, M., Collett, T. E., Vernardos, G., Goldstein, D. A., & Bacon, D. 2018, Monthly Notices of the Royal Astronomical Society, 478, 5081
  * Freedman (2021) Freedman, W. L. 2021, ApJ, 919, 16
  * Freedman et al. (2019) Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34
  * Freedman et al. (2020) Freedman, W. L., Madore, B. F., Hoyt, T., et al. 2020, ApJ, 891, 57
  * Fruchter & Hook (2002) Fruchter, A. S. & Hook, R. N. 2002, Publications of the Astronomical Society of the Pacific, 114, 144
  * Gilman et al. (2020) Gilman, D., Birrer, S., & Treu, T. 2020, A&A, 642, A194
  * Gomer et al. (2022) Gomer, M. R., Sluse, D., Van de Vyvere, L., Birrer, S., & Courbin, F. 2022, A&A, 667, A86
  * Gonneau et al. (2020) Gonneau, A., Lyubenova, M., Lançon, A., et al. 2020, A&A, 634, A133
  * Gonzaga et al. (2012) Gonzaga, S., Hack, W., Fruchter, A., & Mack, J. 2012, The DrizzlePac Handbook
  * Goodman & Weare (2010) Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65–80
  * Graham et al. (2018) Graham, M. T., Cappellari, M., Li, H., et al. 2018, Monthly Notices of the Royal Astronomical Society, 477, 4711
  * Greene et al. (2013) Greene, Z. S., Suyu, S. H., Treu, T., et al. 2013, ApJ, 768, 39
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Jeans (1922) Jeans, J. H. 1922, Monthly Notices of the Royal Astronomical Society, 82, 122
  * Jee et al. (2019) Jee, I., Suyu, S. H., Komatsu, E., et al. 2019, Science, 365, 1134
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., & Others. 2001, SciPy: Open source scientific tools for Python
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt (IOS Press BV, Amsterdam, Netherlands), 87 – 90
  * Knox & Millea (2020) Knox, L. & Millea, M. 2020, Phys. Rev. D, 101, 043533
  * Kochanek (2020) Kochanek, C. S. 2020, MNRAS, 493, 1725–1735
  * Kourkchi et al. (2020) Kourkchi, E., Tully, R. B., Eftekharzadeh, S., et al. 2020, ApJ, 902, 145
  * Krajnović et al. (2006) Krajnović, D., Cappellari, M., de Zeeuw, P. T., & Copin, Y. 2006, Monthly Notices of the Royal Astronomical Society, 366, 787
  * Li et al. (2018) Li, H., Mao, S., Cappellari, M., et al. 2018, The Astrophysical Journal, 863, L19
  * Millon et al. (2020a) Millon, M., Courbin, F., Bonvin, V., et al. 2020a, Astronomy and Astrophysics, 640, A105
  * Millon et al. (2020b) Millon, M., Galan, A., Courbin, F., et al. 2020b, A&A, 639, A101
  * More et al. (2017) More, A., Suyu, S. H., Oguri, M., More, S., & Lee, C.-H. 2017, The Astrophysical Journal, 835, L25
  * Morrissey et al. (2012) Morrissey, P., Matuszewski, M., Martin, C., et al. 2012, 8446, 844613
  * Morrissey et al. (2018) Morrissey, P., Matuszewski, M., Martin, D. C., et al. 2018, The Astrophysical Journal, 864, 93
  * 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
  * Oke et al. (1995) Oke, J. B., Cohen, J. G., Carr, M., et al. 1995, Publications of the Astronomical Society of the Pacific, 107, 375
  * Oliphant (2015) Oliphant, T. E. 2015, Guide to NumPy, 2nd edn. (USA: CreateSpace Independent Publishing Platform)
  * Pesce et al. (2020) Pesce, D. W., Braatz, J. A., Reid, M. J., et al. 2020, ApJ, 891, L1
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A6
  * Raftery (1995) Raftery, A. E. 1995, Sociological Methodology, 25, 111
  * Refsdal (1964) Refsdal, S. 1964, MNRAS, 128, 307
  * Riess et al. (2022) Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, The Astrophysical Journal, 934, L7
  * Robertson (2013) Robertson, J. G. 2013, Publications of the Astronomical Society of Australia, 30, e048
  * Rusu et al. (2017) Rusu, C. E., Fassnacht, C. D., Sluse, D., et al. 2017, MNRAS, 467, 4220
  * Rusu et al. (2020) Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440
  * Schneider et al. (1992) Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses
  * Schneider & Sluse (2013) Schneider, P. & Sluse, D. 2013, A&A, 559, A37
  * Schneider & Sluse (2014) Schneider, P. & Sluse, D. 2014, A&A, 564, A103
  * Shajib (2019) Shajib, A. J. 2019, MNRAS, 488, 1387–1400
  * Shajib et al. (2020) Shajib, A. J., Birrer, S., Treu, T., et al. 2020, MNRAS, 494, 6072
  * Shajib et al. (2019) Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649
  * Shajib et al. (2022a) Shajib, A. J., Glazebrook, K., Barone, T., et al. 2022a, LensingETC: a tool to optimize multi-filter imaging campaigns of galaxy-scale strong lensing systems
  * Shajib et al. (2018) Shajib, A. J., Treu, T., & Agnello, A. 2018, MNRAS, 473, 210
  * Shajib et al. (2022b) Shajib, A. J., Vernardos, G., Collett, T. E., et al. 2022b, Strong Lensing by Galaxies
  * Sluse et al. (2007) Sluse, D., Claeskens, J.-F., Hutsemékers, D., & Surdej, J. 2007, Astronomy and Astrophysics, 468, 885
  * Sluse et al. (2003) Sluse, D., Surdej, J., Claeskens, J.-F., et al. 2003, A&A, 406, L43
  * Sonnenfeld et al. (2015) Sonnenfeld, A., Treu, T., Marshall, P. J., et al. 2015, ApJ, 800, 94
  * Suyu et al. (2013) Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70
  * Suyu et al. (2010) Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201
  * Suyu et al. (2014) Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, Monthly Notices of the Royal Astronomical Society, 371, 703
  * Sérsic (1968) Sérsic, J. L. 1968, Atlas de Galaxias Australes
  * Taubenberger et al. (2019) Taubenberger, S., Suyu, S. H., Komatsu, E., et al. 2019, A&A, 628, L7
  * Tewes et al. (2013) Tewes, M., Courbin, F., Meylan, G., et al. 2013, A&A, 556, A22
  * Tihhonova et al. (2018) Tihhonova, O., Courbin, F., Harvey, D., et al. 2018, MNRAS, 477, 5657
  * Treu et al. (2018) Treu, T., Agnello, A., Baumer, M. A., et al. 2018, MNRAS, 481, 1041
  * Treu & Koopmans (2002) Treu, T. & Koopmans, L. V. E. 2002, MNRAS, 337, L6
  * Treu & Marshall (2016) Treu, T. & Marshall, P. J. 2016, A&A Rev., 24, 11
  * Treu et al. (2022) Treu, T., Suyu, S. H., & Marshall, P. J. 2022, Strong lensing time-delay cosmography in the 2020s
  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., Singh, H. P., & Bell, D. J. 2004, The Astrophysical Journal Supplement Series, 152, 251
  * Van de Vyvere et al. (2022a) Van de Vyvere, L., Gomer, M. R., Sluse, D., et al. 2022a, Astronomy and Astrophysics, 659, A127
  * Van de Vyvere et al. (2022b) Van de Vyvere, L., Sluse, D., Gomer, M. R., & Mukherjee, S. 2022b, Astronomy and Astrophysics, 663, A179
  * Verde et al. (2019) Verde, L., Treu, T., & Riess, A. G. 2019, Nature Astronomy, 3, 891–895
  * Waskom et al. (2014) Waskom, M., Botvinnik, O., Hobson, P., et al. 2014, seaborn: v0.5.0 (November 2014)
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, Astronomy and Astrophysics Supplement Series, 143, 9
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420
  * Yahalomi et al. (2017) Yahalomi, D. A., Schechter, P. L., & Wambsganss, J. 2017, A Quadruply Lensed SN Ia: Gaining a Time-Delay…Losing a Standard Candle
  * Yıldırım et al. (2021) Yıldırım, A., Suyu, S. H., Chen, G. C.-F., & Komatsu, E. 2021, arXiv:2109.14615 [astro-ph] [arXiv:2109.14615]
  * Yıldırım et al. (2020) Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, Monthly Notices of the Royal Astronomical Society, 493, 4783