© The Authors 2026

1. Introduction

The so-called Hubble tension – the statistical difference between the Hubble constant, H₀, measured by several late-Universe probes and that inferred from the early-Universe probes (Abdalla et al. 2022; Di Valentino et al. 2025) – is arguably the most pressing issue in observational cosmology. If the tension is confirmed to be real and not due to unknown systematic uncertainties in multiple measurements, solving it would require new physics beyond the standard flat Λ cold dark matter cosmology. For example, one class of solutions would require modifying the expansion history before recombination to adjust the sound horizon at recombination, either by introducing new relativistic particles or early dark energy (Knox & Millea 2020). Other solutions are also possible. However, it is not trivial to find a solution that solves the tension while satisfying all the other cosmological probes (Vagnozzi 2023), even with models featuring a late-time evolution in the dark energy favored by recent measurements (Lynch et al. 2024).

Time-delay cosmography (Treu et al. 2016, 2022; Treu & Shajib 2024; Birrer et al. 2024) measures the Hubble constant and other cosmological parameters in a way that is complementary and fully independent of traditional methods. It is a “local” measurement in the sense that it measures distances out to the lensed sources (typically z ∼ 1 − 3), thus comparable to the local distance ladder, but it reaches the Hubble flow directly in one step. Furthermore, by measuring the absolute time-delay distance, which is a ratio of angular diameter distances, its degeneracies are complementary to those of other probes (Linder 2011). As a cosmological probe, it is most similar to baryonic acoustic oscillations, but it does not require knowledge of the sound horizon at recombination for absolute calibration.

As the culmination of two decades of effort, in 2020 the Time-Delay COSMOgraphy (TDCOSMO) collaboration and its forerunners (COSMOGRAIL, H0LiCOW, SHARP, STRIDES; Eigenbrod et al. 2005; Lagattuta et al. 2012; Suyu et al. 2017; Treu et al. 2018) published a measurement of the Hubble constant with 2% precision based on seven lenses (Millon et al. 2020, hereafter, TDCOSMO-I). The measurement agreed very well with the local distance ladder measurement by the Supernova H₀ for the Equation of State (SH0ES) team (Riess et al. 2022), strengthening the tension with the early-Universe probe.

The TDCOSMO-I measurement achieved its precision by breaking the mass-sheet degeneracy (MSD; Falco et al. 1985) implicitly via the assumption of simply parameterized mass density profiles for the deflector galaxies. These models – power-law mass density profile, or stars plus Navarro–Frenk–White dark matter density profiles (Navarro et al. 1997) – had been shown to describe the mass profiles of elliptical galaxies by studies based on stellar kinematics (Cappellari et al. 2013) and X-ray measurements (Humphrey & Buote 2010).

Given the importance of the tension, the TDCOSMO collaboration decided to repeat the measurements under virtually no assumptions about the mass density profile. By choosing a parameterization of the mass density profile that is maximally degenerate with H₀ via the MSD, Birrer et al. (2020, TDCOSMO-IV) showed that the uncertainties increase from 2% to 9% for the same systems, where the information content of the unresolved stellar kinematics available at that time was the limiting factor. TDCOSMO Collaboration (2025) provide the most recent measurement from time-delay cosmography, increasing the sample size of time-delay lenses to eight and improving, among other aspects, all the stellar kinematic measurements from new datasets, leading to a 4−5% measurement of the Hubble constant ($H_0=71.6_{-3.3}^{+3.9}$ km s⁻¹ Mpc⁻¹ for the TDCOSMO lenses in combination with a prior on Ω_m from Pantheon+ supernovae and $74.8_{-3.4}^{+3.5}$ km s⁻¹ Mpc⁻¹ from TDCOSMO+SLACS+SL2S in combination with a prior on Ω_m from the DESI baryon acoustic oscillation)¹. Stellar kinematics, in particular spatially resolved information, is a crucial element in further improving the precision of time-delay cosmography under these weaker assumptions, as it breaks the MSD and the mass-anisotropy degeneracy (MAD; Binney & Mamon 1982; Treu & Koopmans 2002; Shajib et al. 2018; Birrer et al. 2020; Birrer 2021; Yıldırım et al. 2020, 2023).

Shajib et al. (2023) have presented a spatially resolved velocity dispersion measurement of RXJ1131−1231, the first of its kind. This measurement is based on integral field spectroscopy (IFS) from the Keck Observatory’s Keck Cosmic Web Imager (KCWI) instrument. In combination with lens models and time delays previously published by our collaboration (Tewes et al. 2013; Suyu et al. 2013, 2014), Shajib et al. (2023) measured $H_0=77.1_{-7.1}^{+7.3}$ km s⁻¹ Mpc⁻¹, reaching a precision on H₀ comparable to that achieved by Birrer et al. (2020) with seven lenses, demonstrating the power of spatially resolved kinematics. However, the limited spatial resolution of the seeing-limited ground-based data meant that the MAD still left a significant uncertainty. Simulations have shown that with the spatial resolution of the NIRSpec integral field unit (IFU) on board JWST, one can dramatically improve the precision of the measurement (Yıldırım et al. 2020, 2023; Wang et al. 2025).

We present here the first space-based stellar kinematic map for a time-delay gravitational lens system. The data for RXJ1131−1231 are from the NIRSpec IFU on JWST and represent a sixfold improvement in angular resolution with respect to the Keck data used by Shajib et al. (2023). Furthermore, the data probe the rest frame Ca II triplet region, which is believed to be the most accurate probe of stellar kinematics in elliptical galaxies (Barth et al. 2002). To achieve maximum precision and accuracy from the data, we modeled the three spectral components (lens galaxy, quasar, and quasar host) simultaneously, accounting for the effects of JWST’s point spread function (PSF). We implemented the improved methodology from Knabel et al. (2025a) to account for systematic uncertainties arising from the choice of stellar template libraries, as well as other sources, achieving a 1.58% systematic covariance. A companion paper will present the lensing and dynamical models, along with the related inference of the Hubble constant (Wang et al., in preparation).

The paper is organized as follows. In Sect. 2, we describe the data acquisition and reduction. In Sect. 3, we describe how we obtained additional pieces of information required for the kinematic extraction. In Sect. 4, we describe the measurement in detail. In Sect. 5, we present the measured kinematic maps and uncertainty estimates, and compare them with the previous measurement of resolved kinematics for the same system. In Sect. 6, we conclude the paper with a discussion and comparison with previous measurements.

2. Observations and data reduction

In this section, we first provide a brief description of the lens system RXJ1131−1231 in Sect. 2.1. Then we describe the spectroscopic observation with JWST-NIRSpec in Sect. 2.2 and the data reduction procedure in Sect. 2.3.

2.1. Description of the lens system

The quasar lens system RXJ1131−1231 was first discovered by Sluse et al. (2003). It comprises an elliptical galaxy as the deflector with a redshift of z_d = 0.295, while the quasar’s (i.e., the lensed source) redshift is z_s = 0.654 (Sluse et al. 2007). Due to the relatively low redshift of the deflector, the system appears brighter and larger in angular size than most other such systems. The intricate features within the Einstein ring offer valuable insights for constraining the lens mass model (see Fig. 1).

Fig. 1. Imaging of the system RXJ1131−1231. Left panel: HST-ACS image in the F814W band. The four quasar images are labeled 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. The North and East directions, along with a 1″ scale, are also illustrated. The $3\stackrel{″}{.}2\times 3\stackrel{″}{.}1$ yellow squares represent the dithered 2 × 2 mosaic pattern of the NIRSpec IFU, which covers a $6\stackrel{″}{.}1\times 5\stackrel{″}{.}55$ field of view over the system. Right panel: “White-light” image of the system from the JWST-NIRSpec datacube, summed within 8255−8890 Å in the lens galaxy’s rest frame. The white bar represents 1″, and the red and yellow arrows point to the east and north, respectively. The blue contours show the region within which the spaxels are summed to obtain a high-S/N spectrum of the host galaxy. The black regions around the quasar images were used to extract summed spectra for modeling the quasar’s emission line components. HST and JWST logo credits: NASA.

This system has been extensively studied due to its early discovery and the abundance of lensing information. Time delays for this system were measured by Tewes et al. (2013). Suyu et al. (2013, 2014) conducted lens modeling and cosmographic analysis of this system. These authors employed simply parameterized lens models based on high-resolution imaging from the Hubble Space Telescope’s (HST) Advanced Camera for Surveys (ACS; HST-GO 9744; PI: Kochanek), the measured time delays, single-aperture velocity dispersion, and an external convergence estimate to infer $H_0=80.0_{-4.7}^{+4.5}$ km s⁻¹ Mpc⁻¹. However, as discussed in the introduction, these simply parameterized lens models implicitly break the MSD. The two independent mass profiles used in the analysis by Suyu et al. (2014) nonetheless account for additional systematic uncertainties in H₀, partly due to the MSD. Birrer et al. (2016) independently modeled the mass of this system while accounting for the MSD with a prior on the quasar-host-galaxy size. The findings from these authors revealed that the choice of prior on the dynamical modeling anisotropy dominated the systematic uncertainties in inferring H₀. These obstacles are overcome by the new data presented in this work.

2.2. JWST-NIRSpec spectroscopy

The JWST-NIRSpec IFS for RXJ1131−1231 was obtained through Cycle 1 program JWST-GO-1794 (PI: Suyu, co-PIs: Yıldırım, Treu). The observation was carried out on May 9, 2023. We observed the system with G140M grating and the F100LP filter, which covers the observed wavelength range of 0.97−1.84 μm with a nominal resolution R ∼ 1000. We performed a 2 × 2 mosaic mapping pattern to obtain a total field of view of $6\stackrel{″}{.}1″\times 5\stackrel{″}{.}55$ encompassing the entire lensing system, including the lensed quasar images and the arcs from the quasar host galaxy. For individual tiles in the mosaic, we performed a four-point dithering with sub-arcsecond shifts. We used the “IRS2” readout mode with the “NRSIRS2” readout pattern to achieve smaller noise levels for the outer regions with lower surface brightness of the lensing galaxy. For each of the dither positions, we exposed for 1050.4 s with two integrations and seven groups to prevent saturation of the pixels from the bright quasar images. For each mosaic tile, we also obtained a multi-shutter array (MSA) leakage calibration exposure of 1050.4 s. Additionally, we obtained a background exposure from a nearby empty patch of the sky with an exposure time of 1050.4 s to mitigate potential contamination from stray light.

2.3. Data reduction

We developed a custom data reduction pipeline REGALJUMPER², which is built on the standard JWST data reduction pipeline³. Our custom pipeline includes additional steps to clean up 1/f noise, cosmic rays (CRs), and artifacts, in addition to the standard three-stage data reduction procedure implemented by the standard pipeline. Below, we describe how these additional steps are interleaved within the standard reduction stages.

• In Stage 1, the standard pipeline performs detector-level corrections to generate countrate images (i.e., rate files) per exposure. We included the standard steps: group scale correction, data quality initialization, saturation detection, superbias subtraction, reference pixel correction, linearity correction, dark current subtraction, jump detection with default settings for CR shower flagging, ramp fitting, and gain scale correction. After these steps, we manually identified outliers in the produced files, which are most likely to be unidentified hot pixels, and marked them, along with their adjacent pixels, as bad pixels. We then cleaned the 1/f noise using the software package NSCLEAN (Rauscher 2024). In this step, we used a hand-drawn trace mask⁴ to achieve a robust cleaning, as recommended by Rauscher (2024). We also added “snowball” regions detected previously in this stage to this mask for enhanced robustness. Next, we performed additional CR cleaning using the Python package LACOSMIC (van Dokkum 2001). Next, we subtracted the median of the obtained “leakcal” exposures to manually perform the MSA leakage correction.

• In Stage 2, the standard pipeline performs additional instrument-level and observing-mode corrections to produce flux-calibrated exposure files (i.e., cal files). We included the standard steps: assigning the world coordinate system (WCS), MSA flag-open correction, source type determination, flat-field correction, path-loss correction, and photometric calibration. After these steps, we manually detected and flagged outliers and their adjacent pixels again, which are likely to be unidentified hot pixels in the flat field reference frame.

• In stage 3, the standard pipeline combines the individual exposures to construct the 3D datacube. We used the drizzling method to build the cube with a pixel scale of $0\stackrel{″}{.}05$ and kept the coordinate system aligned with the IFU frame. After that, we manually performed background subtraction, following the procedure recommended by the TEMPLATES team⁵ to subtract a simulated background level using the JWST backgrounds tool (JBT)⁶. Although we did not use the dedicated background exposures to subtract the mean background level, we constructed the background datacube following the same reduction procedure and used it to estimate the background noise level to add to the total noise level of the data during further analysis. The reduced spectra have a reciprocal dispersion of 6.36 Å per pixel.

After constructing the datacube, we correct the “wiggles” in the spectra, which are a manifestation of the resampling noise due to undersampling of the PSF (Law et al. 2023; Perna et al. 2023). We used the software package RACCOON⁷ (Shajib 2025) to clean the wiggles from the spectra. We describe the settings we used for this step in Appendix A and illustrate an example of the wiggle-cleaned spectra.

3. Ancillary ingredients for kinematic extraction

To achieve an accurate absolute calibration of the stellar velocity dispersion, we measured the effective line spread function, as described in Sect. 3.1. We then constructed a lens model of the white-light image obtained from the data cube to properly reconstruct the effective PSF and account for cross-contamination by the multiple spectral components (Sect. 3.2).

3.1. Line spread function

Shajib et al. (2025) provide the wavelength-dependent formula for the instrumental dispersion σ_inst assuming a Gaussian line spread function (LSF) as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{inst}}(\lambda )={[\frac{{{\sigma }_{\mathrm{piv}}^{\prime }}^2}{{\{1+\alpha (\lambda -{\lambda }_{\mathrm{piv}})/[\mathrm{\mu }\mathrm{m}]\}}^2}-{\sigma }_{\mathrm{PN}}^2]}^{1/2},\end{array}$$(1)

where $\sigma {\prime }_{\mathrm{piv}}=117.01\pm 2.25$ km s⁻¹, α = 0.668 ± 0.035, σ_PN = 6.90 ± 0.49 km s⁻¹, and fixed λ_piv = 13539.59 Å. This LSF was derived by fitting narrow emission lines from the planetary nebula SMP LMC 58 observed with the same grating and filter configuration of NIRSpec (Program JWST-CAL-1492, PI: T. Beck). At the observed wavelength of the Ca II triplets (corresponding to rest-frame mean wavelength of 8570 Å), the instrumental dispersions are 140.6 ± 3.0 km s⁻¹ and 113.0 ± 2.0 km s⁻¹ for the lens and the quasar-host galaxies, respectively. The corresponding LSF full widths at half maximum (FWHMs) are 12.16 ± 0.26 Å and 12.48 ± 0.22 Å, respectively.

3.2. Lens modeling of the white-light image

We performed lens modeling of the image shown in Fig. 2, which was obtained by summing the datacube between 8700 and 8800 Å in the lens rest frame. This wavelength range is immediately redward of the Ca II region, and we chose this range to obtain the lens galaxy’s light profile close to the Ca II triplet to deblend the lens galaxy’s light fraction relative to the other components, thereby accurately extracting the lens galaxy’s S/N for use in Voronoi binning. Additionally, the constrained lens galaxy light profile would be used as the tracer distribution in the Jeans dynamical modeling, along with the effective PSF at that wavelength that is reconstructed simultaneously with the lens modeling. From this lens model, we obtained the lens galaxy’s light profile for estimating the signal-to-noise ratio (S/N) at each spaxel. During lens modeling, we also reconstructed the PSF iteratively, which is necessary for dynamical modeling with kinematic maps. We added the background noise level estimated from the background exposure’s datacube to the noise level provided by the reduction pipeline, which is a combination of the propagated detector-level Poisson noise, the read noise, and the flat field’s uncertainty levels.

Fig. 2. Lens modeling with the image obtained by summing the datacube within 8700−8800 Å in the lens galaxy’s rest frame. This wavelength range covers the Ca triplets of the lens galaxy. Top row, from left to right: Illustration of the data, the optimized-model-based reconstruction of the data, and the normalized residuals. Bottom row, from left to right: Illustration of the initial PSF model from STPSF, the iteratively reconstructed PSF, and the fraction of the lens light compared to the total light in the image. The reconstructed PSF will be necessary for dynamical modeling, and we use the lens light fraction to obtain the S/N of the lens galaxy for Voronoi binning.

We performed lens modeling with the publicly available software package LENSTRONOMY (Birrer & Amara 2018; Birrer et al. 2021). For the mass model, we used an elliptical power-law (EPL) model (Tessore & Metcalf 2015) and a residual shear field (also called external shear). We modeled the satellite galaxy’s mass distribution with a singular isothermal sphere with a fixed Einstein radius ${\theta }_{\mathrm{E}}=0\stackrel{″}{.}2$, which is the best-fit value from Suyu et al. (2013). We modeled the light distribution of the lens galaxy as a superposition of two concentric elliptical Sérsic profiles (Sérsic 1968). We also modeled the satellite’s light profile with a circular Sérsic profile, with fixed $R_{\mathrm{eff}}=0\stackrel{″}{.}01$ and n_Sérsic = 1 following Suyu et al. (2013). We modeled the quasar-host galaxy’s light distribution with a basis set including one elliptical Sérsic profile and a set of exponential shapelets with n_max = 20 (Refregier 2003; Birrer et al. 2015; Berge et al. 2019), with their centroids joined. We modeled the quasar point images with the PSF, with their positions freely varied on the image plane.

We obtained an initial PSF estimate using the STPSF software package⁸ (Perrin et al. 2014), which is a mean of the generated PSFs at six uniformly spaced wavelengths within the wavelength range 8700−8800 Å, that was used to construct the white-light image. We then iteratively reconstructed the PSF along with the lens model optimization (Chen et al. 2016; Shajib et al. 2019). The reconstructed image is shown in Fig. 2. We used ASTROPY⁹ (Astropy Collaboration 2022) to fit a 2D elliptical Gaussian profile to the reconstructed PSF. We find agreement with previous results (e.g., D’Eugenio et al. 2024) that the NIRSpec PSF is slightly more elongated along the x-axis (FWHM $0\stackrel{″}{.}1666$) than along the y-axis (FWHM $0\stackrel{″}{.}1308$). The geometric-mean or circularized FWHM of the PSF is $0\stackrel{″}{.}1476$.

4. Extraction of the kinematic maps

In this section, we describe our methodology to extract the resolved stellar kinematics from the JWST-NIRSpec datacube. We generally follow the methodology presented by Knabel et al. (2025a) to account for the systematic uncertainties stemming from the choice of stellar template libraries, in addition to other sources of systematic uncertainties. For this purpose, we developed the SQUIRREL pipeline¹⁰, which is built on the penalized PiXel Fitting (PPXF) software package¹¹ (Cappellari 2017, 2023). We describe several individual steps in our kinematic fitting procedures in Sects. 4.1–4.5 and then summarize the steps involved in the entire fitting procedure in Sect. 4.6.

4.1. Voronoi binning

For appropriate Voronoi binning (Cappellari & Copin 2003), we first set a target S/N for each bin. We define a “specific” S/N (sS/N) as sS/N ≡ (S/N)/$\sqrt{\mathrm{\Delta }L}$, where S is the summed flux within a wavelength range ΔL and N is the noise summed in quadrature within the same range. The $\sqrt{\mathrm{\Delta }L}$ term standardizes the sS/N as it cancels out the improvement in the S/N solely due to an increase in the summed wavelength range. We took the wavelength range 8700−8830 Å in the lens rest frame, slightly redward of the Ca II triplets, to compute the sS/N. We set our target sS/N for each Voronoi bin to 90 Å^−1/2, which we determined through an trial-and-error experiment to be sufficient for providing stable kinematic measurements (i.e., devoid of extremely high or low spurious values that are indicative of amplified systematic effects due to suboptimal S/N)¹².

We estimated the sS/N of the lens galaxy spectrum for each spaxel by considering the fractional contribution of the lens galaxy’s light, as determined by the best-fit model described in Sect. 3.2. This involved scaling the total sS/N, which includes signals from all components, by the lens galaxy’s light fraction at each spaxel (see Fig. 2). For Voronoi binning, we included only those spaxels that are within $1\stackrel{″}{.}5$ of the lens galaxy center and have more than half of the total flux contributed by the lens galaxy. We also manually exclude from the initial Voronoi binning an aperture around the satellite galaxy, where the satellite’s light is more than 5% of the lens galaxy’s light. We used the software package¹³ VORBIN to perform Voronoi-binning with the minimum target sS/N per bin. We then manually combined a few single-spaxel bins near the center so that a Voronoi bin includes at least two spaxels and a few Voronoi bins in the outer periphery for which the initial Voronoi binning did not lead to a total sS/N ≥ 90 Å^−1/2. Finally, we included back the aperture containing the satellite galaxy as the last bin, yielding a total of 32 Voronoi bins. The Voronoi binning map and the sS/N in each bin are illustrated in Fig. 3.

Fig. 3. Voronoi binning map. Left panel: Voronoi bins outlined over the white-light image of the system. Middle panel: Bin numbering. Right panel: Total sS/N for each Voronoi bin. The horizontal dashed line marks 90 Å−1/2, which we set as the minimum target for each bin. Only the last Voronoi bin, which contains the satellite galaxy, is slightly below this sS/N target.

4.2. Log-rebinning of spectra

We rebinned the spectral data, which is linearly sampled in wavelength, into a logarithmically sampled one using PPXF. Such log-rebinning is a standard procedure for kinematic measurements, facilitating numerical fitting with a constant velocity scale per wavelength-pixel (Cappellari 2017). We chose the default rebinning settings of PPXF that maintain the same pixel numbers before and after rebinning, resulting in a velocity scale of 171.92 km s⁻¹. We also created a covariance matrix for the log-rebinned spectra, as log-rebinning induces covariances within the spectra, with stronger covariances between more closely adjacent bins. We created the covariance matrix using the Monte Carlo method, by generating 5000 realizations of the spectra that are statistically consistent with the data within its noise levels and then taking the covariance of these samples after log-rebinning them. This covariance enabled us to fully propagate the uncertainty provided by the reduction pipeline into the final analysis by providing it to the PPXF kinematic fitting routine.

4.3. Spectral templates

We chose two stellar template libraries: the Indo-US (Valdes et al. 2004) and the X-shooter Spectral Library (XSL) DR3 (Verro et al. 2022). We used the “cleaned” libraries provided by Knabel et al. (2025a)¹⁴. We did not use the MILES library, which was also cleaned by Knabel et al. (2025a), since the original MILES library does not cover wavelengths up to the Ca II triplet. Although the Ca II triplet library from Cenarro et al. (2001) is referred to as the MILES-CaT library, because it contains a subset of the MILES stars, we did not use this library as it was not cleaned by Knabel et al. (2025a).

From these libraries, we further removed spectra that have gaps within the fit wavelength ranges around the lens and quasar-host galaxy’s Ca II triplet regions, either by using the flags provided Knabel et al. (2025a) or by directly checking for such gaps within the relevant wavelength range. Ultimately, we had 609 stellar templates in the Indo-US library and 489 stellar templates in the XSL. The Indo-US library has a constant resolution FWHM of 1.35 Å (Beifiori et al. 2011), which corresponds to instrumental dispersion (standard deviation) σ_temp = 18 km s⁻¹. The XSL has an instrumental dispersion σ_temp = 11 km s⁻¹ (for the VIS arm), corresponding to a resolution FWHM of 0.74 Å at the wavelength of the Ca II triplet. The stellar template libraries are broadened by an additional (σ_inst² − σ_temp²)^1/2 before being used in the kinematic fitting to account for the instrumental broadening in the observed absorption line widths.

The Ca II triplet region of the lens galaxy coincides with the Hα, [N II], and [S II] lines from the quasar and the host galaxy. Therefore, we fit the spectra with multiple components (similar to Turner et al. 2024): the lens galaxy’s stellar continuum and absorption features, the host galaxy’s stellar continuum and gas lines, and the quasar emission lines. Below, we describe how we built a set of optimal templates for each of these components.

4.3.1. Template for the host galaxy’s stellar continuum

To create a template for the quasar-host galaxy’s stellar continuum, we first extracted a summed spectrum within 8200−8850 Å in the host galaxy’s rest frame from a region containing the lensed arcs (shown with blue contours in the right panel of Fig. 1). Then, the spectra were fit with the Indo-US and XSL template libraries, while also accounting for the contribution from the lens galaxy’s stellar continuum at the corresponding wavelength range. This capability of PPXF to fit multiple components with different kinematics is crucial for modeling the complex spectrum analyzed in this paper, where emission lines from both the lens galaxy and the quasar plus host (at significantly different redshifts) can have multiple kinematic components (e.g., broad or narrow profiles or different velocities). We followed the methodology of Knabel et al. (2025a) to select the combination of an additive polynomial degree of four and a multiplicative polynomial degree of one that yields a subpercent difference in the best-fit velocity dispersion between the Indo-US library and XSL, while orders any higher do not significantly affect the results. Fig. 4 shows the best-fit model to the host galaxy’s Ca II triplet.

Fig. 4. Fitting of the Ca II region of the quasar host galaxy. The gray bars represent the data points, with the width of each bar corresponding to the pixel size and the height representing the original ±1σ noise level. The gray lines attached to the bars represent the boosted uncertainty levels to achieve χred2 = 1. The red line illustrates the best-fit model. The Ca II triplet features are marked with blue lines.

Using the above fit, we constructed a single template for the host galaxy’s stars by combining the Indo-US library with the corresponding best-fit weights. The wavelength range around the lens galaxy’s Ca II triplet which is fit later in this analysis, corresponds to 6463−6960 Å in the quasar-host galaxy’s rest frame. The Indo-US stellar library, and hence the single template constructed above, has coverage down to this wavelength range. The quasar-host galaxy does not have any prominent kinematic lines around the lens galaxy’s Ca II triplet, and thus any wavelength-dependent difference between the quasar-host galaxy’s Ca II triplet wavelength and that corresponding to the range around the lens galaxy’s Ca II triplet would have a negligible impact on our analysis. Furthermore, although we fit the host galaxy spectra summed from a region primarily containing the disk, rather than the bulge, which may have different kinematic and stellar population properties, this difference is unlikely to affect our analysis for the same reasons mentioned above.

4.3.2. Templates for the quasar spectra

To build the templates for the quasar images, we extracted summed spectra from circular apertures of $0\stackrel{″}{.}2$ radius around the four quasar images (illustrated in Fig. 1). To each of these spectra, we fit multiple emission lines: narrow [N II] and [S II] doublets, and Hα with both narrow and broad components. While each narrow line was fit with a single Gaussian profile with a common velocity shift and dispersion for all of them, the broad Hα line was fit with two Gaussians with independent velocity shifts and dispersions, following Guo et al. (2018)¹⁵. Here, we aimed to construct the templates for the quasar images’ emission lines independently of the physical interpretation of the components. We used an additive polynomial of order four and a multiplicative polynomial of order one to account for the quasar’s continuum. Additionally, the host galaxy’s combined stellar template from Sect. 4.3.1 was also incorporated as an additional component in the fit, with its associated best-fit velocity and dispersion kept fixed as obtained in that section. We combined the best-fit emission lines for a given quasar image to construct its associated template (Fig. 5). In our fitting of the lens galaxy spectra, we used all four templates as independent components because they have different profile shapes, likely due to microlensing (e.g., Hutsemékers et al. 2024), and a given spaxel can receive contributions at varying levels from the four quasar images carried through the PSF spikes. The constructed quasar templates excluded the quasar continuum because it would be accounted for by the polynomials adopted in the kinematic fitting of the lens galaxy’s spectra. We illustrate the expected and effective levels of quasar contribution in each Voronoi-binned spectrum in Appendix B.

Fig. 5. Fitting of the quasar lines at images A, B, C, and D, respectively, from top to bottom. The observed spectra are shown with gray bars, with the width corresponding to the bin width and the height representing the ±1σ noise range. The red lines show the best fit. The individual narrow and emission lines from the best fit are shown in purple, and the lens galaxy’s stellar continuum is shown in emerald. We adopted the best-fit gas lines with fixed relative amplitudes for all three images as free linear components in our kinematic fitting of the spectra across the entire field of view.

4.3.3. Templates for the host galaxy’s emission lines

We also construct a set of emission line templates for the narrow Hα, [N II], and [S II] lines from the host galaxy, which are independent of the quasar’s narrow lines described above in Sect. 4.3.2. If these emission line components were allowed to have individually free amplitudes, they would sometimes result in unrealistic line ratios, especially in cases where their fractional contribution to the spectrum is relatively small. Therefore, we built eight emission line templates where the line ratios between them were fixed within each template. The number eight was adopted because the NIRSpec PSF is expected to have eight prominent spikes; thus, a given spaxel can receive contributions from eight different positions of the lensed arc in a nontrivial manner. We divided the regions traced with blue contours in Fig. 1 into eight segments and fit the spectra summed within each of these segments. In addition to the emission lines with freely varying amplitudes, we included all the other components in these fits, namely the quasar templates, the host galaxy’s stellar continuum template, and the Indo-US template library for the lens galaxy’s spectra. We then constructed the emission line templates from each of the fits by fixing the relative line ratios to the best-fit values. When fitting the lens galaxy spectra in subsequent steps, including these eight emission line templates, each template was allowed to have independent kinematic properties (both in mean velocity and dispersion) because kinematic properties from different points of the quasar-host galaxy can mix in a single spaxel due to nontrivial superposition of the PSF spikes, as mentioned above. The mean velocity of the host galaxy’s Ca II triplet lines obtained in Sect. 4.3.1 is 293 ± 7 km s⁻¹, and the mean velocity of the emission lines obtained in this section is 217.0 ± 0.4 km s⁻¹. This difference suggests potentially distinct kinematic properties of the stars and gas that generated the fit absorption and emission features, respectively. Interpreting these velocities as residual systemic velocity would lead to adjusting the redshift of the host galaxy (and its quasar) by approximately +0.001 to z_s ≈ 0.655.

4.3.4. Template for the lens galaxy’s spectra

We constructed the stellar template for the lens galaxy in two regions: one for the Voronoi bins with luminosity-weighted centers within a $0\stackrel{″}{.}5$ radius from the lens galaxy center, and the other for those outside (for bins 10, 14, 19, and 31). We excluded bin 32, which contains the satellite galaxy, from this step in the template-making process. For a given set of settings on the polynomial orders and the fit wavelength range, the summed spectra from these two regions are fit with a given template library from the Indo-US and XSL, and the template library was combined using the best-fit weights to produce one single template associated with the given settings and template library. We used this combined template when fitting to the spectra for all the Voronoi bins within the associated region to mitigate spurious scatter in the fit kinematics, following standard procedure in the literature (Shajib et al. 2023; Knabel et al. 2025b).

4.4. Fitting spurious spikes

Some spurious spike-like features in the spectra are not well-fit by the stellar templates, and they also do not align with any expected emission lines at the redshifts of the lens or quasar-host galaxy. We explicitly modeled these spurious spikes as emission lines with instrumental resolution, since otherwise, they can have a significant impact on the measured kinematics. These spurious spikes may have astrophysical origins, such as emission lines from faint sources at other redshifts with undetected continua or the Zodiac, or may be data-level defects, such as artifacts from unidentified warm pixels or uncleaned CRs. For a given spectrum and a set of fitting ingredients and settings, we first identified the 5σ outliers or spurious features. Then, we refit the spectra with the same ingredients but now with the identified spikes modeled. Next, we identified 3σ outlier pixels and refit the spectra again with all the identified spikes modeled. Finally, we identified 2.7σ outlier pixels. We performed this identification process in multiple steps until we reached our outlier detection threshold of 2.7σ for a more robust detection of the true outlier at our target threshold. Furthermore, we performed the same detection procedure for a given Voronoi-binned spectrum using both template libraries and selected only the spikes identified in both cases. We set our outlier detection threshold at 2.7σ because the statistically expected number of pixels with absolute residuals larger than 2.7σ is less than one, given the total number of pixels within the fit wavelength range. Although most of the identified spikes have positive amplitudes, a fraction of them have negative amplitudes, and we allowed for this in our model.

4.5. Estimating systematic uncertainty

We accounted for the systematic uncertainties stemming from several aspects, where a particular choice needs to be made for the kinematic fitting, namely the template stellar library, the orders of the additive and multiplicative polynomials used to model the continuum, the detection threshold for spurious spikes, and the fit wavelength range. We adopted multiple viable choices for each of these aspects to marginalize over them and account for the associated systematic covariances. We took the following choices for these aspects:

• Template library: Indo-US and XSL;

• Order of additive and multiplicative polynomials: (2, 2), (3, 2), (4, 2), (2, 3);

• Spurious spike detection threshold: 2.7σ and 2.6σ; and

• Fit wavelength range: 8255–8890 Å, 8280–8890 Å, and 8255–8865 Å.

By combining these choices, we had a total of 48 model setups. For each choice, we regenerated combined templates for the two regions described in Sect. 4.3.4 to fit the Voronoi-binned spectrum. We first combine the 16 model setups within each given fit wavelength range choice using the Bayesian information criterion (BIC) weighting following Knabel et al. (2025a). Then, we combined the three sets of kinematic measurements and their covariances (mean velocities and velocity dispersions) with equal weights. We did not use BIC weighting in this latter case, as it would weigh more heavily the case with a smaller number of data points, which is undesirable.

4.6. Kinematic fitting procedure

Here, we summarize the steps involved in our fitting procedure, which combines all the ingredients and individual steps described above.

1. We performed Voronoi binning as described in Sect. 4.1.

2. We log-rebinned the Voronoi-binned spectra and produced the covariance matrix of the rebinned spectra propagated from the detector-level uncertainties as described in Sect. 4.2. We used the full covariance matrix as input to PPXF, a first for kinematic measurements in the literature.

3. We prepared the stellar templates for both the lens and quasar-host galaxies, the quasar templates, and the emission line templates (associated with the host galaxy) described in Sect. 4.3. For each of the 48 model combinations described in Sect. 4.5, we generated the combined templates within two annuli: one that includes all the Voronoi bins with luminosity-weighted centers within $0\stackrel{″}{.}5$ of the center, and the other including the rest (but excluding bin 32 since it is contaminated by the satellite galaxy).

4. We first boost the covariance of each Voronoi-binned spectrum within the range 8441−8570 Å, so that the median S/N in this range matches that of the rest. As this wavelength range contains the Hα and [N II] lines from the quasar and its host, which may have a very high S/N in the spaxels they appear brightly in, we adjust the S/N so as not to allow any strong residual in this region to bias the overall fit. We fit this spectrum using the prepared templates and then boosted the covariance of the spectrum by a factor of χ² before refitting it to obtain the final kinematic measurements and uncertainties. We treat the formal error on the kinematic measurements provided by PPXF from this final fit as the associated statistical uncertainty. The kinematic values of the host galaxy’s stellar component were fixed to the ones obtained from the fit for the host galaxy’s Ca II triplet in Sect. 4.3.1. However, the narrow emission line sets were allowed to have free mean velocities and dispersions, independently for each set. As a result, the fit included a total of nine kinematic components: the lens galaxy and eight sets of narrow emission line templates.

5. We combined the models using BIC weighting separately within the three sets with different wavelength ranges, and then combined the resultant three sets of values with equal weighting, as described in Sect. 4.5.

Fig. 6 illustrates the kinematic fitting of bin 22, as an example. The fits for the rest of the bins are illustrated in Appendix C.

Fig. 6. Best fit (red line) of the Voronoi-binned spectra in bin 22 as an example. The gray bars represent the observed data, with the width indicating the bin width and the height representing the original 1σ noise level. The vertical gray lines attached to the bars represent the boosted 1σ noise level. These 1σ noise levels for each pixel are estimated from random sampling given the covariance matrix. The Ca II triplet wavelengths from the lens galaxy are marked with orange lines, and the wavelengths of the host galaxy’s emission lines are marked with blue lines. The blue-shaded region marks where we boosted the noise levels around the Hα and [N II] lines. The emerald line shows the contribution from the lens galaxy’s stars. The sets of narrow emission lines from the quasar host galaxy (including spurious spikes identified close to 8300, 8700, and 8900 Å in this case) are shown in purple.

5. Kinematic maps

In this section, we provide our measured kinematic maps in Sect. 5.1 and compare them with previous resolved measurements from ground-based data in Sect. 5.2.

5.1. Kinematic maps from JWST NIRSpec data

The 2D maps of the velocity dispersion and mean velocity measurements in the 32 bins are illustrated in Fig. 7. We show the total covariance matrix of the velocity dispersion measurements in Fig. 8. The mean for absolute values of the off-diagonal elements in the “fractional” covariance matrix [i.e., with the ij-th element being ${\{\mathrm{cov}({\sigma }_{\mathrm{los},i},{\sigma }_{\mathrm{los},j})/{\sigma }_{\mathrm{los},i}{\sigma }_{\mathrm{los},j}\}}^{1/2}$] is 1.24 ± 0.04%. Although Yıldırım et al. (2020) forecast a 4% precision on the H₀ measurement based on simulated, JWST-quality 2D kinematics of RXJ1131−1231, the actual S/N in our data is less than that assumed in their simulations. Therefore, this 4% precision should be regarded as a lower limit, with the actualized uncertainty on H₀ to be presented by Wang et al. (in preparation).

Fig. 7. Velocity dispersion (top row) and mean velocity (bottom row) measurements in 32 Voronoi bins. The first column shows the measurement for each case, and the second column shows the measurement uncertainties. The bin containing the satellite galaxy (i.e., bin 32) is outlined in white. The bins closer to the center have smaller uncertainties, thanks to less contamination from the quasar and its host galaxy. The systemic velocity of 213 km s−1 was subtracted from the mean velocity map. The mean velocity map indicates that the lens galaxy is a slow rotator, with small mean velocities (≲50 km s−1) and a lack of a bipolar velocity profile.

Fig. 8. “Fractional” covariance matrix. In this matrix, the element with indices i, j is [cov(σlos, i, σlos, j)/σlos, iσlos, j]1/2. The diagonal elements represent the combined statistical and systematic uncertainties. The mean for the absolute values of the off-diagonal elements is 1.14 ± 0.04%.

We find that the lens galaxy is a slow rotator, with mean velocities ≲50 km s⁻¹ across the kinematic map. This result agrees with the previous, resolved measurements from the KCWI by Shajib et al. (2023).

In the kinematic map, the bin containing the satellite galaxy (outlined with a white border in Fig. 7) is noticeable for having a slightly different mean velocity and dispersion compared to the surrounding bins. The satellite galaxy was previously attributed as the “satellite” of the main deflector due to its similarity in color (Claeskens et al. 2006). Indeed, the satellite galaxy does not exhibit distinctly identifiable spectral features at a separate redshift, suggesting that its redshift is very close to that of the central galaxy. The mean velocity of the bin after subtracting the systemic velocity of the central deflector is −59 ± 23 km s⁻¹. Although this velocity measurement is a weighted combination of both the satellite and main deflector contributions in bin 32, if taken at face value, the −5σ value of −174 km s⁻¹ would only change the satellite galaxy’s redshift by 5.8 × 10⁻⁴, which is below the typical uncertainty levels (i.e., ∼10⁻³) in spectroscopic redshift measurements.

5.2. Comparison with previous measurement

In this section, we compare resolved kinematic measurements from the JWST-NIRSpec IFS with those from the previous study using the Keck KCWI IFS (Shajib et al. 2023). We used the updated values of the KCWI measurements from TDCOSMO Collaboration (2025), which employed the improved methodology of Knabel et al. (2025a). We took the 1D v_rms profile as the measurements to compare, where v_rms ≡ (σ_los² + v_los²)^1/2 with the v_los being the mean line-of-sight velocity after subtracting off the systemic velocity. We created the 1D v_rms profile in annular bins from the 2D kinematics of the NIRSpec data, following the same methodology as Shajib et al. (2023), but with slightly different annulus widths that are suited to the different Voronoi binning map used in this paper (see Fig. 9).

Fig. 9. One-dimensional profile of the rms velocity vrms ≡ (σlos2 + vlos2)1/2, where vlos is the line-of-sight mean velocity after subtracting the mean systemic velocity of 213 km s−1. Left panel: One-dimensional vrms profile for JWST-NIRSpec measurements from this paper. Middle panel: One-dimensional vrms profile for the Keck-KCWI measurements from Shajib et al. (2023). The vertical error bar represents the total uncertainties, including both statistical and systematic contributions, and the horizontal error bars denote the annulus widths. The dashed lines represent the best-fit model predictions from JWST-only data (red), KCWI-only data (blue), and from a joint fit (black), and the associated shaded regions illustrate the 1σ extent of the predictions from the model posterior. The horizontal arrows represent the FWHMs of the corresponding PSF in each panel. Right panel: One-dimensional light profile from double Sérsic fits performed on the NIRSpec white-light image (red line) as part of the lens modeling in this paper, and from that, on the HST imaging in the F814W filter (blue line) from Shajib et al. (2023). The vertical dashed line marks the $1\stackrel{″}{.}3$ transition radius, where we stitch the light profiles with IR in the inner region and the visible one in the outer region that is used in the dynamical model of the NIRSpec 1D profile. From a likelihood-ratio test, we find that the JWST and KCWI measurements are consistent at ∼1.2σ (i.e., the likelihood-ratio test statistic’s p-value is ∼0.22).

In addition to the two datasets having different PSF FWHMs, the tracer populations are slightly different, given that the measurements were taken at different wavelengths. Therefore, we performed the comparison between the two datasets in the model parameter space after performing dynamical modeling of the 1D v_rms profiles. For the tracer distribution of the KCWI measurements, we used the double Sérsic light profile from Shajib et al. (2023), which was fit to a large cutout of the HST F814W imaging after subtracting the lensed quasars and arcs using the lens model from Suyu et al. (2013). In the IR band, imaging data with a larger field of view than that in the NIRSpec data was not available. Therefore, we took the best-fit light profile from the lens modeling performed on the NIRSpec white-light image in Sect. 3.2 and stitched it with the light profile from the HST imaging to cover the outer region. We took $1\stackrel{″}{.}3$ as the stitching radius since that is the largest extent covered by the 1D v_rms profile from the NIRSpec data.

We perform the dynamical modeling of the 1D v_rms profile with NIRSpec-only, KCWI-only, and the joint data. We took a spherical power-law mass model and performed Jeans modeling of the kinematics. We took a spherical model for simplicity, given that the kinematic profile to be modeled is 1D. We also added a central super-massive black hole in our mass model, with its mass M_bh ∝ θ_E, bh² a free parameter (see Wang et al. 2025, for its impact on 2D kinematics). Our model has five free parameters: the Einstein radius θ_E, the power-law mass profile’s logarithmic slope γ, the super-massive black hole’s Einstein radius θ_E, bh, the ratio of the tangential and radial velocity dispersion σ_t/σ_r ≡ (1 − β_ani)^1/2 with β_ani being the amplitude of a spatially constant anisotropy profile, and the mass-sheet transformation (MST) parameter λ_mst (for a brief explanation of the MST, see Treu & Shajib 2024). Here, the λ_mst parameter encapsulates both the contribution from the external convergence (i.e., κ_ext) and the internal MST. We imposed a prior on the total Einstein radius as ${({\theta }_{\mathrm{E}}^2+{\theta }_{\mathrm{E},\mathrm{bh}}^2)}^{1/2}\sim \mathcal{N}(1\stackrel{″}{.}640,0\stackrel{″}{.}015)$ based on the lens model of Suyu et al. (2013), but adjusting the mean value given our inclusion of the central black hole with an associated Einstein radius of $0\stackrel{″}{.}15$. We also imposed a uniform prior on the anisotropy parameter σ_t/σ_r ∼ 𝒰 (0.72, 1.22), which imposes the bounds −0.5 < β_ani < 0.5. We performed the dynamical modeling with LENSTRONOMY’s GALKIN module and sampled from the posterior using the affine-invariant ensemble sampler EMCEE (Foreman-Mackey et al. 2013). Then, we performed a likelihood ratio test using the statistic Λ ≡ 2(logℒ_comb − logℒ_NIRSpec − logℒ_KCWI)≈7.0, which is expected to follow a χ² distribution with five degrees of freedom. The p-value of this Λ-statistic is 0.22 (z-score ∼1.2σ). Therefore, we cannot exclude the null hypothesis that the two datasets are consistent. It can be noticed in Fig. 10 that the NIRSpec measurements provide more constraint on the stellar anisotropy (σ_t/σ_r), breaking the MAD more effectively than the previous KCWI measurements with seeing-limited resolution, which left this parameter largely unconstrained within the prior volume.

Fig. 10. Comparison of the parameter posteriors of the dynamical modeling performed with JWST-NIRSpec-only data (red), Keck-KCWI-only data (blue), and both combined (black). The four free parameters in our dynamical model are the ratio of the tangential and radial velocity dispersion σt/σr ≡ (1 − βani)1/2 with βani being the amplitude of a spatially constant anisotropy profile, the MST parameter λmst, the main deflector’s Einstein radius θE, the power-law mass profile’s logarithmic slope γ, and the super-massive black hole’s Einstein radius θE, bh. We have imposed a prior on the total Einstein radius (θE2 + θE, bh2)1/2 based on the lens model posterior of Suyu et al. (2013), which dominates the θE posterior illustrated here. All of them are consistent at ≲1σ confidence levels. The JWST NIRSpec’s superior spatial resolution also allows for the exclusion of a large fraction of the prior volume in the anisotropy parameter (σt/σr) that is largely unconstrained by the Keck-KCWI measurements.

6. Conclusion

In this paper, we extracted the 2D kinematic map for the lensed quasar system RXJ1131−1231 from the JWST NIRSpec IFS, which is a first for a time-delay lens galaxy. The 2D resolved kinematics will improve the Hubble constant measurement by breaking the MSD and the MAD, after combining them with the measured time delays, lens models, and estimation of the line-of-sight convergence. This combination and Hubble constant measurement will be presented in a future study (Wang et al., in preparation). Resolved 1D kinematics from the same NIRSpec dataset presented here have been used in the recent measurement of H₀ by the TDCOSMO collaboration (TDCOSMO Collaboration 2025).

The JWST NIRSpec dataset analyzed in this paper is the first of its kind, presenting novel challenges in robustly extracting stellar kinematics from this dataset. The unprecedentedly high S/N of this dataset at the corresponding wavelengths necessitated meticulous modeling of all components, namely the lens galaxy, the background quasar, and its host galaxy, to robustly extract the stellar kinematics of the lens galaxy. Our methodology for this joint modeling has been packaged for future use into the publicly available SQUIRREL pipeline, built upon the PPXF software program. Additionally, the JWST dataset required special processing beyond the standard JWST pipeline to mitigate several sources of non-random noise and artifacts, including uncleaned residual CRs and resampling noise, which manifests as “wiggles” in the spectra. Our JWST data reduction pipeline, REGALJUMPER, is also made publicly available. These methodological frameworks will enable rapid analysis of current and future JWST datasets of time-delay lenses to achieve less than 2% precision in the Hubble constant measurement in the near future (Birrer & Treu 2021).

Acknowledgments

We thank Elizabeth Buckley-Geer and Cheryl Pavlovski for useful discussions that improved this work. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1794. The specific observations analyzed can be accessed via https://dx.doi.org/10.17909/wjva-0750. Code scripts and kinematic maps are available from the corresponding author on request. AJS received support from NASA through STScI grants JWST-GO-2974 and HST-GO-16773. Support for this work was also 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 acknowledges support by NSF through grants NSF-AST-1836016, NSF-AST-1906976, and NSF-AST-2407277, and from the Moore Foundation through grant 8548. SHS thanks the Max Planck Society for support through the Max Planck Fellowship. SHS and HW are 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 the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No. 787886; LENSNOVA: grant agreement No. 771776). CDF acknowledges support for this work from the National Science Foundation under Grant No. AST-2407278. This research made use of PPXF (Cappellari 2017, 2023), PAFIT (Krajnović et al. 2006), VORBIN (Cappellari & Copin 2003), LENSTRONOMY (Birrer & Amara 2018; Birrer & Treu 2021), STPSF (Perrin et al. 2014), 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

Appendix A: Correction for wiggles in spectra

In this appendix, we briefly describe the "wiggle" correction we performed on the reduced datacube spectra and provide the relevant settings. We used the software package RACCOON, the detailed mechanism and features of which are presented by Shajib (2025).

The wiggles or modulation on the spectra are artifacts due to resampling noise, as the pixels are not optimally sampled by the PSF, and our four-point dither pattern is insufficient to recover the optimal sampling (Law et al. 2023). In RACCOON, the wiggle is modeled as a sinusoidal chirp function:

$$\begin{array}{c}\hfill W(\lambda )=1+A(\lambda )[sin({\varphi }_{\lambda })+k_1{sin}^2({\varphi }_{\lambda })+k_{2}sin(3{\varphi }_{\lambda })],\end{array}$$(A.1)

where A(λ) is the wavelength-dependent amplitude and ϕ_λ = λf(λ)+ϕ₀ is wavelength-depedent phase term. The wiggle can be made asymmetric with the k₁ ≠ 0 and the peaks can be symmetrically sharpened (or, de-sharpened) with k₂ ≠ 0. In our use case, we set both k₁ = 0 and k₂ = 0.

The wiggles affect the single-spaxel spectra. However, they wash out if the spectra are summed within an aperture with a few spaxel radius. From our tests, we find that the wiggle signal is almost completely washed out with an aperture radius larger than three pixels, which approximately equals the FWHM of the PSF. We obtain the best-fit wiggle model by fitting the single-spaxel spectra d(λ) with a model

$$\begin{array}{c}\hfill m(\lambda )=W(\lambda )T(\lambda ),\end{array}$$(A.2)

with the template T(λ) for the original spectra constructed as

$$\begin{array}{c}\hfill T(\lambda )=c_{1}a(\lambda )+c_{2}s(\lambda )+c_3{\lambda }^b+c_{4}P(\lambda ),\end{array}$$(A.3)

where c₁, …, c₄ are linear coefficients, a(λ) is the aperture-summed spectra around the spaxel, s(λ) is the spectra summed within a shell centered on the spaxel, and P(λ) is a polynomial of chosen degree. The power-law term λ^b and the polynomial P(λ) are expected to model the change in the continuum level between the single-spaxel spectra and the aperture-summed one. The combination c₁ a(λ)+c₂ s(λ) allows accounting for changes in the emission or absorption line shapes between the modeled single-spaxel spectra and the aperture-summed one (Dumont et al. 2025). In our case, we chose four spaxels as the radius of the aperture to extract a(λ), the inner and outer radii of the shell to extract s(λ) to be three and five spaxels, respectively, and the polynomial degree to be two. The functions A(λ) and ϕ_λ were modeled with B-splines with three and seven knots (thus with five and nine free parameters), respectively. We propagate the uncertainty of the above fitting procedure of the wiggle model and add it to the noise levels of the corrected spectra d(λ)/W(λ).

We perform the cleaning for the spaxels within an aperture of $0\stackrel{″}{.}7$ radius centered on the lens galaxy, and for which the wiggle is detected with more than 3σ and the wiggle model explains more than half of the scatter observed in the “wiggle signal” [i.e., d(λ)/T(λ)]. We chose the $0\stackrel{″}{.}7$ radius, as this would be an a priori safe choice for a radius outside of which all the spaxels would belong to Voronoi bins large enough to smooth out the wiggles after summing within. Although we fit within a narrow wavelength range out of the whole observed region, we fit and correct for the wiggles across the whole wavelength range so that the modeled wiggle signal is less affected by local outliers, noise, or scatter within the fit wavelength region, and the fit wiggle pattern is dictated more by a global pattern. In Fig. A.1, we illustrate an example of the wiggle signal and the best-fit wiggle model. In Fig. A.2, we show the spaxels for which wiggle cleaning was performed. In Fig. A.3, we illustrate the impact of the cleaning on the measured kinematics (here, only the systematics from the template library choice were marginalized). Without cleaning, the measured velocity dispersions would shift down by ∼0.6% on average in the central region of the lens galaxy (that excludes the four outermost bins and the one containing the satellite), and the dispersion in the bin with the satellite galaxy would shift up by ∼10%.

Fig. A.1. Fit wiggle model used to correct for the resampling noise. The data points show the signal of wiggle d(λ)/T(λ) for the central spaxel on the lens galaxy, where T(λ) is obtained from the best-fit model m(λ). The orange line shows the W(λ) function in the best-fit m(λ), with the orange shaded region illustrating the 1σ confidence region. The vertical gray dashed lines show the wavelength range containing the Ca II triplets of the lens galaxy that we fit to measure the lens galaxy’s kinematics. The gray data points are rejected outliers using the false discovery rate method.

Fig. A.2. Spaxels cleaned for wiggles (within the regions enclosed with white borders). Only the spaxels within a circle of $0\stackrel{″}{.}7$ radius from the lens galaxy center were considered to be cleaned. Then, we only cleaned the spaxels for which the wiggle was detected at a confidence level of more than 3σ and the best-fit wiggle model explains more than half of the variance or scatter observed in the wiggle signal. The wiggles are more prominent in spaxels with more flux. As a result, the white contours above primarily encompass the central region of the lens galaxy and the satellite galaxy.

Fig. A.3. Impact of wiggle cleaning on the measured kinematics. The illustrated test is for a single model setup out of the 48 combinations described in Sect. 4.5. The first column displays the fractional differences in the measured dispersions (top row) and mean velocities (bottom row) within the Voronoi bins, while the second column shows these differences normalized by the statistical uncertainty.

Appendix B: Quasar contamination levels

In this appendix, we show the expected quasar contamination in each Voronoi-binned spectrum alongside the effective levels recovered from our PPXF fits (Fig. B.1). The expected levels are derived from the decomposition based on our best-fit lens model in Sect. 3.2. The effective contamination measured by PPXF closely matches these expectations and follows the same spatial trend – minimal in the center and increasing toward the quasar image – demonstrating the effectiveness of our multi-component fitting approach.

Fig. B.1. Expected quasar contamination in each Voronoi-binned spectrum from the best-fit lens model (left) and the effective levels from our PPXF fits (right). The effective contamination is consistent across bins and follows the same radial increase – from low at the center to higher outward – highlighting the effectiveness of our multi-component fitting approach.

Appendix C: Fits of individual Voronoi bins

In this appendix, we illustrate the kinematic fits to all the individual Voronoi-binned spectra in Figs. C.1 and C.2.

Fig. C.1. Individual kinematic fits for the first 18 out of the 32 Voronoi-binned spectra. The best-fit gas lines are subtracted for better visualization of the Ca II triplet features (marked with orange lines). The blue-shaded region marks where we boosted the noise levels around the Hα and N II lines. The gray bars represent the observed data, with the width marking the bin width and the height marking the original 1σ noise level. The vertical gray lines represent the boosted 1σ noise level. These 1σ noise levels for each bin are estimated from random sampling given the covariance matrix. The measured velocity dispersion σlos is also annotated in each panel.

Fig. C.2. Individual kinematic fits for the last 14 out of the 32 Voronoi-binned spectra. The best-fit gas lines are subtracted for better visualization of the Ca II triplet features (marked with orange lines). The blue-shaded region marks where we boosted the noise levels around the Hα and N II lines. The gray bars represent the observed data, with the width marking the bin width and the height marking the original 1σ noise level. The vertical gray lines represent the boosted 1σ noise level. These 1σ noise levels for each bin are estimated from random sampling given the covariance matrix. The measured velocity dispersion σlos is also annotated in each panel.

All Figures

Fig. 1. Imaging of the system RXJ1131−1231. Left panel: HST-ACS image in the F814W band. The four quasar images are labeled 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. The North and East directions, along with a 1″ scale, are also illustrated. The $3\stackrel{″}{.}2\times 3\stackrel{″}{.}1$ yellow squares represent the dithered 2 × 2 mosaic pattern of the NIRSpec IFU, which covers a $6\stackrel{″}{.}1\times 5\stackrel{″}{.}55$ field of view over the system. Right panel: “White-light” image of the system from the JWST-NIRSpec datacube, summed within 8255−8890 Å in the lens galaxy’s rest frame. The white bar represents 1″, and the red and yellow arrows point to the east and north, respectively. The blue contours show the region within which the spaxels are summed to obtain a high-S/N spectrum of the host galaxy. The black regions around the quasar images were used to extract summed spectra for modeling the quasar’s emission line components. HST and JWST logo credits: NASA.

In the text

Fig. 2. Lens modeling with the image obtained by summing the datacube within 8700−8800 Å in the lens galaxy’s rest frame. This wavelength range covers the Ca triplets of the lens galaxy. Top row, from left to right: Illustration of the data, the optimized-model-based reconstruction of the data, and the normalized residuals. Bottom row, from left to right: Illustration of the initial PSF model from STPSF, the iteratively reconstructed PSF, and the fraction of the lens light compared to the total light in the image. The reconstructed PSF will be necessary for dynamical modeling, and we use the lens light fraction to obtain the S/N of the lens galaxy for Voronoi binning.

In the text

	Fig. 3. Voronoi binning map. Left panel: Voronoi bins outlined over the white-light image of the system. Middle panel: Bin numbering. Right panel: Total sS/N for each Voronoi bin. The horizontal dashed line marks 90 Å−1/2, which we set as the minimum target for each bin. Only the last Voronoi bin, which contains the satellite galaxy, is slightly below this sS/N target.
In the text

Fig. 4. Fitting of the Ca II region of the quasar host galaxy. The gray bars represent the data points, with the width of each bar corresponding to the pixel size and the height representing the original ±1σ noise level. The gray lines attached to the bars represent the boosted uncertainty levels to achieve χred2 = 1. The red line illustrates the best-fit model. The Ca II triplet features are marked with blue lines.

In the text

Fig. 5. Fitting of the quasar lines at images A, B, C, and D, respectively, from top to bottom. The observed spectra are shown with gray bars, with the width corresponding to the bin width and the height representing the ±1σ noise range. The red lines show the best fit. The individual narrow and emission lines from the best fit are shown in purple, and the lens galaxy’s stellar continuum is shown in emerald. We adopted the best-fit gas lines with fixed relative amplitudes for all three images as free linear components in our kinematic fitting of the spectra across the entire field of view.

In the text

Fig. 6. Best fit (red line) of the Voronoi-binned spectra in bin 22 as an example. The gray bars represent the observed data, with the width indicating the bin width and the height representing the original 1σ noise level. The vertical gray lines attached to the bars represent the boosted 1σ noise level. These 1σ noise levels for each pixel are estimated from random sampling given the covariance matrix. The Ca II triplet wavelengths from the lens galaxy are marked with orange lines, and the wavelengths of the host galaxy’s emission lines are marked with blue lines. The blue-shaded region marks where we boosted the noise levels around the Hα and [N II] lines. The emerald line shows the contribution from the lens galaxy’s stars. The sets of narrow emission lines from the quasar host galaxy (including spurious spikes identified close to 8300, 8700, and 8900 Å in this case) are shown in purple.

In the text

Fig. 7. Velocity dispersion (top row) and mean velocity (bottom row) measurements in 32 Voronoi bins. The first column shows the measurement for each case, and the second column shows the measurement uncertainties. The bin containing the satellite galaxy (i.e., bin 32) is outlined in white. The bins closer to the center have smaller uncertainties, thanks to less contamination from the quasar and its host galaxy. The systemic velocity of 213 km s−1 was subtracted from the mean velocity map. The mean velocity map indicates that the lens galaxy is a slow rotator, with small mean velocities (≲50 km s−1) and a lack of a bipolar velocity profile.

In the text

	Fig. 8. “Fractional” covariance matrix. In this matrix, the element with indices i, j is [cov(σlos, i, σlos, j)/σlos, iσlos, j]1/2. The diagonal elements represent the combined statistical and systematic uncertainties. The mean for the absolute values of the off-diagonal elements is 1.14 ± 0.04%.
In the text

Fig. 9. One-dimensional profile of the rms velocity vrms ≡ (σlos2 + vlos2)1/2, where vlos is the line-of-sight mean velocity after subtracting the mean systemic velocity of 213 km s−1. Left panel: One-dimensional vrms profile for JWST-NIRSpec measurements from this paper. Middle panel: One-dimensional vrms profile for the Keck-KCWI measurements from Shajib et al. (2023). The vertical error bar represents the total uncertainties, including both statistical and systematic contributions, and the horizontal error bars denote the annulus widths. The dashed lines represent the best-fit model predictions from JWST-only data (red), KCWI-only data (blue), and from a joint fit (black), and the associated shaded regions illustrate the 1σ extent of the predictions from the model posterior. The horizontal arrows represent the FWHMs of the corresponding PSF in each panel. Right panel: One-dimensional light profile from double Sérsic fits performed on the NIRSpec white-light image (red line) as part of the lens modeling in this paper, and from that, on the HST imaging in the F814W filter (blue line) from Shajib et al. (2023). The vertical dashed line marks the $1\stackrel{″}{.}3$ transition radius, where we stitch the light profiles with IR in the inner region and the visible one in the outer region that is used in the dynamical model of the NIRSpec 1D profile. From a likelihood-ratio test, we find that the JWST and KCWI measurements are consistent at ∼1.2σ (i.e., the likelihood-ratio test statistic’s p-value is ∼0.22).

In the text

Fig. 10. Comparison of the parameter posteriors of the dynamical modeling performed with JWST-NIRSpec-only data (red), Keck-KCWI-only data (blue), and both combined (black). The four free parameters in our dynamical model are the ratio of the tangential and radial velocity dispersion σt/σr ≡ (1 − βani)1/2 with βani being the amplitude of a spatially constant anisotropy profile, the MST parameter λmst, the main deflector’s Einstein radius θE, the power-law mass profile’s logarithmic slope γ, and the super-massive black hole’s Einstein radius θE, bh. We have imposed a prior on the total Einstein radius (θE2 + θE, bh2)1/2 based on the lens model posterior of Suyu et al. (2013), which dominates the θE posterior illustrated here. All of them are consistent at ≲1σ confidence levels. The JWST NIRSpec’s superior spatial resolution also allows for the exclusion of a large fraction of the prior volume in the anisotropy parameter (σt/σr) that is largely unconstrained by the Keck-KCWI measurements.

In the text

Fig. A.1. Fit wiggle model used to correct for the resampling noise. The data points show the signal of wiggle d(λ)/T(λ) for the central spaxel on the lens galaxy, where T(λ) is obtained from the best-fit model m(λ). The orange line shows the W(λ) function in the best-fit m(λ), with the orange shaded region illustrating the 1σ confidence region. The vertical gray dashed lines show the wavelength range containing the Ca II triplets of the lens galaxy that we fit to measure the lens galaxy’s kinematics. The gray data points are rejected outliers using the false discovery rate method.

In the text

Fig. A.2. Spaxels cleaned for wiggles (within the regions enclosed with white borders). Only the spaxels within a circle of $0\stackrel{″}{.}7$ radius from the lens galaxy center were considered to be cleaned. Then, we only cleaned the spaxels for which the wiggle was detected at a confidence level of more than 3σ and the best-fit wiggle model explains more than half of the variance or scatter observed in the wiggle signal. The wiggles are more prominent in spaxels with more flux. As a result, the white contours above primarily encompass the central region of the lens galaxy and the satellite galaxy.

In the text

Fig. A.3. Impact of wiggle cleaning on the measured kinematics. The illustrated test is for a single model setup out of the 48 combinations described in Sect. 4.5. The first column displays the fractional differences in the measured dispersions (top row) and mean velocities (bottom row) within the Voronoi bins, while the second column shows these differences normalized by the statistical uncertainty.

In the text

	Fig. B.1. Expected quasar contamination in each Voronoi-binned spectrum from the best-fit lens model (left) and the effective levels from our PPXF fits (right). The effective contamination is consistent across bins and follows the same radial increase – from low at the center to higher outward – highlighting the effectiveness of our multi-component fitting approach.
In the text

Fig. C.1. Individual kinematic fits for the first 18 out of the 32 Voronoi-binned spectra. The best-fit gas lines are subtracted for better visualization of the Ca II triplet features (marked with orange lines). The blue-shaded region marks where we boosted the noise levels around the Hα and N II lines. The gray bars represent the observed data, with the width marking the bin width and the height marking the original 1σ noise level. The vertical gray lines represent the boosted 1σ noise level. These 1σ noise levels for each bin are estimated from random sampling given the covariance matrix. The measured velocity dispersion σlos is also annotated in each panel.

In the text

Fig. C.2. Individual kinematic fits for the last 14 out of the 32 Voronoi-binned spectra. The best-fit gas lines are subtracted for better visualization of the Ca II triplet features (marked with orange lines). The blue-shaded region marks where we boosted the noise levels around the Hα and N II lines. The gray bars represent the observed data, with the width marking the bin width and the height marking the original 1σ noise level. The vertical gray lines represent the boosted 1σ noise level. These 1σ noise levels for each bin are estimated from random sampling given the covariance matrix. The measured velocity dispersion σlos is also annotated in each panel.

In the text