Astronomy & Astrophysics manuscript no. main March 4, 2025
©ESO 2025
TDCOSMO: XX. WFI2033–4723, the First Quadruply-Imaged Quasar Modeled with JWST Imaging
D. M. Williams⋆,1, T. Treu1, S. Birrer2, A. J. Shajib,⋆⋆,3, K. C. Wong4, T. Morishita5, T. Schmidt1, and M. Stiavelli6
1 Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095, USA
2 Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA
3 Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA
4 Research Center for the Early Universe, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan
5 IPAC, California Institute of Technology, Pasadena, CA 91125, USA
arXiv:2503.00099v1 [astro-ph.CO] 28 Feb 2025
6 Space Telescope Science Institute, Baltimore, MD 21218, USA
Received 00 Month 0000; accepted 00 Month 0000
ABSTRACT
Gravitational time delays offer unique, independent measurements of the Hubble constant, H0. Precise measurements of H0 stand as one of the most pressing challenges in modern cosmology, and to do so with time delays requires precise lens models. While much work has focused on streamlining the modeling process to keep pace with the erumpent discovery of strongly-lensed systems, a critical step toward reducing uncertainty in H0 comes from increasing the precision of individual lens models themselves. In this work, we demonstrate that the unprecedented imaging capabilities of JWST make this goal attainable. We present the first lens model for time-delay cosmography derived from JWST data, applied to the quadruply-imaged quasar WFI2033–4723. While the primary source of systematic uncertainty in time-delay cosmography overall is currently the mass-sheet degeneracy (MSD), the sensitivity of models to this MSD varies depending how the point spread function (PSF) errors are mitigated. As the PSF is also the primary source of uncertainty in lens modeling, we focus on a comparison of different PSF modeling methods. Within the context of power-law models, we recover results in agreement with previous Hubble Space Telescope (HST)-based models, but with better precision of key lensing parameters through implementation of new PSF modeling techniques. Despite the record-holding precision of this system’s HST modeling, we achieve an additional 22% increase in precision of the Fermat potential difference, directly reducing uncertainties
of cosmological inference. These results would produce a 3% (1σ of the lens modeling error) shift of H0 towards a higher value for this lens, if one were to keep all else constant. This work substantiates the groundbreaking potential of JWST for time-delay cosmography and lays the groundwork for modeling systems previously too faint to provide meaningful constraints on H0.
Key words. Gravitational lensing: strong – Methods: data analysis, statistical – Galaxies: active, distances and redshifts – cosmological parameters – distance scale
-
Introduction
The expansion rate of the universe remains one of the most hotly debated topics in cosmology after two decades of measurements with improving statistical precision. Measurements of the early universe, extrapolating to what the current value of the expan- sion rate should be given today’s standard cosmological model ΛCDM, give values in tension with direct measurements of the late universe. Hope that the tension would dissolve in time is proving to be misplaced—as the precision of both the early- and late-universe measurements have increased, so too has their ten- sion. For a more comprehensive review, see, e.g. Verde et al. (2019); Shah et al. (2021); Abdalla et al. (2022).
The confirmation of the "Hubble tension" would have pro- found implications for cosmology, requiring adjustments to or even a complete revision of the standard model (Knox & Millea 2020; Efstathiou 2021, e.g.). Considerable effort has been spent evaluating sources of systematic uncertainties in both methods, none proving successful (Riess et al. 2019, 2021, 2022; Freed- man et al. 2019, 2020; Dainotti et al. 2021; Mortsell et al. 2022). It is clear that multiple independent methods with comparable
⋆ e-mail: devon@astro.ucla.edu
⋆⋆ NHFP Einstein Fellow
precision will be needed to achieve a final verdict on residual systematics error and settle the Hubble tension once and for all.
One such method is time-delay cosmography, first discussed in Refsdal (1964) (for recent reviews, see, e.g. Treu & Mar- shall 2016; Treu et al. 2022). It provides a one-step measure- ment of H0, independent of methods that rely on the local dis- tance ladder for late-universe measurements or on sound horizon physics, such as CMB measurements, for early-universe mea- surements (e.g., Vanderriest et al. 1989; Keeton & Kochanek 1997; Oguri 2007; Suyu & Halkola 2010). In this work, we focus on time-delay cosmography applied to lensed quasars. Quadruply-imaged quasars (hereafter quads), where light from a quasar is bent around a lensing galaxy to produce four dis- tinct images, provide three independent measurements of the ex- pansion rate along with tight constraints on the lensing galaxy’s mass distribution.
Until the last decade, the precision of time-delay cosmog- raphy was limited by the small number of known systems, and further by the number of those systems with reliable time de- lays measured. Fortunately, there has been an explosion of newly discovered strong lenses, resulting in a commensurate effort to gather the necessary ancillary data and to quicken the lens mod- eling process. Key advances include: systematic time-delay mea-
surements (spearheaded by the COSMOGRAIL project: e.g., Millon et al. 2020; Bonvin et al. 2019), automated modeling techniques (e.g., Shajib et al. 2019; Schmidt et al. 2022), con- sistent analyses of lens environments (e.g., Wells et al. 2024b), and improved kinematic constraints (e.g., Tan et al. 2024). Un- fortunately, ground-based seeing-limited images cannot be mod- eled to the precision required for cosmography. The Hubble Space Telescope (HST), with its stable PSF, has been the main workhorse for these kinds of models (e.g., Suyu et al. 2010a; Birrer et al. 2019), even though some breakthroughs have been achieved in the use of adaptive optics assisted data (e.g., Chen et al. 2019).
In this context, JWST plays multiple important roles: first, it holds the potential to further improve the precision of lens mod- els; second, it provides an entirely new means of verifying sys- tematics in previous models that are based on other telescope data (namely HST); third, with its superior resolution and sen- sitivity it will be crucial to exploit the large number of systems about to be discovered, which will likely be fainter than the ones analyzed so far (Oguri & Marshall 2010); and fourth, the ex- tended wavelength range will enhance the detection of the ex- tended emission (arcs) from the quasar host galaxy, providing tighter constraints on the deflector’s mass distribution.
E
N
1"
Fig. 1. JWST NIRCam color image of WFI2033-4723 (GTO Program #1198; PI: Stiavelli). The F115W band is mapped to the luminosity and blue, F150W to cyan, F277W to yellow, and the F356W band to red. Near the bottom image, there is a bright spot on the lower-right side that maps back to other bluer regions of the ring, indicating it is a multiply-
∼ ∼
~
imaged star-forming region (SFR) contained in the z 1.6 host. Prelim- inary models place the physical size at 150pc located 2kpc from the center of the galaxy (assuming flat CDM cosmology with h = 0.7 and m = 0.3).
In this paper, we present the first JWST-based quadruply- imaged quasar lens model for cosmology. We focus on the lens system WFI2033–4723, discovered by Morgan et al. (2004) and previously modeled based on HST data (Rusu et al. 2020). This system is particularly interesting as the lensed host galaxy is detected by both HST and JWST, but only clearly resolved by
JWST, as shown in Fig. 2. The figure also illustrates the com- plexity of JWST’s PSF, which represents a major challenge of this work. Further, the system is fairly complex, having a satel- lite and multiple nearby galaxy perturbers with non-negligible impact on the lens models (and so require explicit modeling) (Sluse et al. 2019; Rusu et al. 2020). The goals of this paper are thus:
-
To verify whether the HST and JWST models agree for this system, justifying our handling of systematic effects in the HST models
-
To quantify the improvements in precision that JWST pro- vides for lens modeling
-
To demonstrate sufficient handling of the PSF to facilitate
these improvements
-
To determine if the superior angular resolution and im- proved extended-emission/arc detection of JWST yields im- proved constraints on nearby perturbers
As the goal of this paper is to make a direct comparison of strong lens modeling based on HST and JWST, we limit our- selves to the core analysis without performing a full determina- tion of the Hubble constant from this study. While this work can be used for such a study, the full analysis will be performed in a proceeding paper. Among the simplifications, we work within the context of a power-law model and do not consider the mass sheet degeneracy. We also do not consider updates to the external convergence estimates that have been implemented by our col- laboration since the Rusu et al. (2020) analysis. These choices make it more straightforward to perform a direct comparison with Rusu et al. (2020). For these reasons, we compare our re- sults in terms of Fermat potential differences between the images
with the longest expected time delay ∆τBC (for brevity, we will refer to this quantity as "the Fermat potential" in the text). This
is the quantity more directly connected to the lens model and yet close enough to H0 to permit an intuitive understanding of our results: all other things being equal, this quantity is directly pro- portional to H0. However, we stress that all other things are not equal, since other factors contribute to H0, for example the exter- nal convergence, determination of the MSD parameter from the velocity dispersion, and the other time delays of the quad. There- fore, our result should not be taken out of context and plugged into the analysis of Rusu et al. (2020) directly to update their measurement of H0.
An important difference with respect to the previous analysis is that we use a different strong lensing code: Rusu et al. (2020) used glee (Suyu & Halkola 2010), while we use lenstronomy (Birrer & Amara 2018a; Birrer et al. 2021). The two codes have been tested extensively and compared directly (e.g., Shajib et al. 2022). The comparisons show that they give consistent results even though they differ in the way they parametrize important components, such as the host galaxy surface brightness distribu- tion. When possible, we use the same parametrization as Rusu et al. (2020) to facilitate the comparison.
The main result of this paper is that JWST’s PSF can be mod- eled with sufficient accuracy to enable the development of the next generation of lens models for time-delay cosmography. We show that JWST yields better precision than HST, and that the results agree within the statistical errors. We also show that with JWST’s superior resolution and sensitivity, the mass distribution of the satellite galaxy is better constrained than with HST.
The paper is organized as follows. First, in Sect. 2 we outline the analysis steps. Then in Sect. 3 we present the data. In Sect. 4 we describe our PSF reconstruction methodology. In Sect. 5 we describe our modeling choices. In Sect. 6 we discuss how our
models are combined in the final inference. In Sect. 7 we present our results and compare them with the previous HST model. In Sect. 8 we provide a brief summary.
When necessary we adopt a standard flat ΛCDM cosmology with H0 = 70 km s−1, Ω,m=0.3, Ωb = 0.05, although we note that this choice is irrelevant from the point of view of the com- parison presented here, since lens modeling is done in angular coordinates.
-
-
The big-picture goal of this paper is to assess how advancements in space-based imaging quality impact the precision and accu- racy of model-derived parameters. Specifically, we investigate how the transition from HST to JWST imaging influences the Fermat potential differences, ∆τ, which are crucial outputs for
cosmographical inference of H0. To this end, we first review the
density of the universe, such that a pointing with 5% more than the cosmological average yields κ = 1.05. In other words,
κ ➟ Σ , (5)
Σcrit
with the critical density given by
c2Ds
Σcrit = 4πGDds Dd . (6)
∇
It is also related to the deflection potential, as 2ψ = 2κ.
The difference in travel time for this lensed ray and an un-
lensed ray is given by the excess time-delay,
t(θ, β) = D∆t τ(θ, β), (7)
c
where c is the speed of light and D∆t is the time-delay distance
(Refsdal 1964, Schneider et al. 1992, Suyu et al. 2010b).
The time-delay distance is defined as
lensing theory and discuss the observables in Sect. 2.1, review
DdDs
D ≡ (1 + z ) , (8)
inherent degeneracies and systematic concerns of the analysis in ∆t
d Dds
Sect. 2.2 and how we address them in Sect. 2.3, and finish with the outline of our Bayesian analysis in Sect. 2.4.
-
This section will sprint through the theoretical groundwork for time-delay cosmography. For more of a stroll, we refer to recent reviews, e.g., Treu & Marshall (2016) or Suyu et al. (2018).
As light is deflected by gravitational fields, we can occasion- ally observe multiple rays of light traversing different paths yet originating from the same object. This special class of objects are dubbed “multiply-imaged". The combination of these two factors allows one to calculate distances at cosmological scales. How our image is distorted is specified by the lens equation,
β = θ − α(θ), (1)
where zd is the redshift of the deflector (lens) and Dd, Ds, and Dds are the angular diameter distances to the deflector, source, and between the lens and source, respectively. The significance of the time-delay distance is that within each distance term lies
an intrinsic H0−1. This leaves us
D∆t ∝ H0−1. (9)
D∆t has units of distance, and it has weak dependence on other cosmological parameters.
Now, we assume the background source and foreground de- flector galaxy are sufficiently aligned such that multiple images are observed. As different images require light to travel different paths, their excess time-delays will differ. The time-delay be- tween two images, say θi and θ j, is thus the difference of their excess time-delays,
c
c
∆ti j = D∆t .τ(θi, β) − τ(θj, β). = D∆t ∆τi j. (10)
in the source plane (β = (β , β
1 2
with the lensed position in the lens plane (θ = (θ , θ )), position
α).
Therefore, constraints on the observed time-delay (∆ti j) and
1 2)), and the deflection angle (
the predicted Fermat potential difference (∆τi j) can be converted
If we assume the deflection occurs in a single lens plane, we can quantify the observed delay between images with the Fermat potential,
into constraints on the time-delay distance (D∆t), which is an inference on H0 through Eq. (9).
If a source is variable on sufficiently short timescales, it is
τ(θ, β) ≡
(θ β)2
" −
#
2 − ψ(θ)
, (2)
possible to measure a time-delay, ∆ti j, by extended monitoring of the image fluxes (e.g., Vanderriest et al. 1989; Schechter et al. 1997; Fassnacht et al. 1999, 2002; Kochanek et al. 2006; Courbin
with the lens potential (ψ(θ)). The Fermat potential is, up to an affine transformation, the light travel time along a ray starting at
β, passing through θ in the lens plane, and arriving at the ob- server. The time required for light to reach the observer from one of these images is influenced by two factors. First is the spa- tial path length the light must traverse, and second is the grav- itational potential the light experiences along that path (via the ‘Shapiro delay’ Shapiro 1964).
The lens potential (ψ(θ)) is defined such that it solves two conditions. First, for the deflection angle
α(θ) = ∇ψ(θ), (3)
and second, for the convergence (κ)
et al. 2011; Bonvin et al. 2019).
The Fermat potential difference, ∆τi j, requires three addi- tional parameters for Eq. (2) (as θi and θ j are simply the positions observed on the sky). Thankfully, neglecting MSD, the source position, β, and the lens potentials at the image positions, ψ(θi)
and ψ(θj), can be determined by accurately modeling the mass of
the system in imaging data; that data we denote by dJWST. Then, with both ∆ti j and ∆τi j, Eq. (10) provides a measurement of D∆t. With an assumption of a cosmological model, D∆t can be turned into an inference on H0.
Thus, the necessary ingredients for cosmographic inference of H0 from strong lensing can be broken into four components. First is the modeling of imaging data, dJWST, which is the focus of this paper. Second is the measured time-delays, ∆ti j, observed
1 2 from combinations quasar’s images (i and j). Third is the stellar
κ(θ) = 2 ∇ ψ(θ). (4)
The convergence is just the surface mass density in the im- age/lens plane at some pointing (Σ) scaled by the critical surface
kinematic information, σP, from spectroscopic data of our lens. The fourth, and final ingredient, is wide-field spectroscopy and imaging, denv, providing environmental information. These are
-
Additional sources of error and degeneracies
-
Internal mass sheet degeneracy
Currently, the primary source of uncertainty in the Fermat potential—and time-delay cosmography at large—is the mass- sheet degeneracy (MSD) discussed in Falco et al. (1985) (more recently Wong et al. 2020, or Gilman et al. 2020 for substructure discussion, or Birrer et al. 2020 for environment discussion). In summary, by adding an infinite “sheet of mass” with constant uniform surface density, one can find the exact same image po- sitions given a different source position, size, and luminosity, completely changing the physical geometry of the system and its associated Fermat potential.
While the time-delay (∆ti j) is an observable, the relative Fer- mat potential (∆τi j) is not. Instead, the Fermat potential is in- ferred from constraints on two particular features in the imaging data: the quasar image positions and the extended distortions in the lensed arcs. Consequently, the MSD allows degeneracy in
our observed image positions to seep into our Fermat potential, thereby imposing limits on its precision and, via Eq. (9) , on H0 (e.g., Falco et al. 1985; Kochanek 2002; Saha & Williams 2006; Schneider & Sluse 2013, 2014; Birrer et al. 2016; Unruh et al. 2017; Birrer 2021).
As the MSD is inherent in the lens model (Saha 2000; Saha & Williams 2006), we first note that the MSD is a mathematical degeneracy as opposed to a physical degeneracy. More formally, the MSD is a multiplicative transform of the lens equation that yields the same image positions despite a linear source displace-
ment and transformation of the convergence field, κ. In fewer words, but more Greek, we first displace the source by some scalar λ:
β → λβ, (11)
and transform the convergence field through:
κ(θ) → λκ(θ) + (1 − λ). (12)
Then, solving the lens equation with these quantities, we will find the image positions at:
θ → θ. (13)
In other words, our direct observable in the imaging data—the image positions—would look the same despite moving the source and changing the convergence field. This is problematic, as the linear source displacement in Eq. (11) would yield a mea- surement of:
H0 → λH0. (14)
This is currently the leading source of uncertainty in time-delay measurements of H0 today, and it is broken with the non-lensing constraints on the mass distribution, namely stellar velocity dis- persion Birrer et al. (2020); Shajib et al. (2023), or by constraints on the absolute luminosity of the source when it is known. As discussed in the introduction, this paper is concerned with in- formation coming from the lens model itself, and therefore we can neglect the mass sheet transformation, keeping in mind that it needs to be considered when converting the inference from Fermat potential to H0.
-
Multiplane lensing
In addition, the observed time-delays will depend on any galaxy (or massive object) having a significant impact on the deflec- tion of our source’s light, and these need not be at the same red- shift as our lens. In this generalized case, we must combine the angular diameter distances of the different lens planes, which is done with the multi-plane lens equation (e.g., Blandford & Narayan 1986; Kovner 1987; Schneider et al. 1992; Collett & Auger 2014; McCully et al. 2014). In our previous assumption, our time-delays were proportional to a single, unique time-delay distance. Under the multi-plane lens equation, this is no longer true. Thankfully, for quadruply-imaged quasars, it very often holds that the mass in a single plane dominates the lensing effect Suyu et al. (2020). This validates the single-lens approximation, as the observed time-delays are dominated by the time-delay dis- tance (Eq. (8)) at the redshift of the lensing galaxy. This has also been found to be the case for this lens, WFI2033-4723 (Rusu et al. 2020). This single time-delay distance is referred to as the effective time-delay distance.
-
External convergence
∆t
An additional correction must be made to validate the single- plane approximation, addressing the lens potential that the light rays experience. We must take into account the additional focus- ing and defocusing that external mass causes, which modifies the measured time-delays (e.g., Seljak 1994). To avoid a biased in- ference of the time-delay distance, we apply a correction to our modeled value, Dmodel, to get the true D∆t with:
Dmodel
D∆t = ∆t . (15)
1 − κext
–
This is equivalent to adding (or removing) an external sheet of mass in the lens plane to account for the focusing (or defocusing) from the mass unaccounted for in the lens plane. This correction is valid for cases where the LOS (line-of-sight) perturbers are small enough to ignore the effect of higher-order terms (Keeton 2003; McCully et al. 2014). With Eq. (9), we note this correction yields a multiplicative factor of (1 κext) for inferences on H0.
Eq. (12), however, tells us that determining the value of κext is generally impossible from the lens model alone due to the MSD.
The approach typically taken includes estimating the mass along the LOS combined with the assumption that the deflector’s mass profile drops to zero at large radii (e.g., Fassnacht et al. 2006; Momcheva et al. 2006, 2015; Williams et al. 2006; Wong et al. 2010; Wells et al. 2024b). Higher-order effects must be explicitly taken into account for perturbers that are extremely massive or particularly close to the lensing galaxy (McCully et al. 2017).
-
Mass density slope/time-delay degeneracy
An additional degeneracy exists between the lensing galaxy’s radial mass slope and the time-delay distance (e.g., Kochanek 2002; Treu & Koopmans 2002; Wucknitz 2002), even within the context of power-law models. This degeneracy is broken in two ways. First, thanks to JWST, we have imaging of the extended arcs that are barely noticeable in HST data (see Fig. 2), whose modeling directly constrains the radial mass slope. Second, we also combine the lensing data with stellar kinematic data (e.g., Treu & Koopmans 2002; Koopmans et al. 2003; Auger et al. 2010; Suyu et al. 2014; Yıldırım et al. 2020).
-
-
Stellar kinematics plays a crucial role in mitigating the effects of degeneracies (Treu & Koopmans 2002; Shajib et al. 2018, 2023; Knabel et al. 2024). With measurements of the stellar velocity dispersion, we can independently probe the 3D mass distribution of our deflector galaxy (after de-projecting from the observed convergence near the center of our lens, κobs, cen). I.e.,
obs, cen
κ
σ −−−−→ Φ(r) → M, (16)
where the observed convergence near the center of the deflector galaxy, κobs, cen, allows us to de-project the velocity dispersion, σ, to the gravitational potential, Φ(r).
To this end, we employ the spherical Jeans equation,
To correct for this effect, we sample κext and weight our mod- els by their ability to recover the observed velocity dispersion of the main deflector, σobs, D (see Sect. 6.2).
We note full cosmographic analyses will also feature an ad- ditional degree of freedom in models with a λ term in this equa- tion. As the purpose of this paper is to evaluate the improve- ments offered exclusively by imaging data, we focus our efforts on a direct power-law comparison, as imaging alone is unable to directly constrain the MSD.
-
To predict a time-delay under a given cosmology, the lens model infers the Fermat potential difference between two image pairs,
2 which we define as ∆τ. We aim to describe Pr(∆τ | O), the con-
d l(r)σr(r)
{ }
dr +
2βani(r)l(r)σr(r)2
r = −l(r)
dΦ(r)
dr , (17)
ditional probability of the Fermat potential difference given our imaging and kinematic observables O ➟ Oimg, Okin . We do not include the spectroscopic and photometric data used for the lens
with the 3D luminosity density, l(r), the radial velocity disper- sion, σr(r), and the anisotropy parameter that relates σr to the tangential velocity dispersion σt. We define it as
σ (r)
2
βani(r) ➟ 1 − t . (18)
environment, as the external convergence inference can be ac- counted for with a separate, explicit prior, Pr(κext).
|
For a full description of Pr(∆τ O), we must first address the functional dependencies of the Fermat potential difference
∆τ(ξ, κext). First is the external convergence κext and second is
r
σ2(r)
our set of model parameters ξ ➟ {ξlens
, ξlight
, rani
}. Given our ob-
The quantity we actually observe is the luminosity-weighted LOS velocity dispersion. We get this from the Jeans equation, as
servables, we can apply Bayes’ theorem to rewrite
,2
Pr(ξ, κext | O)
σ
2
los
(R) = I(R)
Kβ
2G , ∞ r l(r)M(r)
R
R
∝ Pr(O | ξ, κext) × Pr(ξ, κext)
=
K
r dr, (19)
with the gravitational constant, G, the surface brightness, I(R), and the enclosed 3D mass within a radius r, M(r) (Mamon & Lokas 2005). Our function, β, depends on the parameterization
of βani(r). We adopt the Osipkov-Merritt parameterization, given
by
r2
× Pr(ξ, κext | S ) dS dDs/ds . (24)
Pr.O | ξ, κext, S , Ds/ds
We have introduced S , the set of hyper-parameters used by the lens model to describe Oimg (e.g., the number of parameters de- scribing the source, the set of pixels upon which the image likeli- hood is evaluated), and Ds/ds, which represents the distance ratio
βani(r) = r2
+ r
2
ani
, (20)
Ds/ds ➟ Ds/Dds. We then marginalize on these two parameters. We note a relatively recent change to this method of Bayesian
dent data components, rewriting Eq. (24) as
with the scaling radius, rani (Osipkov 1979; Merritt 1985a,b). This sets our function to
analysis compared to similar works, where the external conver- gence is now (once again) independent of the model-predicted external shear (see Sect. 3.4). This splits the joint likelihood of
Kβ
u ➟ R
=
r
κext into its own prior. We also separate O into its two indepen-
ani
u2 + 1/2
u2 + u2 , u2 − 1
, Pr.O | ξ, κext, S , D
s/ds
u
,
ani tan−1
(uani + 1)3/2 u
2
ani
+ 1
– u
1/2 1
2
× Pr(ξ | S ) × Pr(κext) dS dDs/ds
2
ani + 1
1 − u2 , (21)
= ,2
Pr Oimg | ξ, S
with uani ➟ rani/R (Mamon & Lokas 2005). We can then find the observed aperture-averaged velocity dispersion via
× Pr(Okin | ξ, κext, Ds/ds)
× Pr(ξ | S ) × Pr(κext) dS dDs/ds . (25)
2 ,ap
.I(R)σ2
(R). ∗ S dxdy
(22)
Next, we convert the following sub-integral
σap =
,ap
,
los
I(R) ∗ S dxdy
, Pr O
Pr ξ | Oimg, S
ξ, S Pr(ξ
S ) Pr(S ) dS
,
where we integrate over the aperature, , , and convolve with the
img | |
ap
We ∗nSote the external convergence will modify our model-
Pr Oimg | S
Pr(S ) dS . (26)
seeing, .
ap, model
=
predicted velocity dispersion, σap, model through
This form is much more convenient, as sampling it gives us the
σ
2
ap, true
= (1 − κext)σ2
. (23)
model evidence, Pr Oimg | S
JWST NIRCam/F115W
[1804s @ 1.15 m]
E
N
1"
HST ACS/F814W
[2085s @ 0.814 m]
E
N
1"
G3
(z = 0.6575)
B
D
(z = 0.6575)
X
(z = 0.6575?)
A2
Fig. 2. Comparison of the JWST data used in this analysis, NIRCam F115W, to the HST imaging, ACS/F814W. The ring structure is maxi- mized in each visualization, as it provides the tightest constraints on the deflector’s radial mass profile.
A1
C
G2
(z = 0.7450)
-
-
As the focus of this paper is to explore the improvements that JWST imaging has to offer for cosmography-grade lens mod- eling, we use the same data products as the first cosmography model of this system (Rusu et al. 2020, hereafter H0LiCOW XII), replacing the HST imaging with that of JWST. First, we
G7
(z = 0.6575)
E
N
5"
–
review the lens system WFI 2033–4723 (J2000: 20h33m41.9s, 47◦23′43′.′4; hereafter WFI2033) in Sect. 3.1, followed by a discussion of the JWST NIRCam data used in Sect. 3.2. Next,
we summarize the ancillary data, starting with the spectroscopic data used for measurements of velocity dispersions in Sect. 3.3. Finally, we go over the wide-field imaging data and how it is used to constrain the impact of the environment through κext in Sect. 3.4.
-
±
The quadruply-imaged quasar WFI2033 was first discovered by Morgan et al. (2004), wherein they applied color cuts to the Wide-Field Imager (WFI) on the MPG-ESO telescope. The red- shift of the source was measured to be zs = 1.662 by Sluse et al.
(2012), with the deflector at zd = 0.6575 0.0002 (Sluse et al.
2019) (hereafter H0LiCOW X). Its environment is relatively
±
complicated, with at least six galaxies within 20 arcsec, three of which needing to be explicitly included in the lens model1. The deflector is also in a group at z = 0.6588, with a scale ra-
dius predicted to be rs,g = 32.0 8.0 arcsec, calculated by as-
suming its virial mass and radius (Rusu et al. 2020) (hereafter
H0LiCOW XII). In addition, there exists a second group in its
foreground (z = 0.4956) with its scale radius similarly measured
Fig. 3. Environment of WFI2033 with black labels denoting the main deflector (D), its nearby satellite (X) and the three galaxy perturbers, (G2, G3, and G7) along with their redshifts. The quasar images are la- beled in blue.
resolution than otherwise possible. Two other quadruply-imaged systems were also observed, whose NIRCam models are planned in following publications.
~
The NIRCam observations were done with 4 primary INTRAMODULEX dithers with the SHALLOW2 readout pattern, pro- ducing a slightly higher angular resolution ( 0′.′0307 in the F115W band) than the default pixel size. The bright quasar im- ages were not saturated and as such did not require masking.
~
This is extremely beneficial for our point source modeling, later discussed in Sect. 4. There were 9 groups per iteration, with a total exposure time of 1800 sec. The uncertainty on the flux in each pixel was estimated from the science image, and is a combination of the resampled Poisson, read noise, and flat-field variances added in quadrature. While this data also includes the F150W, F277W, and F356W bands, one of the largest bottle- necks of modeling this high resolution data is computational time. For this reason, we choose to exclusively model the F115W band, as it offers the highest spatial resolution which surpasses
at r
s,g
= 34.8 ± 9.3 arcsec.
HST at equivalent wavelength. Our lens model is based on a 4.6′′ (150 pixel) square cutout of the reduced data, with a 0′.′15 mask applied to the very center of the deflector.
-
WFI2033 was observed on JWST with both NIRCam imaging and NIRSpec IFU spectroscopy on September 19th, 2022 (GTO Program #1198; PI: Stiavelli). The program aims to leverage the gravitational lensing effect of the deflector, which spatially stretches the quasar’s image, to study the relationship between the quasar and its host galaxy with greater rest frame angular
1 The updated analysis in H0LiCOW XII found only one of these three galaxies to be above the standard flexion shift cutoff of ∆3 x > 10−4 asec, but the two other galaxies were included to be safe. For a fair comparison, we also choose to include all 3 perturbing galaxies.
-
±
The spectroscopic observations are presented in via the ESO- MUSE integral field spectrograph H0LiCOW X, which we sum- marize here. The lens lies within a galaxy group at z = 0.6588
with a velocity dispersion of σ = 500 80km s−1 measured
from 22 member galaxies (Wilson et al. 2016; Momcheva et al.
2006, 2015). The velocity dispersion for the main deflector was measured to be σLOS = 250km s−1 with a total uncertainty of σσLOS = 19km s−1. Additional multi-object spectroscopic data was collected for constraints on galaxies in the vicinity (< 2′)
of the lens, using the ESO FORS (Appenzeller et al. 1998) and Gemini GMOS (Hook et al. 2004) instruments.
JWST NIRSpec data was also taken for WFI2033, which we plan on modeling for spatially resolved kinematics in future analyses, which will enable tight constraints for more general mass models of the deflector. As the focus of this paper is the power-law model, we use the higher resolution of NIRCam to this end.
-
Wide-Field Photometric Data and External Convergence
Measurements utilizing time-delay cosmography are sensitive to whether the system resides in a relatively over- or under-dense line-of-sight (defined as κext, see e.g. Eq. (15)), but these den-
sities are not directly observable. In response, κext is estimated
by matching observable tracers in the wide-field data to those
in simulations—a process that has undergone multiple develop- ments over time (Suyu & Halkola 2010; Fassnacht et al. 2011; Suyu et al. 2013; Greene et al. 2013; Wells et al. 2023, 2024a).
Therefore, wide-field photometric data is needed for pho- tometric redshift measurements, galaxy-star classification, and measurements of the stellar masses of all galaxies i < 23 mag within a 120 arcsec radius. These data are described in
H0LiCOW X, with the resulting values of κext computed in H0LiCOW XII, which we keep consistent for this analysis. While the cosmological predictions of WFI2033 will require up- dated sampling of κext2, this has a negligible impact on a one- to-one comparison of lens models due to the resulting model-
independence of κext. Because of this, any comparison of models must be done using the same sample of κext. In this work, our aim is to compare models of HST- and JWST-based data. So, for a fair comparison, we use the same sample of κext as Rusu et al. (2020) in our lens modeling and note any inference on cos-
mological parameters will require an update in the sampling of external convergence, which is left for future work.
-
-
Careful consideration must be taken to ensure the complex, ex- tended PSF structure does not bias measurements of our param- eters. Shajib et al. (2022) and Ding et al. (2021) demonstrated that, even for the same imaging data, models with different PSFs can have significant discrepancies on crucial model parameters, such as the power-law slope, γ, or the external shear, γext. As
this is the first cosmographic strong lens model with JWST’s
challenging PSF, we will discuss a new method to model the PSF, and compare it to the method adopted by recent studies to describe the HST PSF.
While the extent of the PSF may be similar to HST, Fig- ure Fig. 2 illustrates how the complexity of the PSF has in- creased. A crucial step in lens modeling is the separation of the source’s light into the extended component, contributing to the ring-like structure of the lens, and the point source compo- nent, contributing to the bright images seen at image positions A, B, C, and D. A model incorrectly disentangling these compo- nents can fit non-physical–or even non-existent–light structure,
2 Recent discussions have raised questions on whether the model- predicted external shear (γext, mod) should be included as a prior for solv- ing for κext, as it may not correlate one-to-one with the physical shear in these environments (e.g., Etherington et al. 2023). As this topic is
not the focus of this work, we decide to sample the same distribution of
potentially biasing the model. For example, a bias in the power- law slope, closely tied to the source light profile, will bias the final inference of H0.
-
To generate an initial PSF, we followed the standard procedure of extracting stars in the field as reference. Ideally, a selection of stars should be as close to the lens galaxy in angular separation as possible to minimize the effects of PSF distortion across the field. However, a larger constraint on star selection proved to be finding stars of similar flux to the lensed quasar images. Having similar flux is crucial to avoid discrepancies caused by saturation and non linear effects at the bright end and by noise at the faint end. Across the entire 2.5′ square field, 13 stars have desirable
flux, with three saturating the centers of their PSFs, and four
being on the edge of an exposure. Our selection is composed of the remaining 6 stars.
Previous strong lens modeling approaches typically utilize either the astropy3 core package photoutils (Astropy Col- laboration et al. 2022) or the more recently developed PSFr (Bir- rer et al. 2022). We will be comparing to PSFr, which imple- ments the following algorithm:
-
Get an initial guess of the PSF by stacking the cutout stars (without any recentering)
-
Calculate the sub-pixel centroid of each star’s cutout
-
Fit for the sub-pixel centroid of the PSF model
-
Line up the sub-pixel centroids of your PSF model and each of the cutout stars
-
Interpolate the PSF model “into the frame" of each star with sub-pixel interpolation (keeping their centroids lined up)
-
Calculate the residuals of each interpolated PSF model to its star’s cutout data
-
To combine the residuals from each of these different sub- pixel-interpolated frames, transform these residuals back “into the frame" of your PSF model (an “inverse sub-pixel shift")
-
With all the residuals in the PSF model’s frame, stack them (e.g., median) to inform a correction to your PSF model
-
Repeat steps (4)-(8), optionally repeating step (3)
In this work, we implement STARRED (Michalewicz et al. 2023; Millon et al. 2024) as an alternative method to generate an initial PSF for strong lens modeling. STARRED directly ad- dresses a few of the novel challenges that JWST proposes for strong lens modeling in particular.
First, wavelet regularization increases the population of suit- able stars for PSF generation by reducing the influence of vary- ing noise levels among the stars on the resulting PSF model. In- stead of representing our PSF in a spatial pixel basis, STARRED exploits the isotropic wavelet basis of starlets (Starck et al. 2006). Our PSF is sparsely represented in this starlet space, meaning it can be reconstructed with a linear combination of very few of these starlet bases. This is not true of the noise. As a result, one can simply select and add up the few coefficients required to reconstruct the PSF, dropping the low-power coeffi– cients composing noise. This denoising process is called image thresholding, and is one of the main benefits of wavelet regular- ization.
Second, sparse regularization has been shown to improve photometric accuracy and astrometric precision to better than a
κext—still informed by γext, mod—as H0LiCOW XII to guarantee a fair
comparison. 3 http://www.astropy.org
hundredth of a pixel for high signal-to-noise data (Fig. 4 of Mil- lon et al. 2024), the latter of which being crucial for tight con- straints on image positions and thus H0 (Birrer & Treu 2019).
~
We estimate uncertainty in the PSF based on residuals of our initial fitting of field stars. Residuals without the additional un- certainty due to the PSF are shown in Fig. C.1 and discussed in Appendix C. We additionally crop the initial PSF models based on a minimum threshold, removing the bottom 1% of the ini- tial kernel. This removes the extended background noise from the initial PSFs (due to the differing fluxes of the quasar vs the stars used in the initial fit). Without this additional cropping, the point-sources used in the modeling would have extended light profiles that add a nonphysical sheet of flux across the image.
-
-
Iterative PSF Updating
After acquiring our initial PSF models, we iteratively update each model based on the quasar images during the image for- ward modelin (result shown in Fig. 4). By first subtracting the best lens light and extended source surface brightness profiles, we can update the PSF with the additional four remaining quasar point-sources, which is repeated until the PSF model converges to a solution (or reaches 100 iterations). This entire process is re-
-
-
This section describes the modeling procedures used to fit the JWST NIRCam data for the final inference on lens parameters. For our analysis, we used Lenstronomy4, a python software package developed by Birrer & Amara (2018b) (Birrer et al. 2015, 2021). The accuracy of Lenstronomy has been verified through methods such as the Time-Delay Lens Modelling Chal- lenge, where “true parameters" of a simulated system were hid- den from different analysis teams (Ding et al. 2021). Since then, it continues seeing use as one of the primary modeling methods for time-delay cosmography today.
Priors on all parameters are uniform with values based on initial modeling of the system, unless otherwise stated, and are summarized in Table 1.
-
Main Deflector Galaxy D
The main deflector is a massive elliptical galaxy. Following stan- dard practice and to facilitate comparison with the previous anal- ysis (Rusu et al. 2020) we characterize its mass profile via the power-law elliptical mass distribution (PEMD; Barkana 1998) specified by
peated six times during the initial fitting process for each model, to avoid biasing the PSF model based on our initial parame-
κPEMD (θ1, θ2) ➟ 3 − γ , θE
γ−1
. (27)
PSFr-based iteration method for both initial PSF models. We use a completely asymmetric PSF, as we found the image residuals
1
2
ter values or overall model parameterization. We use the same
2
qmθ2 + θ2/qm
Free parameters are the logarithmic slope (γ), the Einstein radius (θE), and the axis ratio (qm). If the angle of the major axis is given by φm with reference to the RA-Dec frame, then we define the coordinates, (θ1, θ2), aligned with the major and minor axes.
For the deflector’s light profile, we use two (elliptical) Sérsic profiles,
10 2
10 4
10 6
0
10 6
10 4
10 2
Initial STARRED
PSF Difference (Row – Column)
,
1/n
qLθ2 + θ2/qL
Updated STARRED
I(θ1, θ2) = A exp−k
1 2
r − 1, (28)
eff
Initial PSFr
with the amplitude (A), the effective (half-light) radius (r
eff
), the
Updated PSFr
Initial STARRED
Updated STARRED
Initial PSFr
Updated PSFr
axis ratio (qL), the Sérsic index (n), and a corrective constant (k) (Sérsic 1968).
The two Sérsic light profiles struggle to fit the very center of the deflector, but this has little impact on our lens model. Since fitting the lens light serves two roles—providing a light model for the kinematic analysis and removing the lens light contribu- tion from the extended source’s ring structure—we assess how a mask would impact these goals. Preliminary tests showed that masking this region had little effect on model-predicted velocity dispersions, and the ring structure lies sufficiently far from the deflector’s center. Therefore, we mask the central 0′.′15 to pre-
vent this region from driving the fit (shown in Fig. 5).
We include an additional external shear term to account for distortion by the large-scale structure of the lens environment, such as the group’s halo. We further constrain the ellipticity and
Fig. 4. Comparison of the PSF kernels used. First we show the STARRED initial PSF kernel, followed by an example updated PSF from one of the STARRED models. Then we show the PSFr initial PSF kernel, again followed by an example updated PSF from a PSFr model.
position angle of the mass profile based on the observed ellip- ticity and position angle of the light profile following Schmidt et al. (2022). And, to mitigate degeneracies caused by the exter- nal shear, we apply a Gaussian prior on the ellipticity parameters.
4 https://github.com/sibirrer/lenstronomy
2
Observed
1
0
E
N
1"
2
Reconstructed
1
0
E
N
1"
5
Normalized Residuals
2
0
E
N
1"
log10 flux
log10 flux
(fmodel – fdata)/
2
1 1 5
0.5
Reconstructed Source
0.0
0.
1.
0.1"
2.0
Source Light
1.0
0.0
E
N
1"
10
Magnification Model
B
5
A2
A1
0
C
E
N
5
1"
log10 flux
log10 flux
det (A 1)
5
0
1.5 1.5 10
Fig. 5. An example model output from the best performing model, seen in Table B.1, though the other models are qualitatively indistinguishable. Top left shows the observed image, top middle shows the model-predicted reconstruction of the image, and top right shows the fit’s residuals normalized by the estimated uncertainty of each pixel (see Appendix C for more details). Bottom left is the reconstructed source plot, featuring more detail than any parametric source modeled for time-delay cosmography thus far. The star symbol denotes the location of the quasar host galaxy’s centroid. Bottom middle shows the lensed, unconvolved extended-source light, and bottom right shows the magnification map of the system.
-
Quasar and Host Galaxy
The quasar is modeled as point sources in the image plane con- volved with the PSF, whose generation is further detailed in Sect. 4. We constrain the image positions within an astrometric precision of 0.34 mas, well below the astrometric requirements for cosmography (Birrer & Treu 2019).
In addition to the point-source structure generated by the quasar, we also fit the extended light produced by the quasar’s host galaxy. To account for its complexity, we add standard lin- ear shapelet functions (Refregier 2003; Birrer et al. 2015) to our quasar host galaxy’s Sérsic light profile. The host galaxy’s pa- rameterization is explicitly left scale invariant, as forcing a fixed source reconstruction scale can assume an incorrect solution to the mass-sheet transformation, underestimating uncertainties in H0 (Birrer et al. 2016).
We find the optimal number of shapelet bases to reconstruct the source by finding the BIC turnover point, where the addi- tional complexity of more shapelets becomes incommensurate
-
Galaxy-Scale Perturbers
Of the nearby galaxy perturbers, H0LiCOW X determined three had a flexion shift large enough such that their direct inclusion in the lens model is necessary, due to potential higher-order effects. H0LiCOW XII redid the flexion shift analysis, finding only one of these perturbers (labeled G2) to be above the threshold, how- ever they included the other two perturbers (labeled G3 and G7) to be safe5. In addition, there is a small perturber near the main deflector, which is assumed to be a satellite at the same redshift due to lack of a spectrum, labeled X in Fig. 3.
G2 has the largest flexion shift of the galaxy perturbers, due to its close proximity and comparable size to the main deflector
−8
D. We also include a light profile for G2 due to its proximity, representing it with a single Sérsic. G2 also lies outside of the lens plane, at a redshift of zG2 = 0.7450 (while zD = 0.6575). We take this into account when converting its velocity disper- sion (σLOS, G2 = 223+14km s−1) into the Einstein radius assuming
with the goodness-of-fit as determined by the Bayes Information Criterion (BIC). All initial conditions are kept the same between these models, with the only distinction being the maximum poly- nomial order permitted by the model. Notably, we find that the largest polynomial order tested (nmax = 30) still recovers bet- ter BIC values than lower-order models; however, the compu- tational time is nonlinear with nmax. Because of this, we adopt
∈ { }
the significantly faster polynomial orders of nmax 18, 20, 22 , with a minor drop in accuracy of the source reconstruction.
that the lensing velocity dispersion is equal to the stellar velocity
5 In addition to galaxies, H0LiCOW X also measured and made pre- dictions on the impact of galaxy groups on the lens potential. While no group was found to have a flexion shift above the typical threshold of
∆3 x > 10−4, H0LiCOW XII also included two group halos in the model for extra caution . While these models did recover the best BIC values,
the improvement was within the typical BIC variance. Because of this, we do not explicitly include these groups in the model to minimize the number of nonlinear parameters and computational time.
dispersion within the scatter (Treu et al. 2006). Thus, effectively, we model the system with the multi-plane lens equation.
≈
G3 lies significantly far ( 7.2′′) from the lens as to not require modeling of its surface brightness profile, however G3
−19
may have a flexion shift large enough to require direct model- ing of the mass distribution. Therefore, we model its mass with an SIS profile. H0LiCOW X measured a velocity dispersion of
way that is inconsistent with their measured redshifts and veloc- ity dispersions.
To accurately account for G2, we must solve the multi-plane lens equation, which requires adopting a cosmology. Since we are concerned with comparing two lens models and not with do- ing cosmographic inference, this does not affect our conclusions.
We join the source-plane position of the quasar with the cen-
σLOS, G3
= 79+23km s−1, which we convert to a prior on the Ein-
ter of the Sérsic and its shapelets, the latter two comprising the
stein radius, with the same redshift as the deflector as their red- shifts are consistent within the range allowed by peculiar veloc- ities.
Similarly, G7 is also far enough away such that a surface brightness profile does not need to be included. However its gravitational effects may not be negligible, and therefore we model G7 with an SIS profile anchored to its measured veloc-
ity dispersion of σLOS, G7 = 166+7km s−1 at the redshift of the
host galaxy’s extended light. We also join the centers of the two Sérsics describing the bulge and disk components of D. For both the satellite and G2, we bind the center of their mass profiles with their light profiles.
The point-source positions of the quasar are constrained by the pixels of the JWST image, assuming the PSF model de- scribed in Sect. 4. This allows us to solve the lens equation to constrain six of the lens mass parameters, one for each of the in-
deflector. −6
-
In the previous HST power-law models, H0LiCOW XII found the mass of the satellite X to be unconstrained by the imaging
−0.002
data, with an SIS Einstein radius of θE = 0.018+0.002. In this work, we also fit an SIS profile centered on a Sérsic light profile, due to its proximity to the deflector. To get a rough, initial estimate
of the satellite’s Einstein radius, we assume it follows the Faber- Jackson relation (Faber & Jackson 1976) and estimate its radius based on a preliminary model’s radius of the deflector, yielding θE,X, init ≈ 0.15′′.
-
Flexion
Models for time-delay cosmography feature higher-order correc- tions of the lens potential for any mass affecting the system that is not explicitly included in our lens models. The second order corrections for an external mass’ affect on the lensing potential is given by the external convergence correction (Sect. 3.4), which changes the scale/magnification of the lensed image, and the ex- ternal shear term (Sect. 7.2), which results in a stretching of the lensed image.
One of the questions for higher spatial resolution data is whether these typical corrective terms used in the past are suf- ficient for this data, or if the increase in resolution requires ad- ditional higher-order corrections. After the second-order shear and convergence terms, the next term would be the third-order correction called flexion (Goldberg & Bacon 2005). To test this, one of our systematic choices is the addition of a free flexion term, which allows the model to additionally vary a triangu- lar bending effect on the lensed light. This additional term has only been applied to one time-delay cosmography system thus far (J1721+8842, Schmidt et al., in-prep).
-
Additional Constraints
The velocity dispersions of the galaxy-scale perturbers (G2, G3, and G7) are used to estimate their Einstein radii. The ratios of these radii are fixed and scaled together using a single scal- ing parameter. To encode maximum information, we convert the lower/upper 1σ estimates of the velocity dispersions into our
1σ estimates in θE. This mitigates potential degeneracies aris-
ing from the perturbers’ positions (and the additional external
shear), and prevents the model from optimizing the radii in a
dependent, relative positions of the 4 images. This is done with the PROFILE_SHEAR solver in Lenstronomy, which solves for the following non-linear parameters of D: its centroid (RA and Dec), axis ratio and position angle (represented as e1 and e2), Einstein radius (θE), and the angle of external shear (ϕext).
-
Image Forward Modeling
Given a proposed set of these parameters, the linear response functions in the data (such as amplitude parameters, point sources, and coefficients of shapelets) are rendered in the image plane and optimized with a linear minimization method based on the imaging likelihood.
We supersample the central regions of the point sources by a factor of 3 for ray-tracing and perform the PSF convolution on this subset of the data using the adaptive compute mode. This provides improved constraints on the image positions, as well as avoiding aliasing and other artifacts resulting from modeling sub-pixel features in the PSF. We do not apply the PSF convo- lution on the entire supersampled image, as this quickly domi- nates the computational cost. To benefit from the supersampled grid, the PSF kernel must also be supersampled at a scale of 3. Without a supersampled PSF, the sub-pixel scale features of the PSF would not be accurately accounted for in the convolution, leading to poor fitting of the image positions. By confining this process to only the regions nearest the PSFs, we still gain the increase in accuracy without suffering the computational cost.
Our uncertainty in the data is given by the 2-D resampled Poisson, readout noise, and flat-field variance estimates summed in quadrature, which we then scale by the weight map to account for the drizzling procedure. We add additional uncertainty in the data to account for uncertainty in the PSF models, discussed in Appendix C.
-
Likelihood Sampling
When computing the image likelihood, we punish models that do not map the multiple images match back to the same posi- tion in the source plane with a 0′.′004 tolerance level. We also punish models with negative Sérsic or point-source amplitudes. We only compute the image likelihood on a subset of the data,
as light at the center of the main deflector or far outside the ring contains little information as long as the light profile of the main deflector is accurately represented near the ring itself. The MASK systematic test (in Sect. 6.1) extends the size of the outer radial mask, to ensure this choice does not bias our results.
Table 1. Model parameter priors used, with some values rounded.
Object Component Parameter Distribution Initial Position Step Size (1σ)
–
–
–
θE, PEMD N(0.879, 0.133) 0.879 0.133
Main Deflector (D)
z = 0.6575
External Shear
z = 0.6575
Satellite (X)
z = 0.6575
Galaxy G2
z = 0.7450
Galaxy G3
z = 0.6575
Galaxy G7
z = 0.6575
Free Flexion
z = 0.6575
Quasar Images
z = 0.6575
Quasar Host Galaxy
z = 1.662
Mass:
PEMD
Light:
Bulge Sérsic
Light:
Disk Sérsic
Mass:
Shear
Mass:
SIS
Light:
Sérsic
Mass:
SIS
Light:
Sérsic
Mass:
SIS
Mass:
SIS
Mass:
Flexion
Light:
Point Sources
Light:
Sérsic
Light:
Shapelets
γPEMD U(1.5, 2.5) 1.95 0.02
(e1, e2) N(0, 0.2) (-0.047, 0.107) 0.2
(x, y) U(−3, 3) (-0.031, 0.021) 0.003
Rsersic U(0.001, 5) 0.296 0.01
nsersic Fixed 4
(e1, e2) U(−0.2, 0.2) (-0.048, 0.089) 0.1
(x, y) U(−1, 1) (-0.009, 0.009) 0.003
Rsersic U(0.001, 5) 1.752 0.05
nsersic Fixed 1
(e1, e2) U( 0.5, 0.5) (-0.164, 0.154) 0.1
(x, y) Joined with Bulge Sérsic
γext, 1 U( 0.5, 0.5) 0.109 0.1
γext, 2 U( 0.5, 0.5) 0.025 0.1
(x, y) Fixed 0
θE, X U(0.0001, 0.25) 0.15 0.01
γext, 2 U( 0.5, 0.5) 0.025 0.1
(x, y) Joined with Light
Rsersic U(0.0001, 0.1) 0.044 0.01
nsersic U(1, 10) 2 0.2
(e1, e2) Fixed (0.015, 0.018)
(x, y) U(Initial ± 0.07) (0.228, 2.046) 0.003
θE, G2 N(0.622, 0.062) 0.622 0.1
(x, y) Joined with Light
Rsersic U(0.001, 5) 1 0.05
nsersic U(0.1, 2) 0.5 0.1
–
(e1, e2) Fixed (0, 0)
(x, y) Fixed (-3.998, -0.034)
θE, G3 N(0.088, 0.048)†
(x, y) Fixed (2.823, 7.006)
θE, G7 N(0.388, 0.030)†
(x, y) Fixed (-4.578, -11.719)
–
(g1, g2, g3, g4) U( 0.1, 0.1) 0 0.01
(x, y) Fixed (0, 0)
A1(x, y) U(Initial ± 0.3) (-0.760, 0.958) 0.007
A2(x, y) U(Initial ± 0.3) (-0.050, 1.076) 0.007
B(x, y) U(Initial ± 0.3) (1.434, -0.302) 0.007
C(x, y) U(Initial ± 0.3) (-0.686, -0.579) 0.007
Rsersic U(0.04, 0.3) 0.2 0.01
nsersic U(0.15, 4) 1 0.25
(e1, e2) U(−0.15, 0.15) (0, 0) 0.03
(x, y) U(Initial ± 1) (-0.45, 0.1) 0.03
β U(0.031, 0.3) 0.031 0.025
{ }
nmax Fixed 18, 20, 22 0.1
(x, y) Joined with Sérsic
Notes. Parameters grouped by components (Deflector, Sérsic 1, Shapelets) with their corresponding types (Mass or Light). The bounds are represented as uniform distributions, with initial guesses and sampling step sizes (σ).
† While the Einstein radii are given priors, their actual values are scaled with G2. This ensures the constraints provided by their observed velocity dispersions are passed to the models, while also minimizing the degeneracies between their parameters.
We introduce a likelihood penalty for models where the main deflector’s mass profile is significantly more elliptical than, or
misaligned with, its light profile, inspired by the approach of Schmidt et al. 2022. Without this additional constraint, we find
the large number of perturbers lead to a nonphysical degeneracy in the deflector’s mass ellipticity. We add a flat prior to ensure the axis ratio of the mass profile, qL, agrees within 0.2 of the lens-light’s axis ratio, qLL, ensuring that
qLL − qL < 0.2. (29)
Similarly, we enforce a constraint to ensure that the position an- gles of the mass profile, ϕL, and the lens-light profile, ϕLL, are aligned. Specifically, we penalize configurations where the angu- lar difference between these two position angles exceeds a cer- tain threshold, which is a function of the mass profile’s axis ratio, qL. The alignment criterion is defined as:
180
-
To ensure our parameterization does not bias our results, we give each model a weight based on its ability to predict the data, and then combine these weighted models to form the final inference. As expected, we observe significant variability in the models’ ability to reconstruct the imaging data. Further, goodness-of-fit statistics on imaging data alone cannot distinguish models with diverging physical interpretations, as profiles can become mass- sheet transforms of each other (Sect. 2.2). To address this, mod- els are assigned two separate weights that are combined for the final inference: the first reflecting their ability to reconstruct the imaging data, and the second their ability to predict the deflec-
∆PA = |ϕL − ϕLL| ×
, (30)
π
tor’s velocity dispersion. The models are then combined to form a final weighted prediction. The resulting model weights can be
where ∆PA is the difference in position angles in degrees. Con-
figurations are rejected if
5
∆PA > 10 − qL − 1.0001 . (31)
This condition ensures that the mass and light profiles remain well-aligned, particularly for more elliptical systems where mis- alignment is less tolerable. A small offset is introduced to prevent numerical instabilities in the computation.
-
-
-
This section describes how the systematic tests were combined to form the final inference on the Fermat potential. We first re- view the systematic tests considered for the final inference in
-
Imaging-Based Weighting
To prevent overly-complex models from being unjustly favored, we employ the Bayesian information criterion (BIC) to statisti- cally weight our 24 models. Compared to similar criteria such as the Akaike information criterion (AIC), the BIC is more pe- nalizing to model complexity, encouraging more parsimonious solutions. Formally,
BIC = ln(n) k − 2 ln Lˆ , (32)
where n is the number of data points, k is the number of model
Sect. 6.1 (a complete list of preliminary tests can be found in
parameters, and Lˆ is the maximum likelihood of the model.
Appendix A). Then, we discuss how these models were weighted and combined in Sect. 6.2.
-
Here, we summarize the model combinations introduced in Sect. 5:
-
Two choices for the the inclusion of free flexion,
-
No Flexion
Within a 2.1′′ mask, n = 14, 592 pixels, and for 2.3′′ it is n = 15, 930 pixels. The number of nonlinear parameters is ei- ther 32 or 36, based on choice of Flexion, while the number of
≫
linear parameter ranges from 199 to 285. The BIC remains a valid criterion here as n k.
Formally, the BIC comparison is only valid for models fit to the same data set. One of our systematics is mask size, which inherently changes the size of the data. As an approximation, following Wong et al. (2017) and Birrer et al. (2019), we use the
same mask (the smaller of the two) to calculate χ2 and the BIC,
– Flexion
-
-
Three choices with increasing flexibility in source surface
-
SERSIC+18nmax
allowing consistent model comparisons despite slight differences in the masked region.
When comparing any two models M1 and M2 with BIC val-
-
SERSIC+20nmax
-
SERSIC+22nmax
-
-
Two different pixel masks for the imaging likelihood
ues BIC1 and BIC2, their relative probabilities can be written as
p(M1) " BIC1 − BIC2 #
-
MASK 2.1” (or 68 pixel) radius
-
MASK 2.2” (or 71 pixel) radius
-
-
Two different initial PSF generation methods
p(M2) ∝ exp − 2 . (33)
-
STARRED PSF
-
PSFr PSF
-
This gives 24 separate model configurations—and thus 24 separate Fermat potential difference posteriors—which are com-
Defining BICmin as the lowest BIC among all models, we weight each model n by
1 if x ≤ BICmin,
optimize the lens model, we first use a Particle Swarm Optimiza-
bined through a weighted sampling, discussed in Sect. 6.2. To
fBIC(x) = exp .− x−BICmin . if x > BIC
(34)
2
min
,
tion (PSO) algorithm to identify high-likelihood regions, then re-
fine those solutions with a Markov Chain Monte Carlo (MCMC) sampler to fully map out the posterior distributions. We also it- eratively refine the PSF for each model between the PSO steps, making it more robust to any bias introduced by model choices. Each model takes approximately one week to run on 36 cores.
6 We verify our choices are above the minimum data-supported thresh- old in Appendix D.
where x is the model’s BIC. This piecewise definition penalizes
large BIC values.
However, following Birrer et al. (2019) and TDCOSMO IX, we additionally convolve each model’s nominal BICn with a Gaussian of width σintrinsic, reflecting our uncertainty in BIC
due to limited sampling and systematics. We estimate σintrinsic
by comparing BIC values across our shapelet tests (nmax =
18, 20, 22), fixing the other tests (No Flexion PSF, STARRED
PSF, MASK 2.1”). That convolution ensures that small fluctua- tions in BICn do not spuriously dominate the weighting. Specif- ically, each model (n) receives a single BIC-based weight:
We find models including free flexion introduce a new de- generacy between the Fermat potential and the Einstein radius of the main deflector, which the non-flexion models do not have. This implies the additional degrees of freedom granted by free
√
2 σ
Wabs,n =
−∞
2π σintrinsic
exp −
2
intrinsic
fBIC(x) dx. (35)
, ∞ 1 .BICn − x 2
Flexion are not constrained by the data alone, and so we discard
flexion models for the final inference. Similar to H0LiCOW XII, we also find the kinematic weighting had little effect on the final
The final relative weight Wn is then normalized so that max(Wn) = 1.
Because each model has exactly one BIC-based weight, all individual MCMC draws from that model share that weight. Models with lower (better) BIC values receive proportionally more draws in the final aggregated parameter space. We find that resampling the MCMC chain for each model according to Wn offers an efficient combination of our systematically differ- ent lens models, ensuring each model is represented fairly in the final posterior distribution.
Parameter
This work (JWST)
H0LiCOW XII‡ (HST)
∆τBC
0.469+0.007
−0.009
0.455+0.010
−0.009
∆τBA1
0.311+0.005
−0.007
0.298+0.010
−0.007
∆τBA2
0.330+0.006
−0.008
0.31+0.04
−0.02
ΘE,PEMD
0.977+0.002
−0.002
0.93+0.01
−0.02
γPEMD
1.86+0.01
−0.01
1.95+0.02
−0.01
qPEMD
0.774+0.010
−0.007
0.79+0.01
−0.01
ϕPEMD
30.3+0.7
−0.5
33.0+0.8
−0.9
γext
0.088+0.003
−0.003
0.112+0.006
−0.004
ϕext
95+5
−5
84+7
−2
ΘE,X
0.100+0.005
−0.005
0.0009+0.0010
−0.0007
ΘE,G2
0.78+0.02
−0.02
0.93+0.03‡‡
−0.06
σV
236+20
−9
232+10
−7
As our models are a relatively sparse sample of all possible configurations, the BIC-based weights help mitigate overfitting tendencies while still including sufficiently flexible models. For
inference of the Fermat potential. A corner plot of our posteri- ors is shown in Fig. 7, with final parameter estimates reported in Table 2.
Power-Law Model Parameters
Table 2. BIC+Kinematic-weighted posterior parameter estimates.
parameters such as the Fermat potential differences, our sam-
pling is dense enough that the weighted posterior remains robust. For the broader BIC distribution itself, the convolution integral accounts for scatter, ensuring no single model with an artificially low BIC is overemphasized.
-
-
-
Kinematics-Based Weighting
.
Our second weight assesses consistency with the observed ve- locity dispersion of the main deflector. For each MCMC sample, we compute the model-predicted velocity dispersion, σap, model, then assign a kinematic weight based on a Gaussian likelihood centered on the observed velocity dispersion, σobs,D:
fkin = N σap, model | σobs, D, σobs,σ , (36)
taking into account the uncertainty in the observed velocity dis- persion, σobs,σ. This effectively penalizes models that predict
Notes. Reported values are medians with errors corresponding to the 16th and 84th percentiles.
‡ We note our weighting process may differ slightly for rounding discrepancies to values reported in H0LiCOW XII.
significantly different σap, model from the measured dispersion.
We then combine the kinematic weight with the BIC-based
imaging weight to form a total weight:
Wtotal = WBIC × fkin. (37)
After computing Wtotal for each model sample, we normalize it such that max(Wtotal) = 1. For the final posterior combination in the corner plots, each MCMC sample is drawn in proportion to Wtotal. For the final parameter estimates, we instead calculate the weighted median and 16th/84th percentiles.
-
-
We recover a 22% tighter estimate of the Fermat potential differ- ence after incorporating both kinematic and imaging constraints (Fig. 6). While our results remain consistent with previous mod-
‡‡ The GLEE parameterization used an SIE model for G2, while we fit an SIS model due to degeneracies with the other per- turbers.
-
The PEMD models in H0LiCOW XII found satellite X to have mass consistent with zero, indicative of poor constraints from the HST-based data. In addition, their models found an offset between the mass and light centroids for the main deflector of
– ∼
0.02 0.03′′ ( 200 pc, assuming flat ΛCDM cosmology with h = 0.7 and Ωm = 0.3). The influence of satellite X (lying to the north-west) was proposed as a partial potential explanation, as
the mass’ centroid was found to be south-east of the light.
−0.005
Our models place the Einstein radius of satellite X at θE, X =
els, our measured 3.1% increase in the Fermat potential suggests
0.100+0.005, far closer to the theoretical value of θ
E,X, init
≈ 0.15,
a corresponding 3.1% increase in this lens’ inferred value of H0, assuming all other parameters remain unchanged. This adjust- ment would shift the median measured H0 for this lens from
71.6 to 73.8 kms−1 Mpc−1, with uncertainties likely similar to the HST result. A comprehensive cosmographic analysis, including a resampling of the external convergence (κext), will be presented in a forthcoming paper, which could further adjust the overall H0 estimate in either direction.
estimated in Sect. 5.4. However, our models still identify the off– set of the deflector’s mass, though smaller in scale.
-
The previous HST-based models were unable to explain the large shears required to model the system, even after including two large galaxy groups along the LOS (which we do not explicitly
_M_o_d_e_l_W__e_i_g_h_t_in_g
Unweighted
Kinematics
3.73%
4.15%
3.29%
3.48%
Imaging
Imaging
+
Kinematics
2.16%
2.23%
1.74%
1.73%
+19.44% Precision
Rusu et al. 2020 (HST)
This Work (JWST)
-
|
|
+22.42% Precision
0.44 0.46 0.48 0.50 0.52
Fermat Potential Difference
BC(arcsec2)
(Numbers in % represent fractional uncertainty)
Fig. 6. Comparison of the Fermat Potential differences of images B and C between the HST models (Rusu et al. 2020, orange) and this work (blue). The fractional uncertainty (measured as the standard deviation over the median) is given to the right of the measurements. Vertically, we display the different weighting schemes, starting with the unweighted models at the top. In the second row we have the kinematic weighting, which preferentially weights both estimates slightly higher and wider. This is due to the model-predicted velocity dispersions being lower than the spectroscopically-observed velocity dispersion measured in H0LiCOW X. However, the velocity dispersion estimates are consistent across the HST and JWST models (see Fig. 7). The imaging weighting in the third row is far more constraining than the kinematics for this system, and
~
it indicates the best-fitting models tend to be on the lower end of the models tested. This holds for both HST and JWST results. The estimate recovered by our best-fit models is 3.1% higher, most likely due to tighter constraints on the satellite’s mass provided by the NIRCam imaging. The last row combines both kinematics and imaging weights for the final estimate, with a 22% increase in precision compared to HST imaging.
model here, due to their minor impacts on the HST models, but will consider for the following cosmographic analysis paper).
In this work, leveraging the improved resolution and sen- sitivity of JWST imaging, we identify and robustly constrain the satellite galaxy X, which was only marginally detected in HST data. By constraining this satellite, we prevent the lens model from “absorbing” unmodeled mass distributions into a
large external shear term. As a result, our models prefer a sig- nificantly smaller external shear, γext = 0.088, approximately 21% lower than the values inferred by Rusu et al. (2020)’s HST- based power-law models. This value is in close agreement with the observed shear estimates from Wong et al. (2010)
-
Astrometry
In Fig. 8, we compare the astrometric precision of our mod- els to the H0LiCOW XII results. The expected relative posi- tions of the images, centered at (0, 0), are derived from HST observations with the brightest image (A1) serving as the ref- erence point. Compared to these positions, our models achieve significantly smaller discrepancies than those constrained by the previous HST-based models. Specifically, we find deviations of
3.8 mas for the distance between images B-A1, 3.1 mas for A2- A1, and 4.8 mas for C-A1. These are of the order of a tenth of a NIRCam pixel.
The overall relative astrometric uncertainty across all JWST models is 0.34 mas, obtained by adding the uncertainties of in- dividual images in quadrature. This is approximately 1/100 of a NIRCam pixel. Our estimated errors do not include poten-
tial residual systematics in the distortion correction of NIRCam, which should be in any case very much sub-pixel, especially over a range of a few arcseconds covered by our target.
In this work, we aimed to extend the methods of time-delay cosmography to the first JWST imaging of quadruply-imaged quasars; discovering the benefits offered by the improved imag- ing while exploring potential systematic effects across tele- scopes. We modeled the system in the NIRCam F115W band, utilizing the spectroscopic data from Sluse et al. (2019), and ap- plied the same wide-field data (external convergence) as Rusu et al. 2020 to enable a fair, one-to-one comparison with the pre- vious HST models. To accurately model the system, we imple- mented new tools to model the PSF of the system, via STARRED. Our resulting Fermat potential difference is in agreement with
the previous HST models at the 1σ level, well within expected uncertainties. The primary source of the discrepancy was the satellite, which was previously unconstrained in HST imaging.
|
|
|
|
E, PEMD
0.95
|
0.90
|
|
|
|
|
PEMD
2.0
1.9
|
|
|
|
|
|
1.8
|
Unweighted Weighted
(HST: Rusu et al. 2020)
|
qPEMD
0.80
0.75
|
|
|
|
|
PEMD
35
|
30
|
|
|
|
|
|
0.12
ext
|
0.10
|
|
|
|
|
|
0.08
|
|
ext
100
80
|
|
|
|
|
|
|
|
|
|
E, X
0.10
|
0.05
0.00
|
|
|
|
|
|
|
|
|
E, G2
1.0
|
0.8
|
|
|
|
|
|
|
|
|
|
|
|
V
250
200
|
0.45 0.50
BC
0.90 0.95
|
E, PEMD
1.8 2.0
|
PEMD
0.75 0.80
qPEMD
30 35
|
PEMD
0.075 0.100 0.125
|
ext
80 100
|
ext
0.0 0.1
|
E, X
0.8 1.0
|
E, G2
200 250
|
V
Fig. 7. Parameter distributions from our lens model results, compared to H0LiCOW XII (HST-based data) in gray. Unweighted samples are shown in pink, with the final BIC+Kinematic-weighted selection is shown in green. The contours represent the 68.3%, and 99.7% quantiles.
With JWST data, our models recover a value in agreement with theoretical mass-to-light expectations. In summary:
-
Modeled JWST’s PSF with STARRED, cutting 71% of the astrometric uncertainty compared to typical PSF-modeling techniques (Fig. 8) and enabling 17% smaller scale of source reconstructions (Appendix D)
-
Enhanced accuracy in the Fermat potential, driven by im- proved constraints on the satellite’s mass (Sect. 7.1). This also provides a direct observational demonstration of how PEMD models with external shear can overcompensate for small-scale effects—below the resolution limit of the data—by adjusting the external shear (Etherington et al. 2023).
We note this analysis was not done blinded to the model pa- rameters in order to verify our handling of the PSF. However, the analysis was effectively blinded to H0 as the system requires an updated sampling of the external convergence (κext), a key in-
Astrometry Comparison: JWST (STARRED vs PSFr) vs HST
|
|
|
|
|
|
|
|
|
|
PSFr STARRED |
B-A1 B-A1 |
PSFr STARRED |
A2-A1 A2-A1 |
PSFr STARRED |
C-A1 C-A1 |
|
|
H0LiCOW XII Predicted |
B-A1 |
H0LiCOW XII Predicted |
A2-A1 |
H0LiCOW XII Predicted |
C-A1 |
|
|
H0LiCOW XII Observed |
B-A1 |
H0LiCOW XII Observed |
A2-A1 |
H0LiCOW XII Observed |
C-A1 |
JWST [
|
|
|
|
HST [
20
15
10
ΔDec (mas)
5
0
−5
−10
−15
−20 −15 −10 −5 0 5 10 15
ΔRA (mas)
Fig. 8. Relative image position discrepancies, centered based on the GLEE observed image positions. All distances are measured relative to the brightest image, A1. We note Lenstronomy requires the observed and predicted image positions to align. Image A2 also lies closest to the satellite, which is unconstrained in HST data, likely causing their largest discrepancy in image position. Our models demonstrate astrometry well within the requirements for time-delay cosmography (Birrer & Treu 2019), offering much tighter constraints than the previous HST-based GLEE models.
gredient in the final inference on H0. In our next paper, we will make the complete inference on H0 for this system, along with two other quadruply-imaged quasars imaged with JWST.
For other future work, this analysis expands the sample of time-delay cosmography lenses to systems previously too faint to provide constraints on H0. This is extremely complementary to the explosion of new lenses from future surveys like LSST and Euclid, which will increase our sample of strong lens systems by two orders of magnitude (Collett 2015). Furthermore, develop-
ing methods to reliably constrain the PSF in fields with limited field stars will be critical, as the six suitable reference stars used in this field is unusually large.
Acknowledgements. DW and TT acknowledge support by NSF through grants NSF-AST-1906976 and NSF-AST-1836016, and from the Moore Foundation through grant 8548. 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 Uni- versities for Research in Astronomy, Inc., for NASA, under contract NAS5- 26555. This work is also supported by JSPS KAKENHI Grant Numbers
JP24K07089, JP24H00221. MS is partly supported by NASA through grant 80NSSC21K1294.
References
Abdalla, E., Abellán, G. F., Aboubrahim, A., et al. 2022, arXiv e-prints, arXiv:2203.06142
Appenzeller, I., Fricke, K., Fürtig, W., et al. 1998, The Messenger, 94, 1, aDS Bibcode: 1998Msngr..94. 1A
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, The Astrophysical Journal, 724, 511, aDS Bibcode: 2010ApJ…724..511A
Barkana, R. 1998, The Astrophysical Journal, 502, 531, arXiv:astro-ph/9802002 Bergé, J., Massey, R., Baghi, Q., & Touboul, P. 2019, Monthly Notices of the
Royal Astronomical Society, 486, 544
Birrer, S. 2021, The Astrophysical Journal, 919, 38, aDS Bibcode: 2021ApJ…919. 38B
Birrer, S. & Amara, A. 2018a, Physics of the Dark Universe, 22, 189
Birrer, S. & Amara, A. 2018b, Physics of the Dark Universe, 22, 189, arXiv:1803.09746 [astro-ph]
Birrer, S., Amara, A., & Refregier, A. 2015, The Astrophysical Journal, 813, 102, aDS Bibcode: 2015ApJ…813..102B
Birrer, S., Amara, A., & Refregier, A. 2016, Journal of Cosmology and Astropar- ticle Physics, 2016, 020, aDS Bibcode: 2016JCAP…08..020B
Birrer, S., Bhamre, V., Nierenberg, A., Yang, L., & Van de Vyvere, L. 2022, Astrophysics Source Code Library, ascl:2210.005, aDS Bibcode: 2022ascl.soft10005B
Birrer, S., Shajib, A. J., Galan, A., et al. 2020, Astronomy & Astrophysics, 643, A165, arXiv:2007.02941 [astro-ph]
Birrer, S., Shajib, A. J., Gilman, D., et al. 2021, Journal of Open Source Software, 6, 3283
Birrer, S. & Treu, T. 2019, Monthly Notices of the Royal Astronomical Society, 489, 2097
Birrer, S., Treu, T., Rusu, C. E., et al. 2019, Monthly Notices of the Royal Astro- nomical Society, 484, 4726
Blandford, R. & Narayan, R. 1986, The Astrophysical Journal, 310, 568, aDS Bibcode: 1986ApJ…310..568B
Bonvin, V., Millon, M., Chan, J. H.-H., et al. 2019, Astronomy & Astrophysics, 629, A97
Chen, G. C.-F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS[1907.02533]
Collett, T. E. 2015, The Astrophysical Journal, 811, 20, arXiv:1507.02657 [astro- ph]
Collett, T. E. & Auger, M. W. 2014, Monthly Notices of the Royal Astronomical Society, 443, 969, aDS Bibcode: 2014MNRAS.443..969C
Courbin, F., Chantry, V., Revaz, Y., et al. 2011, Astronomy & Astrophysics, 536, A53, publisher: EDP Sciences
Dainotti, M. G., Simone, B. D., Schiavone, T., et al. 2021, The Astrophysical Journal, 912, 150, publisher: The American Astronomical Society
Ding, X., Treu, T., Birrer, S., et al. 2021, Monthly Notices of the Royal Astro- nomical Society, 503, 1096, arXiv:2006.08619 [astro-ph]
Efstathiou, G. 2021, Monthly Notices of the Royal Astronomical Society, 505, 3866, publisher: OUP ADS Bibcode: 2021MNRAS.505.3866E
Etherington, A., Nightingale, J. W., Massey, R., et al. 2023, Strong gravitational lensing’s ‘external shear’ is not shear, arXiv:2301.05244 [astro-ph]
Faber, S. M. & Jackson, R. E. 1976, ApJ, 204, 668
Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, The Astrophysical Journal, 289, L1, aDS Bibcode: 1985ApJ…289L. 1F
Fassnacht, C. D., Gal, R. R., Lubin, L. M., et al. 2006, The Astrophysical Journal, 642, 30, arXiv:astro-ph/0510728
Fassnacht, C. D., Koopmans, L. V. E., & Wong, K. C. 2011, MNRAS, 410, 2167 Fassnacht, C. D., Pearson, T. J., Readhead, A. C. S., et al. 1999, The Astrophys-
ical Journal, 527, 498, aDS Bibcode: 1999ApJ…527..498F
Fassnacht, C. D., Xanthopoulos, E., Koopmans, L. V. E., & Rusin, D. 2002, The Astrophysical Journal, 581, 823, publisher: IOP Publishing
Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, The Astrophysical Journal, 882, 34, arXiv:1907.05922 [astro-ph]
Freedman, W. L., Madore, B. F., Hoyt, T., et al. 2020, The Astrophysical Journal, 891, 57, publisher: The American Astronomical Society
Gilman, D., Birrer, S., & Treu, T. 2020, Astronomy & Astrophysics, 642, A194 Goldberg, D. M. & Bacon, D. J. 2005, ApJ, 619, 741
Greene, Z. S., Suyu, S. H., Treu, T., et al. 2013, ApJ, 768, 39
Hook, I. M., Jørgensen, I., Allington-Smith, J. R., et al. 2004, Publications of the Astronomical Society of the Pacific, 116, 425, publisher: IOP Publishing Keeton, C. R. 2003, The Astrophysical Journal, 584, 664, publisher: IOP Pub-
lishing
Keeton, C. R. & Kochanek, C. S. 1997, The Astrophysical Journal, 487, 42, aDS Bibcode: 1997ApJ…487…42K
Knabel, S., Treu, T., Cappellari, M., et al. 2024, Spatially Resolved Kine- matics of SLACS Lens Galaxies. I: Data and Kinematic Classification, arXiv:2409.10631 [astro-ph]
Knox, L. & Millea, M. 2020, Phys. Rev. D, 101, 043533
Kochanek, C. S. 2002, The Astrophysical Journal, 578, 25, aDS Bibcode: 2002ApJ…578…25K
Kochanek, C. S., Morgan, N. D., Falco, E. E., et al. 2006, The Astrophysical Journal, 640, 47, publisher: IOP Publishing
Koopmans, L. V. E., Treu, T., Fassnacht, C. D., Blandford, R. D., & Surpi, G. 2003, The Astrophysical Journal, 599, 70, aDS Bibcode: 2003ApJ…599…70K Kovner, I. 1987, The Astrophysical Journal, 316, 52, aDS Bibcode:
1987ApJ…316…52K
Mamon, G. A. & Lokas, E. L. 2005, Monthly Notices of the Royal Astronomical Society, 363, 705, arXiv:astro-ph/0405491
McCully, C., Keeton, C. R., Wong, K. C., & Zabludoff, A. I. 2014, Monthly Notices of the Royal Astronomical Society, 443, 3631
McCully, C., Keeton, C. R., Wong, K. C., & Zabludoff, A. I. 2017, The Astro-
physical Journal, 836, 141, arXiv:1601.05417 [astro-ph]
Merritt, D. 1985a, Monthly Notices of the Royal Astronomical Society, 214, 25P, aDS Bibcode: 1985MNRAS.214P..25M
Merritt, D. 1985b, The Astronomical Journal, 90, 1027, aDS Bibcode: 1985AJ. 90.1027M
Michalewicz, K., Millon, M., Dux, F., & Courbin, F. 2023, Journal of Open Source Software, 8, 5340
Millon, M., Courbin, F., Bonvin, V., et al. 2020, A&A, 640, A105
Millon, M., Michalewicz, K., Dux, F., Courbin, F., & Marshall, P. J. 2024, Image deconvolution and PSF reconstruction with STARRED: a wavelet-based two- channel method optimized for light curve extraction
Momcheva, I., Williams, K., Keeton, C., & Zabludoff, A. 2006, The Astrophysi- cal Journal, 641, 169, publisher: IOP Publishing
Momcheva, I. G., Williams, K. A., Cool, R. J., Keeton, C. R., & Zabludoff, A. I.
2015, The Astrophysical Journal Supplement Series, 219, 29, publisher: The American Astronomical Society
Morgan, N. D., Caldwell, J. A. R., Schechter, P. L., et al. 2004, The Astronomical Journal, 127, 2617, arXiv:astro-ph/0312478
Mortsell, E., Goobar, A., Johansson, J., & Dhawan, S. 2022, The Astrophysical Journal, 933, 212, arXiv:2105.11461 [astro-ph, physics:gr-qc]
Oguri, M. 2007, The Astrophysical Journal, 660, 1, aDS Bibcode: 2007ApJ…660. 1O
Oguri, M. & Marshall, P. J. 2010, MNRAS, 405, 2579
Osipkov, L. P. 1979, Soviet Astronomy Letters, 5, 42, aDS Bibcode: 1979SvAL….5. 42O
Refregier, A. 2003, Monthly Notices of the Royal Astronomical Soci- ety, 338, 35, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1046/j.1365- 8711.2003.05901.x
Refsdal, S. 1964, Monthly Notices of the Royal Astronomical Society, 128, 307, aDS Bibcode: 1964MNRAS.128..307R
Riess, A. G., Casertano, S., Yuan, W., et al. 2021, The Astrophysical Journal, 908, L6, aDS Bibcode: 2021ApJ…908L. 6R
Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, The Astrophysical Journal, 876, 85, arXiv:1903.07603 [astro-ph]
Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, The Astrophysical Journal Letters, 934, L7, arXiv:2112.04510 [astro-ph]
Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 1440, arXiv:1905.09338 [astro-ph]
Saha, P. 2000, The Astronomical Journal, 120, 1654, aDS Bibcode: 2000AJ. 120.1654S
Saha, P. & Williams, L. L. R. 2006, The Astrophysical Journal, 653, 936, aDS Bibcode: 2006ApJ…653..936S
Schechter, P. L., Bailyn, C. D., Barr, R., et al. 1997, The Astrophysical Journal, 475, L85, publisher: IOP Publishing
Schmidt, T., Treu, T., Birrer, S., et al. 2022, arXiv e-prints, arXiv:2206.04696 Schmidt, T., Treu, T., Birrer, S., et al. 2022, Monthly Notices of the Royal As-
tronomical Society, 518, 1260, arXiv:2206.04696 [astro-ph]
Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses, publication Title: Gravitational Lenses ADS Bibcode: 1992grle.book. S
Schneider, P. & Sluse, D. 2013, Astronomy and Astrophysics, 559, A37, aDS Bibcode: 2013A&A…559A..37S
Schneider, P. & Sluse, D. 2014, Astronomy and Astrophysics, 564, A103, aDS Bibcode: 2014A&A. 564A.103S
Seljak, U. 1994, The Astrophysical Journal, 436, 509, aDS Bibcode: 1994ApJ…436..509S
Shah, P., Lemos, P., & Lahav, O. 2021, Astronomy and Astrophysics Review, 29, 9, aDS Bibcode: 2021A&ARv..29. 9S
Shajib, A. J., Birrer, S., Treu, T., et al. 2019, Monthly Notices of the Royal As- tronomical Society, 483, 5649, arXiv:1807.09278 [astro-ph]
Shajib, A. J., Mozumdar, P., Chen, G. C. F., et al. 2023, A&A, 673, A9
Shajib, A. J., Mozumdar, P., Chen, G. C.-F., et al. 2023, Astronomy & Astro- physics, 673, A9, arXiv:2301.02656 [astro-ph]
Shajib, A. J., Treu, T., & Agnello, A. 2018, Monthly Notices of the Royal Astro- nomical Society, 473, 210, arXiv:1709.01517 [astro-ph]
Shajib, A. J., Wong, K. C., Birrer, S., et al. 2022, Astronomy & Astrophysics, 667, A123, arXiv:2202.11101 [astro-ph]
Shapiro, I. I. 1964, Phys. Rev. Lett., 13, 789
Sluse, D., Hutsemékers, D., Courbin, F., Meylan, G., & Wambsganss, J. 2012, Astronomy & Astrophysics, 544, A62, arXiv:1206.0731 [astro-ph]
Sluse, D., Rusu, C. E., Fassnacht, C. D., et al. 2019, Monthly Notices of the Royal Astronomical Society, 490, 613
Starck, J.-L., Moudden, Y., Abrial, P., & Nguyen, M. 2006, Astronomy & Astro- physics, 446, 1191, number: 3 Publisher: EDP Sciences
Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70
Suyu, S. H., Chang, T.-C., Courbin, F., & Okumura, T. 2018, Space Science Reviews, 214, 91, arXiv:1801.07262 [astro-ph]
Suyu, S. H. & Halkola, A. 2010, Astronomy & Astrophysics, 524, A94 Suyu, S. H. & Halkola, A. 2010, A&A, 524, A94
Suyu, S. H., Huber, S., Cañameras, R., et al. 2020, arXiv e-prints, 2002.08378 Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010a, ApJ, 711, 201
Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010b, The Astrophysical Jour- nal, 711, 201, publisher: The American Astronomical Society
Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, The Astrophysical Journal, 788, L35, arXiv:1306.4732 [astro-ph]
Sérsic, J. L. 1968, Atlas de Galaxias Australes, publication Title: Cordoba ADS Bibcode: 1968adga.book. S
Tan, C. Y., Shajib, A. J., Birrer, S., et al. 2024, Monthly Notices of the Royal Astronomical Society, 530, 1474, arXiv:2311.09307 [astro-ph]
Treu, T., Koopmans, L. V., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 640, 662
Treu, T. & Koopmans, L. V. E. 2002, The Astrophysical Journal, 575, 87, pub- lisher: IOP Publishing
Treu, T. & Marshall, P. J. 2016, The Astronomy and Astrophysics Review, 24, 11, arXiv:1605.05333 [astro-ph]
Treu, T., Suyu, S. H., & Marshall, P. J. 2022, The Astronomy and Astrophysics Review, 30, 8, arXiv:2210.15794 [astro-ph]
Unruh, S., Schneider, P., & Sluse, D. 2017, Astronomy and Astrophysics, 601, A77, aDS Bibcode: 2017A&A…601A..77U
Vanderriest, C., Schneider, J., Herpe, G., et al. 1989, Astronomy and Astro- physics, 215, 1, aDS Bibcode: 1989A&A…215. V
Verde, L., Treu, T., & Riess, A. G. 2019, Nature Astronomy, 3, 891, arXiv:1907.10625 [astro-ph, physics:gr-qc, physics:hep-ph, physics:hep-th]
Wells, P., Fassnacht, C. D., & Rusu, C. E. 2023, Astronomy & Astrophysics, 676, A95
Wells, P. R., Fassnacht, C. D., Birrer, S., & Williams, D. 2024a, Astronomy & Astrophysics, 689, A87, arXiv:2403.10666 [astro-ph]
Wells, P. R., Fassnacht, C. D., Birrer, S., & Williams, D. 2024b, TDCOSMO XVI: Population Analysis of Lines of Sight of 25 Strong Galaxy-Galaxy Lenses with Extreme Value Statistics, arXiv:2403.10666 [astro-ph]
Williams, K. A., Momcheva, I., Keeton, C. R., Zabludoff, A. I., & Lehár, J. 2006, The Astrophysical Journal, 646, 85
Wilson, M. L., Zabludoff, A. I., Ammons, S. M., et al. 2016, The Astrophysical Journal, 833, 194, publisher: IOP ADS Bibcode: 2016ApJ…833..194W
Wong, K. C., Keeton, C. R., Williams, K. A., Momcheva, I. G., & Zabludoff,
A. I. 2010, The Astrophysical Journal, 726, 84, publisher: The American As- tronomical Society
Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, Monthly Notices of the Royal Astronomical Society, 465, 4895, arXiv:1607.01403 [astro-ph]
Wong, K. C., Suyu, S. H., Chen, G. C.-F., et al. 2020, Monthly Notices of the Royal Astronomical Society, 498, 1420, arXiv:1907.04869 [astro-ph]
Wucknitz, O. 2002, Monthly Notices of the Royal Astronomical Society, 332, 951, arXiv:astro-ph/0202376
Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, Monthly Notices of the Royal Astronomical Society, 493, 4783, arXiv:1904.07237 [astro-ph]
Appendix A: Additional Systematic Tests
Appendix A.1: Polar Exponential Shapelet Bases
We tested representing our source complexity with polar expo- nential shapelet bases, which are theoretically expected to bet- ter encapsulate galaxy image information with fewer basis el- ements (Bergé et al. 2019). While our tests with polar expo- nential shapelets indeed yielded better BIC values, we found it came at the price of incorrect PSF reconstructions. This is likely caused by the increased information density near the origin of the shapelet model. This allowed for the extended source recon- struction to create point-source-like flux values at the center of the extended host’s light distribution. This is problematic, as we need to accurately disentangle the light contributions from the point source and the extended source, otherwise our astrometry will be incorrect, biasing our results.
Therefore, despite the improved performance of polar expo- nential shapelets, we continue using the standard linear shapelets (e.g., H0LiCOW IX, Shajib et al. 2022). We find their spatially- uniform information density effectively prevents them from overfitting the center of the extended quasar host galaxy light.
Appendix A.2: Two Extended Source Sérsics
All the source reconstructions from this system display a bulge- like component in the center of the host galaxy (see Fig. 5). This bulge may have preference for a distinct Sérsic index from the extended disk structure of the host, as we observe in nearby galaxies. Therefore, instead of reconstructing the source with one Sérsic light profile modified by linear shapelets, we at- tempted fitting an additional Sérsic profile (maintaining the same scale of linear shapelets). We explored both leaving the indices free and fixing them to typical locally observed values (one de Vaucouleurs profile at nsérsic = 4 and one exponential profile at nsérsic = 1). We found the BIC slightly, but consistently, preferred the simpler one-Sérsic extended source model.
Appendix A.3: Cutout Size
Although modeling larger cutouts aids in reconstructing the light profiles of the deflector and various perturbers, it comes at the cost of reduced weight on the extended source structure and drastic increases to computational time. To confirm the chosen mask sizes and positions (2.1′′ and 2.2′′ radii) compose a suf-
ficient representation of the lens system, we fit one model on a
significantly larger cutout and found no statistically significant discrepancies.
Appendix A.4: Fitting Likely Star-Forming Regions in Host
In addition, we attempted to add additional shapelets to account for the residuals found south west of image C. We found the additional computational cost outweighed the negligible affect on parameter constraints.
Appendix A.5: Scaling Deflector Einstein Radius with Perturbers
One possible modeling choice would be to scale the main deflec- tor’s Einstein radius with the rest of the galaxy-scale perturbers. We found it had negligible impact on the BIC and model results, but kept it free based on previous cosmography-grade lens mod- els.
Appendix A.6: Modeling the F150W Band
Test models were run on the F150W band to verify the precision and accuracy of the results when moving to larger wavelengths. No significant discrepancies were found in the inferred parame- ters, including the Fermat potential differences.
Appendix A.7: Supersampling the Extended Source Light
As supersampling the data is computationally intensive, regions that are supersampled are typically confined to the centers of the point sources. To test for potential improvements on parameters like γPEMD, we ran tests that also supersampled a large portion of the extended light of the host. We found the computational trade-off outweighed the marginal gain on parameter precision.
Appendix A.8: Freeing Lens Light Sérsic Profiles
We found leaving the Sérsic indices free for the main deflector’s light profiles led to an additional degeneracy with their resulting radii. This slowed model convergence, while having no notice- able impact on the resulting Fermat potential differences.
Appendix B: Model Weights
The model weights used for the final inference are shown in Ta- ble B.1. The models including flexion are weighted separately, as they are not used in the final inference.
Appendix C: Impact of PSF Noise Estimation
It is relatively challenging to estimate the additional uncertainty on the data caused by inaccuracies of the PSF model. However, as long as the residuals of the PSF are sufficiently amplified (without incorrectly amplifying the uncertainty of the extended source structure), this should have minimal impact on model re- sults as shown in previous analyses (e.g., H0LiCOW XII).
The residuals of our models–after removing this additional PSF uncertainty in the data–showed no discernible trends that would raise concern (Fig. C.1). We also note minimal over- lap between regions with amplified noise and the extended host structure, which proved a non-trivial task with JWST’s PSF (Fig. C.2). This is crucial to avoid down weighting the ring struc- ture of the extended source, which provides the tightest con- straints on the power-law slope.
Appendix D: Tests on Source Reconstruction Supersampling
Given that this represents the most detailed source reconstruc- tion (nmax = 18, 20, 22) among time-delay cosmography systems modeled with Lenstronomy, care must be taken to avoid super- sampling beyond a numerically feasible limit. If the number of shapelets (or similarly, the shapelet scale) pass the threshold of
information supported by the data, the model will be attempt- ing to resolve spatial scales smaller than the data’s pixel scale would constrain. This would lead to a nonphysical solution in the source light reconstruction, as the surface brightness of the source would no longer be conserved.
As we use linear shapelets for the parameterization of our source reconstruction, the minimum source scale reconstructed7
7 Technically, the “true" source scale depends on the spatially- dependent source plane magnification. However, this magnification al-
Table B.1. Final models sorted by total weight (kinematics weight × modeling/BIC weight).
|
Main Deflector |
Source (Shapelet Order) |
Perturbers |
PSF Type |
Mask (asec) |
∆BIC |
Total Weight |
|
No Free Flexion |
||||||
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7 |
STARRED |
2.2 |
38 |
1.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7 |
STARRED |
2.2 |
96 |
0.483 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7 |
STARRED |
2.1 |
190 |
0.094 |
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7 |
STARRED |
2.1 |
195 |
0.084 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7 |
STARRED |
2.2 |
303 |
0.005 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7 |
STARRED |
2.1 |
410 |
0.000 |
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7 |
PSFr |
2.2 |
1968 |
0.000 |
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7 |
PSFr |
2.1 |
2014 |
0.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7 |
PSFr |
2.2 |
2538 |
0.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7 |
PSFr |
2.1 |
2164 |
0.000 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7 |
PSFr |
2.1 |
2502 |
0.000 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7 |
PSFr |
2.2 |
2379 |
0.000 |
|
Free Flexion |
||||||
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7+Flexion |
STARRED |
2.1 |
0 |
0.000 |
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7+Flexion |
STARRED |
2.2 |
27 |
0.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7+Flexion |
STARRED |
2.1 |
63 |
0.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7+Flexion |
STARRED |
2.2 |
66 |
0.000 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7+Flexion |
STARRED |
2.2 |
226 |
0.000 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7+Flexion |
STARRED |
2.1 |
257 |
0.000 |
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7+Flexion |
PSFr |
2.1 |
1702 |
0.000 |
|
SPEMD + 2 SERSIC |
18 nmax + SERSIC |
X+G2+G3+G7+Flexion |
PSFr |
2.2 |
1778 |
0.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7+Flexion |
PSFr |
2.2 |
1852 |
0.000 |
|
SPEMD + 2 SERSIC |
20 nmax + SERSIC |
X+G2+G3+G7+Flexion |
PSFr |
2.1 |
2007 |
0.000 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7+Flexion |
PSFr |
2.2 |
2091 |
0.000 |
|
SPEMD + 2 SERSIC |
22 nmax + SERSIC |
X+G2+G3+G7+Flexion |
PSFr |
2.1 |
2218 |
0.000 |
Notes. The weighting is separated between the flexion and non-flexion model families. The ∆BIC values are calculated relative to the best performing (lowest-BIC) model. For more information, see Sect. 6.1.
is given by Appendix E: STARRED vs PSFr Parameter
Lmin =
β
√nmax + 1
. (D.1)
Comparison
Fig. E.1 shows cornerplots of the non-flexion models (as flex-
The shapelet scales reconstructed by our PSFr models were at ion was discarded for our results), comparing parameter posteri-
βPSFr
at β
= 0.089+0.010, while the STARRED models reconstructed
−0.005
= 0.074+0.006 (STARRED model scales were on av-
ors for models using PSFr vs STARRED to generate their initial PSF. In the comparison, the two samples are completely split and
STARRED
−0.004
assigned weights separately. STARRED had around eight mod-
erage 17% smaller than their PSFr counterparts). Ideally, this reconstruction scale should not be smaller than the spatial scale of our pixels, as this would imply that we are reconstructing at scales not constrained by the data. In other words, we require
Lmin ≥ δpix, (D.2)
~
with our pixels in the drizzled F115W band having a spatial res- olution of δpix 0.0307 arcsec /pix.
However, when supersampling the source, this form changes
to reflect the increased source resolution assumed to be con- strained by the supersampled data. For a supersampling factor of 3, the spatial scales in the image plane that constrain the source ray tracing are 3 times smaller than the true pixel scale of the data, yielding the new constraint
Lmin ≥ δpix/ fss, (D.3)
with our models taking a supersampling factor of fss = 3. This yields a minimum reconstruction scale allowed by the data of Lmin = 0.0103, while the smallest value sampled by any of our
models was Lmin, model = 0.0126, still above the threshold for
numerical stability.
ways increases the minimum allowed source scale. In this case, the anal- ysis that follows applies even more strongly.
els with significant weights, while PSFr had around three. Even though both PSFs were updated during the fitting, the disagree- ment in initial PSF leads to a more than 1-σ disagreement in parameters like the power-law slope, demonstrating the signifi- cance of accurate PSF modeling.
No PSF Error
Noise Map Normalized Residuals
B
A2
A1
C
1"
B
A2
A1
C
1"
|
|
B
A2
A1
C
1"
B
A2
A1
C
1"
|
|
|
|
With PSF Error
|
2 0
log10 (
|
data)
|
5 0 5
(fmodel
|
fdata)/
Fig. C.1. Without an additional correction for the uncertainty in the PSF, the fitting of the point sources drives the fit while providing lit- tle information (example is with STARRED PSF model). Top left: The noise map (or uncertainty in the data for each pixel) derived from the imaging data alone. Top right: Model normalized residuals without ad- ditionally accounting for PSF uncertainty. Large residuals near the im- age positions dominate the fit, while providing little information about the lens system. Bottom left: The noise map incorporating the addi- tional uncertainty from PSF modeling, estimated by fitting stars in the field (Sect. 4.1). Bottom right: Model normalized residuals after in- cluding PSF uncertainty are more evenly distributed, allowing the mod- els to focus on accurately describing the ring structure, which offers the
strongest constraints on the mass slope, γPEMD—arguably the most crit- ical parameter in power-law models.
Reduced Residuals
PSF Error
Noise Map
Data
5
0
E
1''
N
|
|
|
0
E
1''
N
|
|
|
0
E
1''
N
|
|
|
E
1''
N
|
|
|
|
|
(fmodel fdata)/
log10 (
data)
log10 (
data)
2 2
data) (radial profile)
2
0
log10 (
|
2
|
4
Image B
5
PSF
Data
equal
Image A2
2
0
|
2
Image A1
4
2
0
|
2
Image C
2
0
|
2
|
4
0.0 0.22 0.43
asec ('') from center
0.0 0.22 0.43
asec ('') from center
0.0 0.22 0.43
asec ('') from center
0.0 0.22 0.43
asec ('') from center
Fig. C.2. Radial breakdown of image uncertainty contributions between the noise directly in the data and the additional PSF uncertainty term, showing the PSF contribution does not down-weight the ring structure of the extended source. Top left: normalized residual map, same as Fig. C.1. Top middle left: error contributions from the PSF uncertainty alone, and top middle right shows only the uncertainty due to the data alone. These are combined to form the final map, shown in Fig. C.1. Top right: the imaging data used for modeling, with circles indicating the largest radial extent of the images where the PSF uncertainty equals the noise uncertainty. Shown on the bottom are the radial uncertainty contributions from
the imaging data vs the estimated PSF error. We find the maximum extent at which the PSF error dominates the uncertainty is at a radius less than 11 pixels (≈ 0.338′′), shown in the top right figure.
|
|
|
E, PEMD
0.95
|
0.90
|
|
|
PEMD
2.0
|
1.9
|
|
|
|
Sample (n=1000) PSF: STARRED PSF: PSFr
(HST: Rusu et al. 2020)
|
qPEMD
0.80
0.75
|
|
|
|
|
PEMD
35
|
30
|
|
|
|
|
|
0.12
ext
|
0.10
|
|
|
|
|
|
0.08
|
|
ext
100
80
|
|
|
|
|
|
|
|
|
|
|
E, X
0.10
|
0.05
0.00
|
|
|
|
|
|
|
|
|
E, G2
1.0
|
0.8
|
|
|
|
|
|
|
|
|
|
|
V
250
|
200
0.45 0.50 0.90 0.95
|
1.9 2.0
0.75 0.80
30 35
|
0.075 0.100 0.125
80 100
|
0.0 0.1
|
0.8 1.0
200 250
|
|
BC E, PEMD
PEMD
qPEMD
|
PEMD
ext
|
ext
E, X
E, G2
|
V
Fig. E.1. Same as Fig. 7, except comparing weighted PSFr models to weighted STARRED models. Both groups receive BIC weights independently for this comparison.