© A. J. Shajib et al. 2022

1. Introduction

The Hubble constant, H₀, is a central cosmological parameter as it sets the expansion rate of the Universe. Consequently, precise knowledge of its value is crucial for our understanding of the Cosmos, and it also has important implications in extragalactic astrophysics. However, different methods have measured the Hubble constant with discrepant values, producing the so-called Hubble tension (e.g. Freedman 2021). Mapping of the temperature fluctuations of the cosmic microwave background allows one to measure the Hubble parameter H(z ≈ 1100) at the last scattering surface, and then the Hubble constant, H₀, at the current epoch is extrapolated using Λ cold dark matter (ΛCDM) cosmology. This early-Universe probe resulted in constraints of H₀ = 67.4 ± 0.5 km s⁻¹ Mpc⁻¹ (Planck Collaboration VI 2020) and H₀ = 67.6 ± 1.1 km s⁻¹ Mpc⁻¹ (Aiola et al. 2020). In the local Universe, H₀ is typically measured by building a cosmic distance ladder up to type Ia supernovae (SNe) in the Hubble flow by calibrating their absolute magnitudes with intermediate distance probes. The Supernova H₀ for the Equation of State of dark energy (SH0ES) team used Cepheids and parallax distances to calibrate the cosmic distance ladder, measuring H₀ = 73.04 ± 1.04 km s⁻¹ Mpc⁻¹ (Riess et al. 2022), which is in 5σ tension with the Planck measurement. The Carnegie–Chicago Hubble Project used the tip of the red giant branch to calibrate the distance ladder, measuring H₀ = 69.6 ± 1.9 km s⁻¹ Mpc⁻¹ (Freedman et al. 2019, 2020), which, interestingly, is statistically consistent with both the SH0ES and Planck measurements. Several other local probes strengthened the Hubble tension, for example the Megamaser Cosmology Project measured H₀ = 73.9 ± 3.0 km s⁻¹ Mpc⁻¹ (Pesce et al. 2020), the Tully–Fisher method calibrated with Cepheids measured H₀ = 75.1 ± 0.2 ± 3.0 km s⁻¹ Mpc⁻¹ (Kourkchi et al. 2020), and the surface brightness fluctuation method measured H₀ = 73.7 ± 0.7 ± 2.4 km s⁻¹ Mpc⁻¹ (Blakeslee et al. 2021). If systematics in these measurements can be ruled out as the source of this Hubble tension, new physics beyond the standard ΛCDM cosmology will be required to resolve the tension (e.g. Poulin et al. 2019; Knox & Millea 2020; Efstathiou 2021). Therefore, thoroughly investigating the potential systematics that are as yet unknown in each of the probes is paramount.

Strong-lensing time delays provide an independent probe of the Hubble constant (Refsdal 1964). The delays between the arrival times of photons corresponding to different images of the background source depend on the cosmological distances involved in the strong-lensing system, and thus these delays allow us to measure a combination of these distances, called the ‘time-delay distance’ (Suyu et al. 2010). The time-delay distance is inversely proportional to H₀ and weakly dependent on other cosmological parameters. Although early implementations of this method in the 1990s and the early 2000s suffered from limitations in data quality and analysis techniques, both of these aspects have improved by a large margin over the past decade (for a review with a historical perspective, see Treu & Marshall 2016). Inferring H₀ from the time delays requires: (i) measuring the time delays, (ii) measuring the redshifts of the deflector and the background quasar, (iii) modelling the mass distribution in the central deflector to compute the Fermat potential differences between the image positions, and (iv) estimating the extra lensing contribution from the line-of-sight (LOS) mass distribution between the background source and the observer. Thanks to breakthroughs in all of these factors, the H₀ Lenses In the COSMOGRAIL’s Wellspring (H0LiCOW) and the Strong-lensing High Angular Resolution Programme (SHARP) collaborations measured km s⁻¹ Mpc⁻¹ from a sample of six strongly lensed quasar systems (Suyu et al. 2017; Bonvin et al. 2017; Birrer et al. 2019; Chen et al. 2019; Rusu et al. 2020; Wong et al. 2020). The STRong-lensing Insights into the Dark Energy Survey (STRIDES) collaboration analysed a seventh lens system to measure km s⁻¹ Mpc⁻¹ (Shajib et al. 2020). It is noteworthy that six out of these seven analyses were performed blindly, with only the first one being a non-blind analysis. The H0LiCOW, STRIDES, Cosmological Monitoring of Gravitational Lenses (COSMOGRAIL), and SHARP collaborations have united under the umbrella of the Time-Delay COSMOgraphy (TDCOSMO) collaboration.

The TDCOSMO collaboration has already performed a number of tests to search for previously unknown systematics. Millon et al. (2020, TDCOSMO-I) checked for systematics arising from the current treatments of the stellar kinematics, LOS mass distribution, and the choice of lens model families, finding no evidence for unaccounted errors. Gilman et al. (2020, TDCOSMO-III) find that dark sub-halos – which are ignored in lens modelling through the assumption of smooth mass profiles – also do not systematically bias the H₀ inference, adding negligible random uncertainty. Birrer et al. (2020, TDCOSMO-IV) relaxed the assumption of the power-law mass distribution in the deflector galaxies to allow maximal degeneracy in the mass distribution under the mass-sheet transformation (MST; Falco et al. 1985). By constraining the mass distribution from the stellar kinematics only, these authors inferred km s⁻¹ Mpc⁻¹, that is, relaxing the power-law assumption leads to an increase in H₀ uncertainty from 2.2% to 7.9% for the sample of the seven analysed systems. To regain the lost precision, TDCOSMO-IV combined an external sample of galaxy–galaxy strong lenses from the Sloan Lens ACS¹ (SLACS) survey to add more information on the galaxy mass distribution, under the assumption that the SLACS lenses and the TDCOSMO lenses belong to the same galaxy population. Adding a sample of 33 SLACS lenses improved the precision to 5.4%. Although the point estimate of H₀ shifted to km s⁻¹ Mpc⁻¹ with the addition of the SLACS lenses, this value is still consistent with all the previous TDCOSMO measurements within 1σ. Birrer & Treu (2021, TDCOSMO-V) forecasted that a future sample of 40 time-delay lenses with spatially resolved stellar kinematics and an external lens sample of 200 non-time-delay lenses will be able to infer H₀ with 1.2−1.3% precision, which is necessary to independently settle the Hubble tension at the ∼5σ confidence level. Van de Vyvere et al. (2022a, TDCOSMO-VII) find that the systematic bias in the measured H₀ arising from the boxy-ness or discy-ness of the deflector galaxy is < 1% and thus insignificant. Blind data challenges are also important tests for the presence of systematics. The Time-Delay Challenge validated the robustness of the methods currently used to measure time delays from quasar light curves (Dobler et al. 2015; Liao et al. 2015). The Time-Delay Lens Modelling Challenge similarly validated the modelling techniques currently used to recover the ground truth when the shapes of the underlying galaxy mass profiles are known (Ding et al. 2021).

In this paper we present the results of an experiment to search for potential systematics in the lens modelling – within specific assumed mass profile families – that may arise from different modelling software programs used by different investigators. In this experiment, two teams using different software programs independently modelled the strongly lensed quasar system WGD 2038−4008 to the level required for cosmographic application (i.e. to the noise level; Agnello et al. 2018). The two modelling software programs being compared are GLEE² and LENSTRONOMY³. The core members of the GLEE team are K. C. Wong and S. H. Suyu; the core members of the LENSTRONOMY team are A. J. Shajib, S. Birrer, and T. Treu. Both of the software programs have previously been used for lens modelling in cosmographic analyses by the TDCOSMO collaboration – five systems with GLEE and two systems with LENSTRONOMY. Although Birrer et al. (2016) performed a cosmographic analysis outside the TDCOSMO umbrella using LENSTRONOMY for the system RX J1131−1231, which was previously analysed by the H0LiCOW collaboration using GLEE, a systematic blind comparison between the two software programs on the same lens system has not been done previously. Both software programs perform parametric modelling of the deflector mass distribution, but they differ in the method used for source reconstruction. Whereas GLEE uses a pixel-based source reconstruction with regularization conditions (Suyu et al. 2006), LENSTRONOMY uses a basis set of parameterized profiles for source reconstruction (Birrer et al. 2015, 2021; Birrer & Amara 2018).

In addition to the software architectures, differences in the lens models may arise from modelling choices made by an investigator in such modelling processes. Our experiment also encompasses this human aspect of the modelling process by having the two teams work independently and blindly. However, to facilitate a fair comparison between the model predictions, we established a baseline model setup with minimal specifications that was agreed upon by the two teams before performing their own analyses. After each team separately completed their internal systematic checks and went through an internal review by the TDCOSMO collaboration, the lens models were frozen and the model predictions were un-blinded to make comparisons between the two teams. As the time delay for this system has not yet been measured with sufficient precision for an H₀ measurement, we leave the H₀ inference from our models to be done in the future. However, we predict the time delays for this system as a function of H₀ after marginalizing over the inferences from the two modelling software programs. As a result, our ‘preemptive’ lens models enforce an additional layer of blindness for the future H₀ measurement from this system.

The baseline models for comparison have two different lens model setups: (i) a power-law mass model and (ii) a two-component mass model that individually accounts for the dark and baryonic components. It is well known that conventional parametric models such as the power-law model impose assumptions that break the mass-sheet degeneracy (MSD; e.g. Birrer et al. 2020; Kochanek 2020). However, a lens model is still useful for extracting the relevant lensing information (i.e. the Fermat potential difference) from the data, which can then be processed to allow the additional freedom along the MSD following TDCOSMO-IV. Although techniques to extract lensing information without relying on parametric models have recently been proposed (e.g. Birrer 2021), they have not yet been applied to real systems for rigorous lens modelling similar to the TDCOSMO analyses. Furthermore, no evidence has so far demonstrated that the simply parametrized models are not an adequate description, and the necessity or physical reality of a mass component that acts as a physical mass sheet has not been demonstrated. For all these reasons, until new evidence is gathered to inform new choices, simply parametrized lens models are going to be the baseline in TDCOSMO analyses. Therefore, it is important to compare the modelling methods based on these software programs to check for systematic differences as performed in this paper.

In this paper we only predict the time delays for WGD 2038−4008 based on our lens models, as the actual time delays for this system are yet to be measured and thus the H₀ cannot be inferred. Measuring H₀ based on the lens models presented in this paper is left for a future paper.

This paper is organized as follows. In Sect. 2 we provide a brief review of the strong lensing formalism to establish the notations and describe the Bayesian inference framework of our model predictions. The observables in our analysis are described in Sect. 3. We present the baseline models that are common to both teams in Sect. 4. The modelling procedures and results are presented by the GLEE and LENSTRONOMY teams in Sects. 5 and 6, respectively. We compare and discuss the results from the two teams in Sect. 7 and conclude the paper in Sect. 8. Sections 1–6 were written prior to the un-blinding. After un-blinding on October 22, 2021, Sects. 7 and 8 were written and no major edits were done to Sects. 1–6, except for minor fixes for typos and grammatical errors.

2. Framework of the lens modelling

In this section we describe the theoretical framework for our analysis. We give a brief overview of the strong lensing formalism in Sect. 2.1, discuss the MSD in Sect. 2.2, explain our modelling of the stellar kinematics in Sect. 2.3, and present the Bayesian inference framework for our analysis in Sect. 2.4.

2.1. Strong lensing formalism

The goal of this section is to provide the necessary definitions in strong lensing and establish the notation. This formalism was developed in multiple previous studies (see e.g. Schneider et al. 1992; Blandford & Narayan 1992) and has been implemented in numerous previous TDCOSMO analyses (e.g. Suyu et al. 2010; Birrer et al. 2019; Shajib et al. 2020).

The delay Δt_XY between arrival times of photons corresponding to images labelled as X and Y is given by

(1)

Here, D_d is the angular diameter distance to the deflector, D_s is that to the source, and D_ds is that between the deflector and the source, z_d is the deflector redshift, c is the speed of light, θ is the image position, β is the un-lensed source position, and ψ is the deflection potential that is related to the deflection angle as ∇ψ ≡ α and the convergence as ∇²ψ = 2κ. The convergence is the surface mass density scaled by the critical density as κ ≡ Σ/Σ_crit with

(2)

The Fermat potential ϕ is defined by combining the geometric delay term with the deflection potential as

(3)

The so-called time-delay distance is defined as

(4)

Each distance term contains a factor of , which cancel out such that . Equation (1) can be written in short form as

(5)

2.2. Mass-sheet degeneracy

The imaging observables of the lensing phenomenon – the image positions and the flux ratios – remain invariant under the transformation

(6)

which is referred to as the MST (Falco et al. 1985). The invariance of the observables under this transformation gives rise to the MSD. We note that the magnifications are not invariant under the MST (although magnification ratios are), and thus strongly lensed standard candles can break the MSD (Bertin & Lombardi 2006).

We can separate all of the mass contributing to lensing of the background source into two components as

(7)

where κ_cen is the convergence from the central deflector and κ_ext is the convergence from all the LOS mass distribution – except the central deflector – projected onto the plane of the central deflector (i.e. the image plane). In some cases, the central deflector may have nearby companions or satellites, or nearby LOS perturbing galaxies that are explicitly accounted for in the lens model, for example RX J1131−1231, HE 0435−1223, and ES J0408−5354 (Suyu et al. 2013; Wong et al. 2017; Shajib et al. 2020). We consider these additional mass components to be included in κ_cen. As the mass distribution of the central deflector goes to zero at very large radius, we have

(8)

Therefore, κ_ext can be interpreted as lensing mass in the 3D space far from or ‘external’ to the central deflector. Let be the model convergence that can reproduce the imaging observables. However, due to the MSD, is not a unique solution and we cannot ascertain that . If we impose the condition , then is a mass-sheet transform of κ_true with the rescaling factor λ = 1/(1 − κ_ext) as

(9)

If the external convergence κ_ext can be independently estimated by studying the lens environment, then the true lensing convergence κ_true can be recovered from through the corresponding inverse MST. However, the lens model κ_model that we actually constrain can be an internal MST of as

(10)

Interestingly, both κ_model and can go to zero at θ → ∞ by construction. In such a case, λ_int is not a constant and it satisfies lim_{θ → ∞} = 1 (Schneider & Sluse 2014). We can combine Eqs. (8)–(10) to write the relation between the true mass distribution κ_true and the modelled mass distribution κ_model as

(11)

To constrain λ_int, we require observables that rescale with the MST, for example the stellar kinematics. Although such observables rescale with λ_int(1 − κ_ext), the external convergence κ_ext is independently estimated from the LOS properties leaving only λ_int to be constrained from those observables. The LOS velocity dispersion rescales with the MST as

(12)

This rescaling is only valid for a pure MST, such as the external MST, and is approximately valid for an internal MST with single aperture kinematics. However, this is not valid for internal MST with spatially resolved kinematics (Chen et al. 2021; Yıldırım et al. 2021). The time delay rescales with the MST as

(13)

As a result, we need to correct the time delays Δt_model predicted by the model κ_model as

(14)

In the next section, we describe our framework for the kinematics analysis.

2.3. Kinematics analysis

The stellar velocity dispersion probes the 3D mass distribution of the deflector galaxy that is deprojected from κ_cen. We adopt the spherical Jeans equation that connects the velocity dispersion with the gravitational potential Φ(r) as

(15)

Here, l(r) is the 3D luminosity density, σ_r(r) is the radial velocity dispersion, and β_ani(r) is the anisotropy parameter that relates σ_r to the tangential velocity dispersion σ_t as

(16)

The observable quantity is the luminosity-weighted LOS velocity dispersion, which we can obtain by solving the Jeans equation as

(17)

where G is the gravitational constant, I(R) is the surface brightness, and M(r) is the 3D enclosed mass within radius r (Eqs. (A.15) and (A.16) of Mamon & Łokas 2005). The function 𝒦_β(r/R) depends on the parameterization of β_ani(r). We adopt the Osipkov–Merritt parameterization given by

(18)

where r_ani is a scaling radius (Osipkov 1979; Merritt 1985a,b). For this parameterization, the form of 𝒦_β(r/R) is given by

(19)

with u_ani ≡ r_ani/R (Mamon & Łokas 2005). The observed aperture-averaged velocity dispersion is

(20)

where ∫_ap denotes integration over the aperture and 𝒮 denotes convolution with the seeing. Thus, the lens-model-predicted LOS velocity dispersion can be written in the form

(21)

where ξ_lens is the set of mass model parameters and ξ_light is the set of light distribution parameters. The internal and external MST parameters modify the lens-model-predicted velocity dispersion as

(22)

The dependence of σ_ap on the cosmology is fully captured in the D_s/D_ds term. The function J is independent of cosmology as all of its arguments are expressed in angular units, but it should be noted that J is directly connected to the model convergence κ_model through the parameters ξ_lens (Birrer et al. 2016).

2.4. Bayesian inference

We denote the set of all the observables as O ≡ {O_img, O_kin}, where O_img is the imaging data of the lens system and O_kin is the measured stellar velocity dispersion. Although data from spectroscopic and photometric surveys of the lens environment are necessary to estimate the external convergence, we fold in the estimated external convergence as the prior p(κ_ext) in our inference. To predict the time delay for a given cosmology, we want to infer the Fermat potential difference Δϕ between the corresponding image pairs. The Fermat potential difference Δϕ(ξ, κ_ext, λ_int) is a function of the set of model parameter ξ ≡ {ξ_lens, ξ_light, r_ani} in a model family M, external convergence κ_ext, and internal MST parameter λ_int. Thus, to obtain p(Δϕ ∣ O), we first aim to infer p(ξ, κ_ext, λ_int ∣ O). Applying Bayes’ theorem, we can write

(23)

Here, S is the set of lens model hyper-parameters that is only relevant for O_img, and D_s/ds is a short notation for the distance ratio D_s/ds ≡ D_s/D_ds. We explicitly separate the hyper-parameters S – that need to be fixed during optimizing a lens model, for example the set of pixels for computing the image likelihood, resolution of the source reconstruction – from the choice of lens model family M. The prior p(κ_ext ∣ M) depends on the model family M, since the model-constrained shear is used to estimate κ_ext corresponding to M. Since O_img and O_kin are independent data, the likelihood term p(O ∣ ξ, M, S, D_s/ds, κ_ext) can be decomposed as

(24)

Then, we can first perform the following sub-integral within the right-hand side of Eq. (23):

(25)

Here, p(O_img ∣ M, S) is the model evidence. We perform this integral in the form of the right-hand side of Eq. (25) for numerical convenience, as it allows us to first obtain the posterior p(ξ ∣ O_img, M, S) using Monte Carlo sampling, and then combine the posteriors weighted by the model evidence to perform the integration in Eq. (25). We use the Bayesian information criterion (BIC) as a proxy for the model evidence in our analysis (Schwarz 1978). The BIC is defined as

(26)

where k is the number of free parameters, N_data is the number of data points, and is the maximum of the likelihood function ℒ. Both the BIC and directly computed model evidence were used in previous analyses for Bayesian model averaging (BMA; e.g. Madigan & Raftery 1994; Hoeting et al. 1999) in the context of lens modelling for cosmographic analysis (BIC: Birrer et al. 2019; Chen et al. 2019; Rusu et al. 2020; model evidence: Shajib et al. 2020).

Specific implementations of the Bayesian inference framework presented in this section through sampling by each team are described in Sects. 5 and 6.

3. Imaging data and ancillary measurements

The system WGD 2038−4008 was discovered from a combined search in the Wide-field Infrared Survey Explorer and Gaia data over the Dark Energy Survey (DES) footprint (Agnello et al. 2018). The deflector redshift is z_d = 0.230 ± 0.002 and the source redshift is z_s = 0.777 ± 0.001 (Agnello et al. 2018). In this section we describe the imaging data and spectroscopic measurements used in our analysis.

3.1. HST imaging

We obtained Hubble Space Telescope (HST) imaging of the system (GO-15320, PI: Treu; Shajib et al. 2019) using the Wide-Field Camera 3 (WFC3). The imaging was taken in three filters: F160W in the infrared (IR) channel, and F814W and F475X in the ultraviolet-visual (UVIS) channel. Four exposures were taken in each filter to cover the large dynamic range in surface brightness of the brighter quasar images and the fainter extended host galaxy. For the IR band, we adopted a four-point dither pattern and STEP100 readout sequence for the MULTIACCUM mode. The total exposure times are 2196.9 s, 1428.0 s, and 1158.0 s, respectively, in the three filters. We show a false-colour red-green-blue (RGB) image of the system created from the HST imaging in Fig. 1.

Fig. 1. False-colour image of the lens systems WGD 2038−4008. This RGB image is created from the F160W (red), F814W (green), and F475X (blue) filters of the HST WFC3. We adjusted the relative amplitudes between the three filters to achieve a higher contrast for better visualization. The four lensed quasar images are marked as A, B, C, and D.

The point spread function (PSF) corresponding to each filter is estimated from stacking 4−6 stars that are within each corresponding HST image. These PSFs are only used as an initial estimate by both teams and they are refined to more accurately match the PSF at the quasar image positions by iterative reconstruction during the lens model optimization (see Sects. 5 and 6 for more details on the iterative reconstruction).

3.2. Stellar velocity dispersion

Buckley-Geer et al. (2020) measure the stellar velocity dispersion of the deflector from spectroscopic observation using the Gemini Multi-Object Spectrograph (GMOS-S) on the Gemini South Telescope. The measured velocity dispersion is σ_los = 296 ± 19 km s⁻¹ from a rectangular aperture, which is in agreement with a more recent measurement from the X-shooter instrument on the Very Large Telescope (VLT; Melo et al. 2021). We used the measurement from Buckley-Geer et al. (2020) instead of the more precise measurement from Melo et al. (2021) because the latter was published after the un-blinding, when the lens models were frozen and utilized the previous measurement. The seeing full width at half maximum (FWHM) is , and the exponent parameter of the Moffat PSF is β = 1.74.

3.3. LOS environment

The LOS environment of the system WGD 2038−4008 was studied by Buckley-Geer et al. (2020). These authors estimated the external convergence based on the weighted galaxy number counts approach (Greene et al. 2013; Rusu et al. 2017, 2020; Birrer et al. 2019). The weighted number counts were obtained in two separate apertures with radii 45″ and 120″ centred on the lens system from the DES multi-band imaging. The magnitude limit of counted galaxies is I = 22.5 mag. The counts are weighted based on simple physical quantities, such as the inverse of the distance to the lens. The spectroscopic redshifts were obtained from Gemini South GMOS-S and the photometric redshifts are based on DES multi-band photometry. Analogous number counts are also obtained within a large number of different apertures with the same sizes along random LOSs in the DES footprint. By comparing the weighted number counts for the LOS around WGD 2038−4008 with those for random LOSs, the over- or under-density is estimated in terms of a weighted number count ratio. The external convergence is then estimated by comparing the weighted number count ratio with that from statistically similar LOSs from the Millennium simulation with computed external convergence (Springel et al. 2005; Hilbert et al. 2009). If no external shear is considered, then the system WGD 2038−4008 was found to be along a LOS with approximately no overdensity within ∼1% uncertainty. We provide the κ_ext re-weighted based on the best-fit external shear magnitudes from our lens models in Sects. 5 and 6. Buckley-Geer et al. (2020) also find that no nearby LOS perturbers are significant enough that they need to be included explicitly in the lens mass modelling.

4. Setup of baseline models

In this section we describe the baseline models that were initially agreed upon by the two teams before performing separate and independent lens modelling. In our baseline models, we use two families of mass models for the central deflector: (i) a power-law profile, and (ii) a composite profile with an elliptical NFW potential for the dark component and a superposition of three Chameleon profiles (hereafter, triple Chameleon profile) in convergence for the luminous component. We also add external shear to both types of mass model. For the light profile of the central deflector, we adopt a triple Sérsic profile in all three bands in the models with the power-law mass profile. In the models with the composite mass profile, however, we adopt a triple Chameleon light profile in the F160W band linked with the triple Chameleon mass profile and a triple Sérsic profile in the UVIS bands.

We adopted three Chameleon profiles to sufficiently account for the complexity in the light profile of the deflector. Moreover, we adopted the triple Chameleon light profile only for the F160W profile, since this is the only band that is connected to the luminous component of the convergence profile.

Although both teams adopted these baseline models, individual teams were allowed to make their own choices – which may not necessarily be identical – pertaining to other model specifications, for example parameter priors and fixing parameter values.

In the next subsections we provide the definitions of the mass and light profiles in the baseline models.

4.1. Mass profiles

The two baseline lens model families we adopt are the power-law mass profile and the composite mass profile.

4.1.1. Power-law mass profile

We adopted the power-law elliptical mass distribution (PEMD; Barkana 1998) defined as

(27)

where γ is the logarithmic slope, θ_E is the Einstein radius, and q_m is the axis ratio. The coordinates (θ₁,  θ₂) are in the coordinate frame that is aligned with the major and minor axes. The position angle of this frame is φ_m with respect to the RA–Dec frame.

4.1.2. Composite mass profile

The composite mass profile consists of two individual mass profiles for the baryonic and the dark components of the mass distribution.

For the dark matter distribution, we adopt a Navarro–Frenk–White (NFW) profile with ellipticity defined in the potential. The 3D NFW profile in the spherical case is given by

(28)

where ρ_s is the density normalization, and r_s is the scale radius (Navarro et al. 1997). We refer to Golse & Kneib (2002) for the expressions of the lens potential and deflection angles associated with the elliptical NFW profile.

For the baryonic mass distribution, we adopt the Chameleon convergence profile. The Chameleon profile matches with the Sérsic profile within a few per cent at 0.5−3θ_eff, where θ_eff is the half-light or effective radius of the Sérsic profile (Dutton et al. 2011). The Chameleon profile is defined as the difference between two non-singular isothermal ellipsoids:

(29)

where a₀ is the normalization and w_c and w_t are the core sizes for the individual non-singular isothermal components in the Chameleon profile (Dutton et al. 2011; Suyu et al. 2014). This profile is numerically convenient for computing lensing quantities using closed-form expressions unlike the Sérsic profile.

4.2. Light profiles of the deflector

4.2.1. Sérsic profile

The Sérsic profile is defined as

(30)

where I_eff is the amplitude, θ_eff is the effective radius along the intermediate axis, and n_s is the Sérsic index (Sérsic 1968). The factor b_n is a normalizing factor so that θ_eff is the half-light radius.

4.2.2. Chameleon light profile

In the composite baseline model, we use the same Chameleon profile from Eq. (29) for the light profile of the deflector, but replacing the convergence amplitude κ₀ with the flux amplitude I₀.

5. GLEE modelling

In this section we describe the GLEE modelling procedure. GLEE is a software package developed by S. H. Suyu and A. Halkola (Suyu & Halkola 2010; Suyu 2012). The lensing mass distribution is described by a parameterized profile. The lensed quasar images are modelled as point sources on the image plane convolved with the PSF. The extended host galaxy of the lensed quasar is modelled on a 50 × 50 pixel grid with curvature regularization (Suyu et al. 2006), spanning the range of source coordinates corresponding to the pixels within a region containing the lensed arcs (hereafter, referred to as the ‘arcmask’). The quasar image amplitudes are independent of the extended host galaxy light distribution to allow for deviations due to microlensing, time delays, and substructure. The lens galaxy light distribution is represented as the sum of three Sérsic (or three Chameleon) profiles with a common centroid.

The lens model is constrained by the positions of the lensed quasar images and the surface brightness of the pixels of the lensed Einstein ring of the quasar host galaxy in the three HST bands that are fit simultaneously. The quasar positions are fixed to the positions of the point sources on the image plane (after they have stabilized) and are given a fixed Gaussian uncertainty of width to account for offsets due to substructure in the lens or LOS. This uncertainty is small enough to satisfy astrometric requirements for cosmography (Birrer & Treu 2019). The quasar flux ratios are not used as constraints, as they can be affected by microlensing, which has been detected in this system (Melo et al. 2021). We use the initial PSF estimate in each band that was created from ∼4−6 bright stars within the HST image (Sect. 3.1). We first model the lens separately in each band to iteratively update the respective PSFs using the lensed active galactic nucleus (AGN) images (Chen et al. 2016; Wong et al. 2017; Rusu et al. 2020). We then keep the ‘corrected’ PSFs fixed and use them in our final models that simultaneously use the surface brightness distribution in all three bands as constraints. We use the positions of the quasar images to align the cutouts in the three HST bands. We do not enforce any similarity of pixel values at the same spatial position across different bands (i.e. the model flux at any position in one band is independent of the model flux in other bands). In our Markov chain Monte Carlo (MCMC) sampling, we vary the light parameters of the lens galaxy and quasar images, the mass parameters of the lens galaxy, and the external shear. The source position is also sampled in the modelling. The quasar image positions are linked across all bands, but the other light parameters are allowed to vary independently.

We create cutouts of the HST images with dimensions of , which corresponds to a 140 × 140 pixel cutout for the UVIS/F475X and UVIS/F814W bands and a 70 × 70 pixel cutout for the IR/F160W band. This conservative cutout size is chosen to include the entire region containing the lensed host galaxy arc light. We define the arcmask around the deflector galaxy in each of the three bands, which encloses the region where we reconstruct the lensed arc from the extended quasar host galaxy. The arcmask is used to calculate the likelihood involving the reconstructed lensed arc light, but the whole cutout is used for calculating the likelihood associated with the lens light. The construction of the weight images and bad pixel masking for each cutout are analogous to the procedure in Wong et al. (2017) and Rusu et al. (2020). In order to avoid biasing the modelling due to large residuals from a PSF mismatch near the AGN image positions, we rescale the weights in those regions by a power-law model such that a pixel originally given an estimated 1σ noise value of σ_img, i is rescaled to a noise value of . The constants A and b are chosen for each band such that the normalized residuals (the residual flux of each pixel normalized by its 1σ uncertainty) in the AGN image regions are approximately consistent with the normalized residuals in the rest of the arc region. We do not rescale the weights outside of the AGN image regions. The arcmask region and the regions around the AGN with rescaled weights are shown in the first column of Fig. 2.

Fig. 2. Fiducial power-law model results for IR/F160W (top row), UVIS/F814W (middle row), and UVIS/F475X (bottom row) from GLEE. The maximum-likelihood model in the MCMC chain is shown. Shown are the observed image (first column), the reconstructed image predicted by the model (second column), the normalized residual within and interior to the arcmask region (defined as the difference between the data and model, normalized by the estimated uncertainty of each pixel; third column), and the reconstructed source (right column). First column: the dotted cyan lines indicate the arcmask (donut-shaped) region used for fitting the extended source, the dotted orange lines indicate the AGN mask region where the power-law weighting is applied, and the region outside the dotted cyan arcmask is used to further constrain the foreground lens light and (partly) the AGN light (but not the AGN host galaxy light since its corresponding lensed arcs are below the noise level in this outer region). The colour bars show the scale in the respective panels. The results shown here are for the fiducial power-law model, but the results for the other systematics tests (Sect. 5.3) are qualitatively similar.

5.1. Power-law model

Our fiducial power-law mass model uses the triple Sérsic parameterization for the lens galaxy light and has the additional free parameters: (i) position (θ₁, θ₂) of the mass centroid (allowed to vary independently from the centroid of the light distribution), (ii) Einstein radius θ_E, (iii) minor-to-major axial ratio, q_m, and associated position angle φ_m (measured east of north), (iv) 3D slope of the power-law mass distribution γ, and (v) external shear γ_ext and associated position angle φ_ext (measured east of north).

We assume uniform priors on the model parameters over a wide physical range. Figure 2 shows the data and the lens model results in all three bands for our fiducial power-law model, as well as the reconstructed sources. Our model simultaneously reproduces the surface brightness structure of the lensed AGN and host galaxy in all bands. The normalized residual in the third column shows the area within the arcmask, as well as the region interior to the arcmask. In the IR/F160W band, there is an excess residual at the inner boundary of the arcmask (as well as outside of the arcmask, not shown in this figure) arising from the technical details of the PSF not being corrected outside of the arcmask. We run a test where the pixels showing excess residual outside of the arcmask are downweighted and find no significant change in the model parameters.

5.2. Composite model

Our composite model consists of a baryonic component linked to the light profile of the lens galaxy, plus a dark matter component. The composite model assumes the triple Chameleon light profile for the lens galaxy in the IR/F160W band scaled by an overall mass-to-light (M/L) ratio. The Chameleon light profiles link to parameters describing the light distribution to those of the mass distribution in a straightforward way, as they are fundamentally just a combination of isothermal profiles. We keep the triple Sérsic model for the lens galaxy light in the UVIS bands to maintain consistent parameterization with the power-law models. The dark matter component is modelled as an elliptical NFW (Navarro et al. 1996) halo with the centroid linked to the light centroid in the F160W band.

The fiducial composite model has the following free parameters in addition to the lens light parameters: (i) mass-to-light ratio (M/L) for the baryonic component, (ii) NFW halo scale radius, r_s, (iii) NFW halo normalization, κ_0, h (defined as ; Golse & Kneib 2002), (iv) NFW halo minor-to-major axial ratio, q_NFW, and associated position angle, φ_NFW, and (v) external shear, γ_ext, and associated position angle, φ_ext.

A Gaussian prior for the M/L of the baryonic component is employed, using the stellar mass constraint from Agnello et al. (2018) of for a Salpeter initial mass function (IMF). Although this value is lower than our estimate derived from the photometry of our models of the lens light profile (see Sect. 6), this prior has little influence on the result, as the model prefers an almost maximal M/L with little dark matter contribution (see Sect. 5.6). We set a Gaussian prior of based on the results of Gavazzi et al. (2007) for a sample of lenses in the SLACS survey (Bolton et al. 2006). These lenses span a redshift and velocity dispersion range that includes WGD 2038−4008, with a mean virial mass of . All other parameters are given uniform priors. The relative amplitudes of the three Chameleon profiles representing the stellar light distribution of the lens galaxy can vary within an MCMC chain. However, their relative amplitudes in the mass model initialization are necessarily fixed (due to the way that the GLEE user interface is set up), even though they share the same global M/L parameter. To account for this, we iteratively run a series of MCMC chains for the fiducial composite model and update the relative amplitudes of the three mass components to match that of the light components after each chain. After several iterations, the predicted Fermat potential stabilizes, and we stop iterating. We subsequently ran a test fiducial model using an updated version of GLEE in which the amplitudes of the mass components are directly linked to the light components and found that the results were unchanged. Figure 3 shows the data and the lens model results in all three bands for our fiducial composite model, as well as the source reconstructions.

Fig. 3. Same as Fig. 2, but for the fiducial composite model from GLEE.

5.3. Systematics tests

In this section we describe a variety of tests of the effects of various systematics in our modelling arising from different assumptions in the way we constructed the model that might affect the posterior. In addition to the basic fiducial models described above, we perform inferences for both the power-law and composite models given the following sets of assumptions: (i) a model where the regions near the AGN images are given zero weight rather than being scaled by a power-law weighting; (ii) a model where the region near the AGN images scaled by the power-law weighting is increased by one pixel around the outer edge; (iii) a model where the reconstructed source plane resolution in all bands is reduced to 40 × 40 pixels; (iv) a model where the reconstructed source plane resolution in all bands is increased to 60 × 60 pixels; and (v) a model with the arcmask region increased by one pixel on both the inner and outer edges. We combined the MCMC chains from all of these tests, weighted by the BIC (similar to Rusu et al. 2020, see Sect. 5.5).

5.4. Kinematics and external convergence

We used the kinematics and external convergence constraints from Buckley-Geer et al. (2020). We combined both LOS velocity dispersion measurements to constrain the lens models. Buckley-Geer et al. (2020) constrain the external convergence for different external shear amplitudes in steps of 0.01. For each model, we use the distribution corresponding to the external shear that is closest to the median amplitude for that model. We use importance sampling (e.g. Lewis & Bridle 2002) to simultaneously combine the velocity dispersion and external convergence distributions in a manner similar to Wong et al. (2017) and Rusu et al. (2020). For each set of lens parameters ν from our lens model chain, we draw a κ_ext sample from the distributions in Buckley-Geer et al. (2020) and a sample of r_ani from the uniform distribution [0.5, 5]θ_eff (θ_eff is calculated from the lens light distribution in the IR/F160W band from the power-law model). From these together with the D_ds/D_d ratio (that is fixed given the fixed Ω_m value of 0.3 in flat ΛCDM), we can compute the kinematics likelihood for the joint sample {ν, Ω_m, κ_ext, r_ani} via Eq. (21) and use this to weight the joint sample. We can then combine the Fermat potential computed from our lens model parameters ν with values of κ_ext and D_Δt to predict the time delays as a function of H₀ (via Eqs. (5) and (14)).

5.5. BIC weighting

We weight our models using the BIC, defined in Eq. (26). We take N_data (the number of data points) to be the number of pixels in the image region across all three bands that are outside the fiducial AGN mask (so that we are comparing equal areas), plus eight (for the four AGN image positions), plus one (for the velocity dispersion). k (the number of free parameters) is taken to be the number of parameters in the model that are given uniform priors, plus two (for the source position), plus one (for the anisotropy radius to predict the velocity dispersion). (the maximum likelihood of the model from the MCMC sampling) is the product of the AGN position likelihood, the pixellated image plane likelihood, and the kinematic likelihood. The image plane likelihood is the Bayesian evidence of the pixelated source intensity reconstruction using the imaging data within the arcmask (which marginalizes over the source surface brightness pixel parameters and is thus the likelihood of the lens parameters excluding the source pixel parameters; see Eqs. (12) and (13) in Suyu & Halkola 2010) multiplied by the likelihood of the lens model parameters within the image plane region that excludes the arcmask. We evaluate the BIC using the fiducial weight image and arcmask, as the majority of the models were optimized with these.

We estimate the variance in the BIC, , by sampling the fiducial model with source resolutions of [47, 48, 49, 50, 51, 52, 53, 54, 56, 58, 60] pixels on a side (the 50 × 50 pixel case is just the original fiducial model), keeping the arcmask the same. Changing the source resolution in this way shifts the predicted time delays stochastically, but there is no overall trend with resolution, and the degree of the shifts are smaller than the scatter among the different models in the systematics tests we run. We calculate the BIC for each of these models with different source resolutions and take the variance of this set of models as . We find for the power-law models and for the composite models.

To avoid biases due to our choice of lens model parameterization, we split the samples into the power-law and composite models and calculate the relative BIC and weighting for each set separately, similar to Birrer et al. (2019) and Rusu et al. (2020). Specifically, we weight a model with a given BIC of value x by a function f_BIC(x), defined as the convolution

(31)

where BIC_min is the smallest BIC value within a set of models (power-law or composite), and h is a Gaussian centred on x with a variance of . The exponential term is a proxy to the evidence ratio. We follow the calculation of Yıldırım et al. (2020) in evaluating the convolution integral in Eq. (31). Once we weighted time delay distributions for the power-law and composite models, we combined these two with equal weight in the final inference.

5.6. Modelling results with λ_int = 1

The marginalized parameter distributions of the power-law model are shown in Fig. 4. We show the combined distributions of all power-law models where each model is given equal weight, as well as the BIC-weighted distribution. Figure 5 shows the similar parameter distribution for the composite models. The point estimates for the mass model parameters from the GLEE models are presented and compared with those from the LENSTRONOMY models later in Sect. 7.2. The reconstructed sources of each model are qualitatively very similar, which is an important consistency check of the two models.

Fig. 4. Marginalized parameter distributions from our power-law lens model results from GLEE. We show the combined results from our systematics tests (dashed red contours) with each model weighted equally, as well as the BIC-weighted model results (shaded red contours). The contours represent the 68.3% and 95.4% quantiles.

Fig. 5. Marginalized parameter distributions from our composite lens model results from GLEE. We show the BIC-weighted model (shaded blue contours) and the combined results from our systematics tests (dashed blue contours). The contours represent the 68.3% and 95.4% quantiles.

The power-law model has a steep mass profile slope of γ = 2.30 ± 0.01, but the parameters are consistent with the previous model of Shajib et al. (2019). The various systematics tests do not show substantial variation. The ‘island’-like feature in Fig. 4 comes from the model with a lower source plane resolution, but this model is downweighted by the BIC, so it does not affect our result. The centroid of the mass and light profiles are consistent to within , and the model is able to fit the quasar positions to an rms of .

The composite model fits the quasar positions to an rms of , slightly worse than the power-law model. We note that the dark matter component contributes a very small fraction of the mass (of order ∼1%) relative to the stellar component, which has a large mass-to-light ratio. While this may appear unusual, the stellar mass enclosed within the Einstein radius determined from stellar population synthesis (SPS) models fit to the imaging data assuming a Salpeter IMF is consistent with the total enclosed mass as constrained by the lensing. In Fig. 6, we show the circularly averaged convergence of both the power-law and composite models. The effective Einstein radii (at which ⟨κ(< r)⟩ = 1) of the two models agree to within less than one UVIS pixel (), which corresponds to ∼2 − 3%. At the Einstein radius, the composite model slope closely matches the slope of the power-law model. The magnitude of the external shear (γ_ext) required for the power-law and composite models differs, resulting in a difference in the external convergence (κ_ext) as determined by Buckley-Geer et al. (2020).

Fig. 6. Radial mass profiles of the central deflector constrained by the GLEE models. Top: circularly averaged convergence ⟨κ(< R)⟩ as a function of radius for the GLEE power-law model (red) and composite model (blue). The shaded regions represent the 1σ credible regions. The stellar (green) and dark matter (black) components of the composite model are plotted separately. The vertical dashed black lines mark the pixel size in the F160W band and the best fit Einstein radius. Bottom: ratio of average convergence of the composite model to that of the power-law model as a function of radius.

The relative BIC weightings of each model are provided in Table 1. The blinded distributions of Fermat potential differences are plotted individually for each model in Fig. 7. The un-blinded illustrations of the BIC-weighted distributions are provided later in Sect. 7.1. Notably, the power-law and composite model have predicted time delays that are offset by ∼13%, indicating a difference in the two models. Contributing to this difference is the larger κ_ext for the composite model. As a result, the combined constraint has a larger uncertainty, reflecting this difference. Without factoring in the different κ_ext distributions, the power-law and composite models would be offset by ∼8%.

Fig. 7. Model-predicted distributions of Fermat potential differences (blinded) for each of the GLEE models tested, with power-law models (top) and composite models (bottom).

6. Lenstronomy modelling

In this section we describe the LENSTRONOMY model setups and modelling results. The software package LENSTRONOMY (Birrer & Amara 2018; Birrer et al. 2021) is a publicly available lens modelling software⁴. In contrast with GLEE, the software LENSTRONOMY uses basis sets to reconstruct the flux distribution of the background source galaxy (Birrer et al. 2015). In this section we describe the specific model settings for LENSTRONOMY on top of the baseline models from Sect. 4, then present our modelling results, and lastly combine the lens models with the measured stellar kinematics and the estimated external convergence.

6.1. LENSTRONOMY specific model settings

We explain particular model settings related to the mass and light profiles of the deflector galaxy in Sect. 6.1.1, the source light profiles in Sect. 6.1.2, and the image region for likelihood computation in Sect. 6.1.3. We summarize the set of all the lens models combining these different settings in Sect. 6.1.4.

6.1.1. Mass and light profiles of the deflector galaxy

We simultaneously model the HST images from all three bands. We join the centroids of the triple Sérsic profiles across the three bands in the power-law model setup, and also the centroids of the triple Chameleon profiles in the composite model setup. We join the ellipticity parameters of the light profiles only between the two UVIS bands. We let the amplitudes I_eff, effective radii θ_eff, and the Sérsic indices n_s in the three bands independently vary to allow for a colour gradient.

In the composite model setup, we adopt a Gaussian prior with mean and standard deviation for the NFW scale radius r_s based on the measurements of Gavazzi et al. (2007) for a sample of SLACS survey lens systems (Bolton et al. 2006). Since the velocity dispersion and the redshift of the central deflector of WGD 2038−4008 fall within the ranges spanned by the SLACS lenses, such a prior is appropriate (Treu et al. 2006). Similar priors were also adopted in previous H0LiCOW and STRIDES analyses (e.g. Wong et al. 2017; Rusu et al. 2017; Shajib et al. 2020). Although the measurement by Gavazzi et al. (2007) are reported in the physical kpc unit, we use the same fiducial cosmology as Gavazzi et al. (2007) to recover the scale in the observable angular unit. We also impose a prior on the concentration parameter using the theoretical M₂₀₀ − c relation from Diemer & Joyce (2019) with an intrinsic scatter of 0.11 dex.

6.1.2. Source light profiles

We adopt a basis set of shapelets and one elliptical Sérsic profile to describe the flux distribution of the quasar host galaxy. The Sérsic profile describes the smooth component of the flux distribution of the host galaxy, and the shapelets account for the non-smooth features (Refregier 2003; Birrer et al. 2015). The number of shapelets n_shapelets depends on the maximum polynomial order n_max as n_shapelets = (n_max + 1)(n_max + 2)/2, and the spatial extent of the shapelets is characterized with a scale size ς. We model the quasar images as point sources on the image plane. We treat the positions of the quasar images as free parameters throughout the model optimization and MCMC procedures. The point source positions are constrained directly through the likelihood of the pixel-level flux values in the imaging data. The four image positions give six independent relative positional parameters. We chose the option within LENSTRONOMY to solve the lens equation to constrain six parameters out of the set of the mass model parameters from these six independent relative positional parameters⁵. These six mass model parameters then have ‘one-to-one’ correspondence with the sampled quasar image positions. Therefore, they are not treated as non-linear parameters anymore in the optimization and sampling procedures. For the power-law model, the six parameters chosen are the PEMD’s centroid RA and Dec, axis ratio q_m, position angle φ_m, Einstein radius θ_E, and the external shear angle φ_ext. For the composite model, the six parameters chosen are the NFW profile’s centroid RA and Dec, axis ratio q_NFW, position angle φ_NFW, density normalization ρ_s, and the external shear angle φ_ext.

We join the ellipticity parameters of the source Sérsic profiles across the three bands. The centroids of all the light profiles are also joint across the three bands. This centroid is set at the quasar position in the source plane that is constrained through solving the lens equations for the four image positions. The effective radii θ_eff, the Sérsic indices n_s, the shapelet scale sizes ς for different bands are independent of each other.

We treat n_max as a hyper-parameter and fix it for a particular model optimization. A minimum number of shapelet components is necessary to describe the complex features in the lensed arcs; however, too many shapelet components will fit the noise in the imaging data. Thus, striking a balance between these two scenarios is necessary when choosing the number of shapelet components. We adopt three choices for : {7, 11}, {8, 12}, {9, 13}.

6.1.3. HST image region for likelihood computation

We chose a circular aperture in each band encompassing the lensed arcs centred on the lens galaxy to compute the imaging likelihood. The radii of these apertures are hyper-parameters in the model. We take two sets of choices for : with r_ℒ. Some nearby objects (stars or smaller galaxies) are masked out if they fall within the likelihood computation region (see Figs. 8 or 9 for the shape and comparative size of the likelihood computation regions).

Fig. 8. Most likely LENSTRONOMY lens model and reconstructed image of WGD 2038−4008 using the power-law model. The top row shows, from left to right, the observed RGB image, the reconstructed RGB image, the convergence profile, and the magnification model. The next three rows show, from left to right, the observed image, the reconstructed image, the residual, and the reconstructed source for each of the HST filters. The three filters are F160W (second row), F814W (third row), and F475X (fourth row). All the scale bars in each panel correspond to 1″. The star symbol in the reconstructed source panels marks the position of the quasar host galaxy’s centroid.

Fig. 9. Most likely LENSTRONOMY lens model and reconstructed image of WGD 2038−4008 using the composite model. The top row shows, from left to right, the observed RGB image, reconstructed RGB image, the convergence profile, and the magnification model. The next three rows show, from left to right, the observed image, the reconstructed image, the residual, and the reconstructed source for each of the HST filters. The three filters are F160W (second row), F814W (third row), and F475X (fourth row). All the scale bars in each panel correspond to 1″. The star symbol in the reconstructed source panels marks the position of the quasar host galaxy’s centroid. In the magnification model, a central image is predicted due to a central core in the triple Chameleon light profile. However, this central image is highly de-magnified, with magnification 0.019 ± 0.02, and thus its presence cannot be ascertained in our imaging data.

6.1.4. Model choice combinations

Summarizing the above sections, we have the hyper-parameter choices (i) for the lens galaxy mass profile: power-law and composite; (ii) for the source light : {7, 11}, {8, 12}, and {9, 13}; and (iii) for the likelihood computation region radii : and . Taking a combination of these choices, we have 12 different model setups. We perform the optimization with the same models setups twice. These twin runs are different due to stochasticity in the PSF reconstruction and MCMC sampling procedures, and help us assess random errors. As a result, we have 24 different optimized models, on which we perform BMA. The light profiles from the deflector, the lensed light profiles from the quasar host galaxy, and the point sources at the quasar image positions form a linear basis set for reconstructing the observed HST imaging. As a result, the amplitudes of these profiles are linear parameters, as they can be obtained through a linear inversion for a sampled set of non-linear parameters that describe all the mass and light profiles. There are 206−281 linear parameters and 51−54 non-linear parameters in our models.

6.2. Modelling workflow

For each model setting, we reconstruct the PSF in each HST band. The reconstruction is initiated from a PSF estimate with a corresponding error map created from ∼4−6 bright stars within the HST image. The PSF reconstruction is carried out in multiple iterations with model optimization having fixed PSFs interlaced in between the PSF reconstruction iterations (see Birrer et al. 2019 for details, and also Chen et al. 2016 for a similar algorithm). There is an offset between the recorded coordinates between IR and UVIS images. After each iteration of PSF reconstruction, we re-align the coordinate system of the IR image with that of the UVIS images using the quasar positions (Shajib et al. 2019). Thus, we have a block of three operations constructing one unit of PSF reconstruction iteration: (i) IR-band image re-alignment, (ii) lens model optimization, and (iii) PSF reconstruction.

We optimized the lens model using the particle swarm optimization (PSO) method (Kennedy & Eberhart 1995), which is implemented in LENSTRONOMY. After the PSF reconstruction procedure, we performed MCMC sampling of the model posterior using EMCEE (Foreman-Mackey et al. 2013), which is an affine-invariant ensemble sampler (Goodman & Weare 2010). We chose the number of walkers to be eight times the number of sampled parameters. We run the chain for 10 000 steps. We check for convergence of the chain by manually inspecting that the median and standard deviations of the parameters within the walkers at each step has reach equilibrium for at least 1000 steps. We take the walker distribution from the last 1000 steps of the chain to be the model posterior.

We illustrate the best-fit model from the model setup with the lowest BIC value among the power-law and composite model families in Figs. 8 and 9, respectively.

6.3. Bayesian model averaging

We have 24 models that make up our set of models {S, M}, with each lens model family from M ≡ {powerlaw, composite} has 12 different hyper-parameter settings in S. We approximate the integral on the right-hand side of Eq. (25) as a discrete summation over S as

(32)

Here, ΔS_n can be interpreted as the model space volume that represents the model S_n. Although the models {S_n} differ from each other by discrete steps, an appropriately chosen expression for ΔS_n can account for sparse sampling from the model space, as we cannot adopt a sufficiently large number of models that are densely populated in the model space due to computational limitation. We use the BIC score of a model as a proxy for the model evidence p(O_img ∣ M, S_n). Thus, exp(−ΔBIC/2) acts as the evidence ratio and provides the relative weight between two models. The ΔS_n term is effectively an additional weighting on top of this BIC weighting (Birrer et al. 2019; Shajib et al. 2020). We take p(S_n) = 1 and therefore need to effectively implement the weighted sum of p(ξ ∣ O_img, M, S_n) in the right-hand side of Eq. (32) through sampling.

We tabulate the BIC values of the models in Table 2. The BIC values are computed from the maximum sampled likelihood in each MCMC chain. We estimate the sparsity of models {S_n} by taking the variance of ΔBIC between ‘neighbouring’ models that differ with each other by one step in only one setting (Shajib et al. 2020). We furthermore accounted for the numeric uncertainty in estimation of ΔBIC by taking the variance of ΔBIC between identical models that we have optimized twice. These twin runs produce slightly different posteriors – and thus BIC values – due to stochasticity in the PSF reconstruction, PSO, and MCMC sampling steps, similar to what was done in Birrer et al. (2019). Thus, our total variance in ΔBIC is

(33)

We compute that and . To implement the ΔS_n weighting through sampling, we first follow Birrer et al. (2019) to obtain the absolute weight W_n, abs of the nth model by convolving the ΔBIC with the evidence ratio function f(x) as

(34)

where the evidence ratio function f(x) is defined using the BIC difference as

(35)

Then, we obtain the relative weight W_n, rel by normalizing the absolute weights by the maximum absolute weight as

(36)

Finally, we combine the individual model posteriors following Eq. (32) as

(37)

In Sect. 6.4, we compare the two mass model families after performing the above model averaging procedure within each model family.

6.4. Lensing-only constraints on the Fermat potential difference

We constrain the image positions with uncertainty for the power-law model, and with uncertainty for the composite model. Given the longest predicted time-delay for this system, these precisions are well below the astrometric requirement of ∼002 uncertainty so that the astrometric uncertainty is subdominant to achieve ≤5% uncertainty in H₀ from this single system (Birrer & Treu 2019).

In Fig. 10 we compare between the model settings – namely, the choices for n_max and the size of the image likelihood computation region. All the posteriors within a particular lens model family (i.e. power law or composite) are statistically consistent with each other within 1σ.

Fig. 10. Comparison between the Fermat potential difference posteriors from different LENSTRONOMY model settings. The posterior for a particular setting is obtained by averaging over models that differ in the other model settings, but within a particular model family using the procedure described in Sect. 6.3. Top two rows: correspond to the power-law mass model families, and bottom two rows: correspond to the composite mass model families. The shaded regions illustrate the 68% credible regions. The posteriors are blinded by , where the ‘ref’ subscript refers to one of the compared models. The potential differences are consistent within 1σ between our adopted choices of model settings.

We compare the combined model posteriors between the power-law and composite models in Fig. 11. The Fermat potential differences deviate by 16−21% between the power-law and composite model setups. However, the MST-invariant quantity is consistent between the two model setups (see Eq. (42) of Birrer 2021 for the full definition of ξ_rad, and also Kochanek 2020). Thus, the difference in the Fermat potential from the two model setups can be interpreted as a manifestation of the internal MSD. We combine the stellar kinematics and estimated external convergence with the lens models to mitigate the internal MSD in Sect. 6.5.

Fig. 11. Comparison of the LENSTRONOMY lens model parameters and Fermat potential difference between the power-law (red) and composite (blue) mass models. The posteriors are obtained by averaging over all the model settings following the procedure described in Sect. 6.3. The parameter γ for the composite model is computed from circularly averaging the quantity 2 − [d log α(r)/d log r]r = θE. Some parameters are blinded as for p ∈ {γ,  Δϕ}, where the subscript ‘pl’ refers to the power-law model posteriors. The composite model-predicted γ is approximately 8% smaller than that from the power law. Consequently, the Fermat potential differences are smaller by approximately 16% for the composite model.

However, we first check for potential unphysical properties in our best fit composite models as the source of the large difference in the Fermat potential in Sects. 6.4.1 and 6.4.2. These checks were performed prior to un-blinding the models.

6.4.1. Halo properties in the composite model

Figure 12 illustrates the M₂₀₀ − c relation posterior for our system in comparison with the adopted prior; the median of the concentration posterior is consistent with the concentration prior within 1σ. Our combined posterior from the composite model setup provides the total halo mass and the total stellar mass is . The total stellar mass is obtained by doubling the enclosed mass within the half-light radius of corresponding to the F160W band. The projected dark matter fraction within the Einstein radius is . The total baryon-to-dark-matter fraction is , which is consistent with the upper limit set by the cosmic baryonic fraction 0.19 (Planck Collaboration VI 2020). In Fig. 13 we plot the azimuthally averaged convergence profiles for the power-law and composite models to illustrate the difference in the convergence slope at the Einstein radius θ_E. The inner region (≲02) of the triple Chameleon profile is flat unlike the singular centre in the power-law model. The flat or cored convergence profile at the centre gives rise to an inner critical curve in the image plane (Fig. 9). This core in the centre of the stellar mass distribution follows from the stellar flux distribution, as the radial flux profile from isophotal fitting also shows a stellar core. We fit a core Sérsic profile – that is defined by Eq. (2) in Dullo (2019) – to the azimuthally averaged light profile from our isophotal fitting to obtain the stellar core radius. We obtain 0.780 ± 0.004 kpc assuming a fiducial flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3. This stellar core radius is consistent with the core radii measured in local elliptical galaxies (0.64−2.73 kpc; Bonfini & Graham 2016; Dullo 2019). In Fig. 14, we illustrate the velocity dispersion profiles predicted by the lens model posteriors assuming isotropic orbit and a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹. The composite-model-predicted velocity dispersion profile decreases towards the centre by ∼20% due to the flattened mass profile. Such a large decrease in the velocity dispersion has not been observed in local massive ellipticals (e.g. Cappellari 2016; Ene et al. 2019). As a result, the composite lens model posterior suggests an inconsistency with kinematic observations of the local ellipticals. Interestingly, an inner critical curve in our composite model predicts a central image. The magnifications of the 5 images are A: , B: , C: , D: 4.8 ± 0.3, and central: 0.019 ± 0.002. The predicted appearance of the central image is invariant under the MST. However, the demagnified and potentially dust-extincted central image is indistinguishable in the present imaging data. Thus, we are unable to distinguish the two profile families on the basis of the presence or absence of the central image.

Fig. 12. Distribution of M200 and c200 parameters for the NFW halos in LENSTRONOMY composite model (blue shaded region). This distribution is averaged over all the model settings within the composite model family following the procedure described in Sect. 6.3. The 2 contours correspond to the 68% and 95% credible regions, respectively. The black solid line traces the theoretical prediction of the M200 − c200 relation at zd = 0.230 from Diemer & Joyce (2019) with the grey shaded region corresponding to the 68% confidence interval. We adopt this M200 − c relation as a prior in our analysis in addition to a M200 prior based on Sonnenfeld et al. (2018).

Fig. 13. Radial mass profiles of the central deflector constrained by the LENSTRONOMY models. Top: circularly averaged convergence ⟨κ(< R)⟩ as a function of radius from LENSTRONOMY power-law (red) and composite (blue) models. The stellar (green) and dark matter (grey) distributions in the composite model are also individually illustrated. The shaded regions encompass the 16th and 84th percentile of the sampled profiles for the corresponding case. The grey points illustrate the surface brightness profile fitted with isophotes as described in Sect. 6.1.1. The amplitude of the isophotal fit is normalized to match with the triple Chameleon profile (green shaded region) at θE for the purpose of this illustration. The vertical black dashed lines mark the pixel size in the F160W band and the best fit Einstein radius. Bottom: ratio of the circularly averaged convergence profiles between the composite and the power-law models. At the Einstein radius the convergence slope deviates by 16−21% between the two model setups.

Fig. 14. Line-of-sight velocity dispersion profile (circularized) corresponding to the lens model posteriors in the power-law (red) and composite (blue) models from LENSTRONOMY. Isotropic stellar orbits are assumed in computing these velocity dispersion profiles. The solid lines correspond to the median and the shaded regions encompass the 16th and 84th percentiles. The composite profile predicts a decrease in the velocity dispersion towards the centre due to the flattened core – which is not observed in local massive ellipticals (e.g. Cappellari 2016; Ene et al. 2019) – thus pointing to the atypicality of the posterior mass profile in the composite model.

6.4.2. Test of potential systematics from modelling choices and priors

We checked if our composite models are robust against potential biases from our particular model choices, for example the likelihood computation region and the prior on the NFW halo scale radius. We first optimized a lens model with the power-law mass profile and the triple Chameleon profile for the light instead of the triple Sérsic profile. We then took these best fit parameters for the triple Chameleon profile and fixed them in the test composite setup. We let the overall scaling of the baryonic mass and light distributions be free, thus effectively allowing for a free mass-to-light ratio (M/L). We adopted a halo mass prior for the NFW profile dependent on the stellar mass given by

(38)

where is the total stellar mass based on the SPS method assuming a Chabrier IMF (Sonnenfeld et al. 2018). Here, the parameters are μ_h = 13.11 ± 0.04, β_h = 1.43 ± 0.15, and σ_h = 0.23 ± 0.04. We measure the stellar mass from the total fluxes in the three HST bands. We fit the surface brightness profile of the deflector separately in three bands from large cutouts that fully contains the light distribution of the deflector (see Fig. 15). First, we subtract lensed arcs and the quasar images from these cutouts using our best fit power-law model. Then, we fit elliptical isophotes using the PHOTUTILS software (Bradley et al. 2020). The method photutils.isophote.Ellipse.fit_image() allows a convenient way to ignore the overlapping objects (i.e. stars and galaxies) by sigma-clipping pixels along an isophote. Thus, we do not need to mask out these overlapping objects. We reconstructed the surface brightness profile of the deflector based on the fitted isophotes, which effectively interpolates through the pixels that are contaminated by overlapping objects. We obtain the total flux in each HST band from the reconstructed surface brightness profile using the fitted isophotes. We used PYGALEXEV⁶ to obtain , which is a PYTHON wrapper for GALAXEV (Bruzual & Charlot 2003). By adopting the Basel Stellar Library (BaSeL; Lejeune et al. 1998), exponentially decaying stellar formation history, and free metallicity, we obtain . Thus, from Eq. (38), our Gaussian prior on the halo mass log₁₀ M₂₀₀ has mean 13.5 and standard deviation 0.3. We additionally adopt a prior for the total stellar mass log₁₀ M_⋆. We take an ad hoc prior that is uniform between 11.51 and 11.88 and drops off like a Gaussian function outside these limits. The range between 11.51 and 11.88 accounts for the unknown IMF and spans the range between light (e.g. Chabrier) and heavy (e.g. Salpeter) IMFs that differ by ∼0.25 dex in the stellar mass. The exact form of this prior is

(39)

Fig. 15. Fitting the lens galaxy light profile from a much larger cutout than that used for lens modelling. Left: large cutout around the deflector galaxy in the F814W band, with the lensed arcs and quasar images subtracted using the best-fit of a power-law lens model. The galaxy extends much further beyond the Einstein radius with θeff/θE ≈ 1.7. Middle: model for the surface brightness profile of the deflector constructed using elliptical isophotes. Right: model residual showing that the model has captured the overall light distribution despite numerous overlapping objects. The model from middle panel was further approximated with an elliptical multi-Gaussian expansion (MGE; Cappellari 2002) to allow deprojection along the LOS for kinematic modelling.

where , A is the amplitude that normalizes the probability distribution to have . The actual value of A is not required for sampling in the MCMC method. We also allow an additional uncertainty of ±0.06 dex in M_⋆ to allow 15% uncertainty on the assumed H₀ in the SPS-based stellar mass estimation.

We furthermore mask out the central region in the deflector galaxy in the test composite model setup so that the optimization does not incentivize the presence of a central image to make up for residuals in the deflector light galaxy model that cannot be fully accounted by the triple Chameleon light profile.

We perform the 12 different model setups also for this test composite model and combine the posteriors based on their BIC values. We compare our primary composite model with the test composite model in Fig. 16. Although the Fermat potential differences are consistent between these two model setups within 1σ, the ones from the test setup are smaller by 4−7% than the primary setup. For the test setup the stellar mass is , the halo mass is , the total baryonic fraction is , and the dark matter fraction within the Einstein radius is . The total halo mass increases in the test setting over the one obtained from our primary setting due to the adopted halo mass prior. As a result, the increase in the halo normalization leads to a decrease in the Fermat potential, or equivalently the shallowing of the convergence profile. This impact of the prior on the halo profile normalization is the same as observed by Shajib et al. (2021) for the SLACS lenses. As adopted priors can systematically shift Fermat potential differences, more physically motivated priors are not sufficient to explain all the differences between the composite and the power-law models.

Fig. 16. Comparison of the Fermat potential posteriors between two different model settings for the LENSTRONOMY composite model. The blue solid distributions correspond to our primary model settings. In the test model setup (black dashed distributions), we fix the triple Chameleon profile for the stellar mass and light distributions to best fit values from a separately optimized model with the power-law mass profile. We also mask the central region of the deflector galaxy in the test model setup. Additionally, the prior on the halo mass profile is different. Whereas we adopt a prior on the NFW scale radius rs in the primary setup, we adopt a combination of priors on M200 and M⋆ in the test setup. Despite multiple differences in the model settings, the Fermat potential differences are consistent within 1σ with each other.

The mass difference of 3.8 × 10¹⁰ M_⊙ within (assuming flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹) between the power-law and composite models could be explained by an ultra-massive black hole (e.g. Mehrgan et al. 2019; Dullo 2019). Another potential solution that would push the Fermat potential differences from the composite models towards the ones from the power-law models is to have stellar mass-to-light ratio gradient (η ∼ 0.27 with M_⋆/L ∝ R^−η) to steepen up the total mass density profile. Although Shajib et al. (2021) find no evidence for a significant mass-to-light ratio on average for SLACS lenses at the similar redshift (⟨z⟩∼0.2) as WGD 2038−4008 (z_d = 0.23), there can still be individual cases with steep mass-to-light ratio gradients. Moreover, the value η ∼ 0.27 will be consistent with the constraints ⟨η⟩ = 0.24 ± 0.04 reported by Sonnenfeld et al. (2018). Adopting a mass-to-light ratio gradient for the luminous component or including a central black hole in the composite lens model is beyond the scope of this study. At this point, however, there is no a priori reason to modify the composite model to push the Fermat potential differences towards the ones from the power-law models. We use the kinematics data to bridge the discrepancy between the two models or to be the decider between them next in Sect. 6.5. As the composite models are related to the true underlying mass distribution through an approximate MST, we rely on the kinematics data to appropriately adjust the Fermat potential differences along the MSD towards the true values.

6.5. Combining with stellar kinematics and external convergence

In this subsection we perform dynamical modelling of the deflector based on our lens models using Eq. (17). We need to estimate the luminosity density l(r) to use in this equation by deprojecting the surface brightness profile of the deflector. The surface brightness of WGD 2038−4008 extends far beyond the size of our likelihood computation region (Fig. 15). Thus, deprojecting the light profile fit from our lens models may potentially produce early truncation in the 3D luminosity distribution along the LOS. Therefore, we use the reconstructed surface brightness profile from the fitted isophotes in the F814W band from Sect. 6.1.1 (see Fig. 15). We chose the light profile from the F814W band, because the velocity dispersion was measured in the optical. From this reconstructed surface brightness profile, we numerically find the circular aperture that contains half of the total light as . The aperture size in this numeric computation can only grow by a size of a pixel, which is . We adopt a 3% Gaussian uncertainty for θ_eff, which combines one pixel size as a systematic uncertainty with a typical 2% uncertainty for θ_eff from fitting surface brightness profiles with continuous parameters from Shajib et al. (2021). We check that adopting 20% uncertainty for θ_eff or using the θ_eff measured from the F160W band does not significantly impact (≲0.1%) the resultant Fermat potential difference (Fig. 17). We approximate the reconstructed surface brightness profile with an elliptical 2D multi-Gaussian expansion (MGE) series (Cappellari 2002). The MGE approximation allows for a straightforward deprojection into a 3D light profile, which we use as l(r) in Eq. (17). Since we only solve the Jeans equation in the spherical case, we adopt the spherical equivalent of the elliptical Gaussian components by taking the Gaussian scales along the intermediate axes. We also apply a self-consistent 2% uncertainty to the MGE scale parameters by letting them vary with θ_eff. We adopt a uniform prior on log a_ani, where a_ani is the anisotropy scaling factor defined by r_ani ≡ a_aniθ_eff. Birrer et al. (2016, 2020) find that the uniform prior for log a_ani is a less informative choice than the uniform prior on a_ani.

Fig. 17. Checking for systematics in the kinematic computation. Our primary settings for kinematic computation adopts the F814W light profile with 3% uncertainty in θeff (red). The orange contours are for the case with 20% uncertainty in θeff, and the black contours are for the case with F160W light adopted in the kinematic computation. The vertical blue line marks the measured LOS velocity dispersion, which has an uncertainty of 19 km s−1. The difference in the computed velocity dispersion between these cases is negligible, and it impacts the Fermat potential differences by ≲0.1%.

We performed a test for systematics in our velocity dispersion modelling. We adopted two test cases: (i) where the θ_eff uncertainty is taken as 20%, and (ii) where the light profile from the F160W band is used in the kinematic computation.

For the external convergence κ_ext, we impose a selection criterion on the P(κ_ext) estimated in Buckley-Geer et al. (2020) by requiring that the selected LOSs also correspond to the combined (through BIC weighting) external shear value from our lens models within a lens model family (i.e. power law or composite). In Fig. 18 we illustrate the two κ_ext distributions consistent with the external shear values for the power-law and composite mass profiles⁷.

Fig. 18. External convergence distribution from Buckley-Geer et al. (2020) with additional weighting applied based on the predicted external shears for the power-law (red lines) and composite (blue lines) models from LENSTRONOMY (solid lines) and GLEE (dashed lines). Each illustrated GLEE distribution is a BIC-weighted combination of multiple κext distributions corresponding to the external shear constraint from individual lens model setups. In contrast, each illustrated LENSTRONOMY distribution is a single κext distribution corresponding to the combined (through BIC weighting) external shear value from all the model setups within a model family (i.e. power law or composite). The κext distributions used by one team were not revealed to the other team before the un-blinding to maintain independence.

We combine the stellar kinematics and external convergence information with the lens model posterior in two different ways: (i) with fixed λ_int = 1 (Sect. 6.5.1), and (ii) with free λ_int constrained by the stellar kinematics (Sect. 6.5.2).

6.5.1. The case with λ_int = 1

Assuming λ_int = 1, the model-predicted velocity dispersion can be written as

(40)

This assumption when combining the stellar kinematics information with the lens model posteriors is the same as done in earlier TDCOSMO analyses prior to Birrer et al. (2020, TDCOSMO-IV), for example in Suyu et al. (2013), Wong et al. (2017), and Rusu et al. (2020). We assume a flat ΛCDM cosmology with Ω_m = 0.3 to compute the fiducial distance ratio D_s/D_ds. To combine the stellar kinematic information, we consider the kinematics likelihood function

(41)

We first combine the lens model posteriors from the power-law and composite models with equal weights, and then importance sample from this combined posterior with weight ℒ_kin (Lewis & Bridle 2002). We note that each of the power-law and composite models are already averaged over the various adopted model settings within each mass model family following our BMA procedure from Sect. 6.3. We illustrate the Fermat potential differences from each of the power-law and composite models in Fig. 19.

Fig. 19. Fermat potential differences Δϕ = Δϕmodel(1 − κext) from LENSTRONOMY for the power-law (red) and composite (blue) model with the kinematics information folded into these individual model posteriors. The folding in of the kinematics information is performed through importance sampling from the posterior weighted by the kinematics likelihood (Eq. (41)). Next, the two model posteriors are joined together with equal weights and then the kinematics information is folded in. The combined posterior (purple) mostly resemble the power-law posterior, as the composite posterior is heavily down-weighted by the kinematics likelihood.

After combining the kinematics information with λ_int, the combined posterior for the Fermat potential differences end up mostly similar to the power-law posterior, as the kinematic likelihood heavily down-weights the posterior from the composite model. Although the composite model was designed with physical motivations to mimic a real galaxy structure, the kinematics data heavily disfavours the composite lens model posterior for λ_int = 1. Furthermore, applying more physical priors to resolve this inconsistency rather makes the kinematics data disfavour the composite model more, which suggest that our composite model is not adequate in describing the true galaxy mass distribution. Further generalization in the composite model (e.g. mass-to-light ratio gradient and generalized NFW halo) may thus be necessary for a composite model to be simultaneously consistent with the lensing data, the kinematics data, and the cosmological expectations for galaxies, e.g. the M − c relation, baryonic fraction. The uncertainties on the combined Fermat potential differences, and thus on the predicted time delays, are approximately 4%, which is comparable with those from the previous TDCOSMO analyses under the same assumption of λ_int = 1 (Wong et al. 2020).

6.5.2. The case with free λ_int

Now, we treat λ_int as a free parameter and constrain it using the stellar kinematics by re-expressing Eq. (22) as

(42)

Such constraining of λ_int from stellar kinematics is the same approach as Birrer et al. (2020, TDCOSMO-IV), albeit these authors achieved a tighter constraint on λ_int from a joint sample of seven time-delay lens systems and 33 non-time-delay lenses through a hierarchical Bayesian analysis.

We obtain the D_s/D_ds distribution to use in Eq. (42) from the relative distance constrained by the Pantheon SN sample (Scolnic et al. 2018). We approximate the luminosity distance up to the Pantheon supernovae using a fourth-order Taylor expansion. The coefficients in the Taylor expansion allow increasing complexity by including the deceleration parameter q₀, the jerk parameter j₀, the snap parameter s₀, and the curvature density parameter Ω_k. We compute the model evidence using nested sampling for different size of the parameter set (Skilling 2004). We select the model that goes up to the jerk parameter j₀ based on its highest model evidence. We use the relation D_A = D_L/(1 + z)² to convert the luminosity distance to angular diameter distance and transform the 2D posterior distribution of (q₀, j₀) to obtain the D_s/D_ds distribution.

In Fig. 20 we compare the MST-corrected Fermat potential differences between the power-law and the composite models. We find for the power-law model and for the composite model. This large median value of falls in the excluded region in Birrer et al. (2020, TDCOSMO-IV) that is based on physical arguments on the mass density distribution. However, the mass profile adopted by Birrer et al. (2020) is a cored power-law profile, which can be interpreted as the presence of a cored component in the NFW profile with its radius being much larger than the Einstein radius. Shajib et al. (2021) demonstrate that deviations from a power-law profile can be explained by shifting the normalization of the NFW profile without any core component. Whereas the MST considered by Birrer et al. (2020, TDCOSMO-IV) allows redistribution of matter only within the dark component, the composite model considered here allows redistribution of matter between dark and luminous components. Thus, large deviations of λ_int from 1 is still physically plausible, so the large λ_int produced by our composite model is not in tension with the exclusion range set by Birrer et al. (2020, TDCOSMO-IV).

Fig. 20. MST-adjusted Fermat potential differences Δϕ = Δϕmodelλint(1 − κext) from LENSTRONOMY for the power-law (red) and composite (blue) model families. The internal MST parameter λint is estimated by combining the lens models with the measured velocity dispersion and the external convergence estimate. After adjusting for the MST, the power-law and composite model predictions for the Fermat potential differences become consistent with each other within 1σ.

The predicted Fermat potential difference between the power-law and composite models are consistent within 1σ after adjusting for the internal MST and the external convergence. Thus, we demonstrate that the two model families we adopted are linked through an approximate MST. Therefore, to predict the time delays or to measure H₀ through constraining λ_int from stellar kinematics, the choice of mass model family is largely irrelevant as the same result can be obtained with any of the conventional model families. In this case, the uncertainty on the Fermat potential difference is dominated by the velocity dispersion uncertainty, which is at 6.4%. As , the uncertainty on the λ_int is thus expected to be twice the uncertainty of the velocity dispersion. Our obtained uncertainties of and are consistent with this expectation. The uncertainties on the predicted time delays from the final combined posterior with free λ_int is 19−22%.

6.6. Discussion on LENSTRONOMY models

The predicted Fermat potential differences from LENSTRONOMY are discrepant between the power-law and composite models, with the predictions from the composite model being 16−21% lower than those from the power-law model. Both model families fit the imaging data almost equally well, with the composite model providing a slightly higher likelihood value (Table 2). The discrepancy between the two model families are caused by the NFW profile normalization in the composite model, which makes the logarithmic slope of the density profile shallower at the Einstein radius (Fig. 13). We check for any potential unphysical properties in the mass profile posterior for the composite model (see Sect. 6.4.1). The size of the central core observed in stellar mass profile is consistent with previous observations (e.g. Bonfini & Graham 2016; Dullo 2019). The halo properties (e.g. the M − c relation, the total baryonic fraction, and the dark matter fraction within the Einstein radius) are consistent with previous observations of galaxy properties and cosmology. However, the predicted velocity dispersion profile (Fig. 14) from the composite mass model shows a decrease towards the centre, which has not been observed in local massive elliptical galaxies (e.g. Cappellari 2016; Ene et al. 2019), thus pointing to a potential inconsistency in the composite profile. We tested with additional physically motivated priors for the halo mass profile; however, that only amplified the discrepancy further. Alternatively, the discrepancy between the power-law and composite mass profiles can be reconciled by including an ultra-massive black hole (M_BH ∼ 3.8 × 10¹⁰ M_⊙), or by incorporating a stellar mass-to-light ratio gradient with exponent η ∼ 0.27, or a combination of both. These values are plausible based on previous observations (e.g. Sonnenfeld et al. 2018; Mehrgan et al. 2019; Dullo 2019).

Furthermore, the predicted central velocity dispersion from the composite model is inconsistent with the observed one. As a result, when no internal MSD is assumed (i.e. λ_int = 1), the kinematics likelihood largely excludes the posterior from the composite model in the final combined posterior. As a result, the final combined posteriors from LENSTRONOMY is almost entirely contributed by the power-law model.

7. Comparison of the two software programs

The lens model posteriors from both teams were un-blinded on October 22, 2021, and no further modification to the lens models was performed afterwards. We only performed tests to investigate the differences or the lack thereof between the two modelling teams; the final time-delay predictions are kept frozen at the values during un-blinding. Table 3 compares the model parameters, derived quantities, predicted time delays between GLEE and LENSTRONOMY. We compare the un-blinded time-delay predictions from the combination of lensing, kinematics, and LOS analyses in Sect. 7.1. Then, we compare the lens model parameters and Fermat potential differences from lens modelling only in Sect. 7.2. We compare the pixelized PSF reconstructions between the software programs in Sect. 7.3 and the computational requirements in Sect. 7.4. Finally, we discuss our findings in Sect. 7.5.

7.1. Predicted time delays

We illustrate the final predicted time delay from both teams as a function of H₀ in Fig. 21, assuming a flat ΛCDM cosmology with Ω_m = 0.3. The predictions for all image pairs are consistent with each other within ∼1σ. We further compare the time delay predictions for combined, power-law-only, and composite-only cases in Fig. 22. The combined time-delay posteriors differ the largest for the AB image pair by 11% (1.2σ). We note that the GLEE team applied kinematics weighting to the lens model posteriors only within the mass families and then combined the mass families with equal weighting, whereas the LENSTRONOMY team weighted the mass families by kinematics. As a result, the combined posteriors from the GLEE team receives equal contribution from both model families giving rise to bi-modalities, where the combined posterior from the LENSTRONOMY team is almost entirely contributed by the power-law model posterior. The power-law-model-predicted time delays agree better between the two teams, with the largest difference appearing for AB image pair by 3.4% (0.6σ). The composite model predictions are more discrepant, with the largest difference appearing for AC image pair by 15% (2.1σ).

Fig. 21. Comparison between the two modelling teams for the predicted time delays (un-blinded) as a function of H0 for the three image pairs, in flat ΛCDM cosmology with Ωm = 0.3. Each posterior is the final combined posterior from the two mass-model setups: power-law and composite, including external convergence and stellar kinematics, and assuming λint = 1.

Fig. 22. Comparison of predicted time delays between the LENSTRONOMY (green) and GLEE (purple) teams. The three rows show the comparison for combined (top), power-law-only (middle), and composite-only (bottom) cases. The combined posteriors are consistent within < 1.2σ, the power-law model’s posteriors are consistent within < 0.6σ, and the composite model’s posteriors are consistent within < 2.1σ. In percentage, the maximum deviation for the power-law model is between AB image pair by 3.4%, and for the composite model is between AC image pair by 15%. The GLEE team combined the power-law and composite model posteriors with equal weighting after applying the kinematics weighting within each model family, whereas the LENSTRONOMY team combined the two model families with kinematics weighting, leading to the final combined posterior being dominated by the power-law model.

7.2. Model parameters and Fermat potential differences

The astrometric uncertainty on the constrained AGN image positions are consistent between the two teams, as both teams constrain the image positions with uncertainty at maximum. This astrometric precision satisfies the requirement for cosmographic measurements (Birrer & Treu 2019).

In Fig. 23 we compare the lens model parameter posteriors and the predicted Fermat potential posteriors from the power-law lens models from both teams. The power-law exponent constrained by the LENSTRONOMY model is γ = 2.21 ± 0.02, whereas that constrained by the GLEE model is γ = 2.30 ± 0.01, which is a discrepancy at 4σ. The external shear magnitude constrained by the LENSTRONOMY model is γ_ext = 0.065 ± 0.004, and that by the GLEE model is γ_ext = 0.078 ± 0.001, which is at a 3.4σ discrepancy. We identify a degeneracy between γ and γ_ext internal to both models, and the discrepancy between the posteriors from the two models lie along this degeneracy (Fig. 23). The external shear magnitudes are typical for quadruply lensed quasar systems (e.g. see Schmidt et al. 2022). We investigated for deviation from the simple ellipticity description in our mass models as a potential source of the external shear. We find that boxy-ness or discy-ness in the luminous component is negligible within the Einstein radius (i.e. ), which implies that allowing boxy-ness or discy-ness in the description of the mass profile is not required (for definitions of a₄ and b₄, and their impact on H₀ measurement, see Van de Vyvere et al. 2022a). As a result, we attribute the LOS galaxies around the central deflector to be the main source of the external shear, with additional potential contribution from the mild isophotal twist beyond the Einstein radius in the central deflector (Van de Vyvere et al. 2022b).

Fig. 23. Comparison of lens model parameter differences for the power-law model between the LENSTRONOMY (green) and the GLEE (purple) teams.

Next in Fig. 24, we compare the Fermat potential differences only from the power-law models of both teams without adjusting for the external convergence (top row) and with adjustment for the external convergence (bottom row). The Fermat potential differences from the lens model are discrepant, for example by 5.5σ (8.9%) for the AD image pair that has the longest predicted time delay. However, after combining the corresponding external convergence – based on selection cuts using the best fit external shear from each model – the Fermat potential differences all become consistent within 1σ, for example by 0.26σ (1%) for the AD image pair. Interestingly, the positive correlation between γ and γ_ext allows the estimated κ_ext selected on γ_ext to bring the time-delay posteriors closer, as higher γ_ext selects higher κ_ext. However, the strength of the positive correlation between γ and γ_ext depends on the particular morphology of the quad lenses (Shajib et al. 2019). Thus, we cannot conclude if this effect – that the γ_ext-selected κ_ext brings the time delay posteriors more into agreement – applies to all lensing systems and is not just a particular occurrence for the lens system WGD 2038−4008. A detailed analysis of a larger sample of lenses is required to reach a conclusion on this matter, and it is left for future work.

Fig. 24. Comparison of Fermat potential differences for the power-law model between the LENSTRONOMY (green) and the GLEE (purple) teams, without including the external convergence (top row), and with including the external convergence (bottom row).

7.3. Reconstructed PSFs

The reconstructed pixelized PSFs by the GLEE modelling procedure have smaller FWHM by ∼2−7% than the ones from LENSTRONOMY: 1.7% in F160W, 3.9% in F814W, and 6.7% in F475X (Fig. 25). Furthermore, the F160W PSF in GLEE is supersampled with a supersampling factor of 3, whereas the F160W PSF in LENSTRONOMY has the same pixel resolution () as the drizzled image. The PSFs for the UVIS filters have the same pixel resolution as the drizzled image () for both GLEE and LENSTRONOMY.

Fig. 25. Comparison of the reconstructed pixelized PSFs by the LENSTRONOMY (first column) and the GLEE (second column) modelling procedures. Third column: illustrates the difference between the LENSTRONOMY and GLEE PSFs. The three rows correspond to the F160W, F814W, and F475X filters from top to bottom. The F160W PSF is supersampled with a supersampling factor of 3. We note that the illustrated supersampled LENSTRONOMY PSF was not used in the pre-un-blinding models and was reconstructed after the un-blinding to perform further tests. The original reconstructed PSF in LENSTRONOMY in F160W was not supersampled. The GLEE PSF FWHMs are smaller by ∼2−7% than the LENSTRONOMY ones.

We test how the differences in the reconstructed PSFs contribute to the differences in the logarithmic slope parameter for the power-law model between the two software programs. The test results are illustrated in Fig. 26. In these tests, we change the adopted PSFs and optimize a fiducial lens model from each software program. For LENSTRONOMY, the fiducial model is the power-law model with the highest BIC score (Table 2). For GLEE, the fiducial model is the ‘power-law fiducial model’ (Table 1). We first interchanged all the PSFs between the software programs. Due to numerical requirements, the GLEE team artificially supersample the LENSTRONOMY-reconstructed PSF through interpolation for this and subsequent tests. With the PSFs interchanged, the power-law slope parameter γ constrained by one software program shifts towards the fiducial constraint from other software program. These shifts bring the γ distributions within ∼1.2−1.6σ consistency between the two software programs given the same PSFs, although still leaving some unexplained deviations.

Fig. 26. Deviations in the logarithmic slope γ of the power-law model with different PSF settings. In all panels, the dashed distributions show the fiducial constraints – LENSTRONOMY in green and GLEE in purple. The LENSTRONOMY fiducial constraint is from the highest BIC value power-law model (Table 2), and the GLEE fiducial constraint is from the ‘power-law fiducial’ setup (Table 1). First panel: when only the reconstructed PSFs are interchanged to optimize the models, the constraint from one software program moves towards the other’s fiducial constraint to be consistent within 1.2 − 1.6σ. We note that GLEE model artificially creates supersampled version of the LENSTRONOMY PSF through interpolation to perform the tests here. Second panel: when both the PSFs and weighted uncertainty maps are interchanged, the resultant γ becomes slightly more consistent between the software programs within 1.0 − 1.5σ. Therefore, we conclude that the differences in the PSF itself significantly contributes to the discrepancy in the power-law γ constraint and not the particular method of weighting the uncertainty map to account for PSF uncertainty. Third panel: when models from both of the software programs are optimized only with UVIS data and interchanged PSFs, the GLEE constraint does not shift significantly; however, the LENSTRONOMY constraint shifts significantly towards the GLEE constraints. This result indicates that the difference between the GLEE and LENSTRONOMY fiducial models are largely created by the difference in the IR PSF. Fourth panel: when both software programs use a supersampled PSF with a supersampling factor of 3 reconstructed by LENSTRONOMY, the resultant γ constraint agrees very well, within 0.6σ.

We further interchange the weighted uncertainty maps in addition to the reconstructed PSFs between the software programs. Originally, the GLEE team creates a weighted uncertainty map by boosting the noise levels around the quasar positions (see Sect. 5), whereas the LENSTRONOMY team adds the PSF uncertainty map of the initial PSF estimate in quadrature with the data uncertainty map at the positions of the quasars. In this test, the resultant γ distributions become slightly more consistent within ∼1.0−1.5σ compared to the previous test. Therefore, the particular choice or method to estimate the weighted uncertainty maps does not significantly contribute to the deviation in the power-law γ parameter distribution from the two software programs.

In the next test, we constrained the lens models from UVIS data only with interchanged PSFs. In this test, the GLEE constraint on γ remained stable; however, the LENSTRONOMY constraint shifted towards the GLEE constraint to be consistent within ∼0.5σ. As a result, we conclude that the discrepancy in the γ distribution between LENSTRONOMY and GLEE is dominated by the difference in the IR PSF.

For the last test, we add a new feature in LENSTRONOMY to reconstruct the IR PSF with a supersampled resolution and reconstruct the PSF with a supersampling factor of 3 following the GLEE team. If this supersampled IR PSF from LENSTRONOMY is used by both software programs, then the resultant γ distribution becomes consistent within ∼0.6σ. However, the γ distribution of the GLEE model with its own supersampled PSF still differs by 2.6σ from that of the LENSTRONOMY model with its own supersampled PSF. It is not possible to evaluate which reconstructed PSF is more accurate a priori. Therefore, it is recommended to marginalize over multiple PSF reconstructions to account for the stochasticity within one particular reconstruction algorithm and as well as different reconstruction algorithms. In addition, supersampled PSFs are recommended especially for the IR band with large pixels; the subsampling factor can be set to the minimal value to produce stable results while keeping computational time low.

7.4. Requirements for computational resources

The entire LENSTRONOMY modelling procedure required ∼4 × 10⁵ CPU hours including initial modelling trials, running full MCMC chains of the adopted models, post-processing of posteriors, and post-un-blinding tests. The estimated usage for GLEE models are 𝒪(10⁵) CPU hours only to run the final models and MCMC chains to produce the final un-blinded result. However, the total usage of CPU hours can be 𝒪(10⁶) CPU hours including modelling trials, robustness tests, and reruns of chains to recover lost progress due to numerical issues.

7.5. Discussion

The final un-blinded time-delay predictions agree within < 1.2σ between the two modelling teams. As a result, the inferred Hubble constants from the two teams based on the observed time delays will be consistent within < 1.2σ. However, the predictions from the composite models only are less in agreement between the two teams. Interestingly, the composite-model predictions deviate from the power-law ones towards the same direction for both teams. We were already aware prior to the un-blinding that the composite model is atypical with respect to previous observations, for example the velocity dispersion profile significantly decreases towards the centre for the one from the LENSTRONOMY (Fig. 14), the one from the GLEE team has a very low inner dark matter fraction. Although such discrepancies between composite and power-law model predictions have been observed in the previously analysed systems (e.g. Suyu et al. 2014), this system WGD 2038−4008 demonstrates the largest discrepancy to date out of the systems analysed by H0LiCOW/TDCOSMO. However, unlike the previously analysed systems, the Einstein radius of this system encompasses only the very central region of the very extended lens galaxy, and thus the imaging observables probe a different region of the elliptical galaxy: at ∼θ_eff/3 instead of ∼θ_eff. This discrepancy in the composite model illustrates that this model is not an adequate description for the mass distribution at the central region of all elliptical galaxies. Rather an appropriate combination of a mass-to-light ratio gradient and a supermassive black hole can be necessary to sufficiently describe the mass distribution at the scales considered here. In such cases when the different mass model families fit the imaging observables equally well, but lead to different predictions in the Fermat potential and kinematics, the observed kinematics is crucial to act as the differentiator between the mass model families through appropriate weighting of the kinematics likelihood. In the future, spatially resolved velocity dispersion from integral field spectra, for example from the Multi Unit Spectroscopic Explorer (MUSE) on the VLT, will be able to constrain such an improved composite model with additional degrees of freedom allowed.

For the composite model setup, the modelling teams had more freedom in choosing particular priors and model settings, allowing for discrepancies between the results from the two teams. However, the power-law models are specified with less room for independent choices to be made by the modelling teams. Therefore, we compare the results from the power-law models between the two teams to identify systematic differences at the level of the software packages.

In particular, we focus on the power-law logarithmic slope parameter γ, as this parameter is the most sensitive lens model parameter to the predicted time delays and thus the inferred Hubble constant. The γ distributions between the modelling teams are discrepant at 4σ. We identify the difference in the reconstructed PSFs, especially in the IR band, to be the dominant source of this discrepancy. Given the same supersampled PSF in the IR band and non-supersampled PSFs in the UVIS bands, both modelling softwares produce γ constraints with differences below 0.5%. The γ distributions from the two modelling software programs are not expected to be identical due to differences in the numeric implementation, for example the likelihood computation region and the source reconstruction method. Thus, we can conclude that the systematic differences in the model fitting part of the software programs are below 0.5% in γ corresponding to ∼1% on H₀. However, the PSF reconstruction method can lead to systematic differences as large as ∼4%. This difference is reconciled in the time-delay predictions from the two teams after combining the lens model posteriors with the kinematics data and the estimated external convergence for the particular system WGD 2038−4008. However, it is inconclusive if such differences can be similarly reconciled for all lens systems in general, and we leave this investigation with a larger lens sample for a future study. If such differences cannot be reconciled for other lens systems and these differences do not average out when a sample of lenses is considered, then these differences would be non-negligible in the long run when a large sample of lens systems are combined to infer the Hubble constant. As it is currently not possible to evaluate the appropriateness of one reconstructed PSF over the other, we recommend to marginalize over different realizations of the PSF reconstruction and also over different reconstruction algorithms to avoid any potential bias. Furthermore, a supersampled PSF in the IR band is also recommended, as the drizzled pixel scale of in the IR does not optimally sample the PSF.

8. Summary and conclusion

In this study two teams independently modelled the lens system WGD 2038−4008 from three-band HST imaging using two different software programs: LENSTRONOMY and GLEE. Two families of models were specified as baseline models before the modelling was performed – one family with a power-law ellipsoidal mass distribution and the other family with a two-component mass distribution that accounts for the dark and baryonic components separately. The baseline settings were pre-specified to allow a fair comparison between the modelling results. Individual teams were allowed to improve upon the baseline settings as they deemed appropriate, for example the choice of priors and the numerical settings specific to the software program being used. The two modelling procedures were carried out blindly with regards to the other team. The models were un-blinded on October 22, 2021, after an internal review by scientists from the TDCOSMO collaboration not directly involved with either modelling team, and no further modifications to the lens models were performed. The predicted Fermat potential differences from both teams were combined with the observed kinematics data and estimated external convergence to predict the time delays between the three image pairs. A future study will infer the Hubble constant by comparing the predicted time delays with the observed ones from ongoing monitoring campaigns. We investigated the observed systematic differences between the model outputs from the two teams and identify that differences in the reconstructed PSF are the dominant source of systematic differences. The main results of this study are as follows:

• The final predicted time delays from LENSTRONOMY are: d, d, and d; and the ones from GLEE are: d, d, and d. These values assume a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3. The negative value of Δt_AB, for example, signifies that image B lags image A. This system is currently being monitored under the COSMOGRAIL programme to measure the time delays (Eigenbrod et al. 2005). Once the time delays are measured, the mutual agreement between the predicted values will result in a mutually consistent inference of the Hubble constant.

• The logarithmic slope, γ, and the external shear, γ_ext, of the power-law model deviate by 4σ (3.9%) between the LENSTRONOMY and GLEE models. This discrepancy is predominantly created by the difference in the reconstructed pixelized PSFs. When the same PSF is used by both modelling programs, then the resultant γ and γ_ext distributions agree within ∼0.6σ (∼0.5%), which is compatible with the ∼1% precision goal in the Hubble constant measurement from a large sample of ∼40 quad lenses (e.g. Shajib et al. 2018; Birrer & Treu 2021). The particular method of weighting the uncertainty map to account for the PSF uncertainty is non-dominant in this discrepancy.

• Our composite model posteriors are not generally in good agreement with the one from the power-law model, and the discrepancy is more prominent for LENSTRONOMY models due to the adoption of more stringent physical priors on the halo mass. This inconsistency points to the inadequacy of our composite model in describing the mass distribution for this particular lens system, WGD 2038−4008; WGD 2038−4008 is atypical compared to previously analysed TDCOSMO lenses in that the θ_E/θ_eff is relatively small and thus the imaging information only probes the inner region of the deflector galaxy where the combination of NFW and mass-follows-light profiles is not an adequate model. We stress that for all the seven systems previously analysed by the TDCOSMO collaboration, the power-law and composite models were in excellent agreement. For such discrepancies between the power-law and composite models, additional data, such as the stellar kinematics, should be used to select the better model. The two models predict significantly different velocity dispersion profiles and will therefore be easily separable by spatially resolved kinematics.

In the context of the recent debate around the Hubble constant, it is paramount to thoroughly investigate for potential systematic biases in the measurement methods (e.g. Freedman 2021; Riess et al. 2022). This study performs one such crucial systematic check for time-delay cosmography to investigate the robustness of lens modelling software programs. By keeping the modelling systematics under control, future large samples of lensed quasars and SNe will robustly measure the Hubble constant to ≲1% precision (e.g. Jee et al. 2016; Birrer & Treu 2021; Birrer et al. 2022).

Acknowledgments

We thank Christopher D. Fassnacht, Veronica Motta, and Abigail Lee for useful comments and discussions that improved this manuscript. Support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51492 awarded to A.J.S. by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. A.J.S., S.B., and T.T. were supported by the National Aeronautics and Space Administration (NASA) through the Space Telescope Science Institute (STScI) grant HST-GO-15320. A.J.S. was also supported by the Dissertation Year Fellowship from the University of California, Los Angeles (UCLA) graduate division. This research was supported by the US Department of Energy (DOE) Office of Science Distinguished Scientist Fellow Program. This work was supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. This work was supported by JSPS KAKENHI Grant Number JP20K14511. S.H.S. thanks the Max Planck Society for support through the Max Planck Research Group. This research is supported in part by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. T.T. acknowledges support by the Packard Foundation through a Packard Research fellowship, by the National Science Foundation through NSF grants AST-1836016, AST-1906976, and by the Gordon and Betty Moore Foundation through grant 8548. AA’s research is funded by Villum Experiment Grant Cosmic Beacons (project number 36225). This work was supported by the Swiss National Science Foundation (SNSF) and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 787886). This work used computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology Group. Additionally, this work was completed in part with resources provided by the University of Chicago’s Research Computing Center. A.J.S. and S.H.S. acknowledge the hospitality of the Aspen Center of Physics (ACP) and the Munich Institute for Astro- and Particle Physics (MIAPP) of the Excellence Cluster “Universe”, where part of this research was completed. A.C.P. is supported by National Science Foundation grant PHY-1607611. This research made use of GLEE (Suyu & Halkola 2010; Suyu 2012), LENSTRONOMY (Birrer et al. 2015; Birrer & Amara 2018; Birrer et al. 2021), FASTELL (Barkana 1999), NUMPY (Oliphant 2015), SCIPY (Jones et al. 2001), ASTROPY (Astropy Collaboration 2013, 2018), PHOTUTILS – an ASTROPY package for detection and photometry of astronomical sources (Bradley et al. 2020), JUPYTER (Kluyver et al. 2016), MATPLOTLIB (Hunter 2007), SEABORN (Waskom et al. 2014), SEXTRACTOR (Bertin & Arnouts 1996), EMCEE (Foreman-Mackey et al. 2013), COLOSSUS (Diemer 2018), GETDIST (https://github.com/cmbant/getdist), PYGALAXEV (Bruzual & Charlot 2003), and DYNESTY (Skilling 2004; Speagle 2020).

References

All Tables

All Figures

	Fig. 1. False-colour image of the lens systems WGD 2038−4008. This RGB image is created from the F160W (red), F814W (green), and F475X (blue) filters of the HST WFC3. We adjusted the relative amplitudes between the three filters to achieve a higher contrast for better visualization. The four lensed quasar images are marked as A, B, C, and D.
In the text

Fig. 2. Fiducial power-law model results for IR/F160W (top row), UVIS/F814W (middle row), and UVIS/F475X (bottom row) from GLEE. The maximum-likelihood model in the MCMC chain is shown. Shown are the observed image (first column), the reconstructed image predicted by the model (second column), the normalized residual within and interior to the arcmask region (defined as the difference between the data and model, normalized by the estimated uncertainty of each pixel; third column), and the reconstructed source (right column). First column: the dotted cyan lines indicate the arcmask (donut-shaped) region used for fitting the extended source, the dotted orange lines indicate the AGN mask region where the power-law weighting is applied, and the region outside the dotted cyan arcmask is used to further constrain the foreground lens light and (partly) the AGN light (but not the AGN host galaxy light since its corresponding lensed arcs are below the noise level in this outer region). The colour bars show the scale in the respective panels. The results shown here are for the fiducial power-law model, but the results for the other systematics tests (Sect. 5.3) are qualitatively similar.

In the text

	Fig. 3. Same as Fig. 2, but for the fiducial composite model from GLEE.
In the text

	Fig. 4. Marginalized parameter distributions from our power-law lens model results from GLEE. We show the combined results from our systematics tests (dashed red contours) with each model weighted equally, as well as the BIC-weighted model results (shaded red contours). The contours represent the 68.3% and 95.4% quantiles.
In the text

	Fig. 5. Marginalized parameter distributions from our composite lens model results from GLEE. We show the BIC-weighted model (shaded blue contours) and the combined results from our systematics tests (dashed blue contours). The contours represent the 68.3% and 95.4% quantiles.
In the text

Fig. 6. Radial mass profiles of the central deflector constrained by the GLEE models. Top: circularly averaged convergence ⟨κ(< R)⟩ as a function of radius for the GLEE power-law model (red) and composite model (blue). The shaded regions represent the 1σ credible regions. The stellar (green) and dark matter (black) components of the composite model are plotted separately. The vertical dashed black lines mark the pixel size in the F160W band and the best fit Einstein radius. Bottom: ratio of average convergence of the composite model to that of the power-law model as a function of radius.

In the text

	Fig. 7. Model-predicted distributions of Fermat potential differences (blinded) for each of the GLEE models tested, with power-law models (top) and composite models (bottom).
In the text

Fig. 8. Most likely LENSTRONOMY lens model and reconstructed image of WGD 2038−4008 using the power-law model. The top row shows, from left to right, the observed RGB image, the reconstructed RGB image, the convergence profile, and the magnification model. The next three rows show, from left to right, the observed image, the reconstructed image, the residual, and the reconstructed source for each of the HST filters. The three filters are F160W (second row), F814W (third row), and F475X (fourth row). All the scale bars in each panel correspond to 1″. The star symbol in the reconstructed source panels marks the position of the quasar host galaxy’s centroid.

In the text

Fig. 9. Most likely LENSTRONOMY lens model and reconstructed image of WGD 2038−4008 using the composite model. The top row shows, from left to right, the observed RGB image, reconstructed RGB image, the convergence profile, and the magnification model. The next three rows show, from left to right, the observed image, the reconstructed image, the residual, and the reconstructed source for each of the HST filters. The three filters are F160W (second row), F814W (third row), and F475X (fourth row). All the scale bars in each panel correspond to 1″. The star symbol in the reconstructed source panels marks the position of the quasar host galaxy’s centroid. In the magnification model, a central image is predicted due to a central core in the triple Chameleon light profile. However, this central image is highly de-magnified, with magnification 0.019 ± 0.02, and thus its presence cannot be ascertained in our imaging data.

In the text

Fig. 10. Comparison between the Fermat potential difference posteriors from different LENSTRONOMY model settings. The posterior for a particular setting is obtained by averaging over models that differ in the other model settings, but within a particular model family using the procedure described in Sect. 6.3. Top two rows: correspond to the power-law mass model families, and bottom two rows: correspond to the composite mass model families. The shaded regions illustrate the 68% credible regions. The posteriors are blinded by , where the ‘ref’ subscript refers to one of the compared models. The potential differences are consistent within 1σ between our adopted choices of model settings.

In the text

Fig. 11. Comparison of the LENSTRONOMY lens model parameters and Fermat potential difference between the power-law (red) and composite (blue) mass models. The posteriors are obtained by averaging over all the model settings following the procedure described in Sect. 6.3. The parameter γ for the composite model is computed from circularly averaging the quantity 2 − [d log α(r)/d log r]r = θE. Some parameters are blinded as for p ∈ {γ,  Δϕ}, where the subscript ‘pl’ refers to the power-law model posteriors. The composite model-predicted γ is approximately 8% smaller than that from the power law. Consequently, the Fermat potential differences are smaller by approximately 16% for the composite model.

In the text

Fig. 12. Distribution of M200 and c200 parameters for the NFW halos in LENSTRONOMY composite model (blue shaded region). This distribution is averaged over all the model settings within the composite model family following the procedure described in Sect. 6.3. The 2 contours correspond to the 68% and 95% credible regions, respectively. The black solid line traces the theoretical prediction of the M200 − c200 relation at zd = 0.230 from Diemer & Joyce (2019) with the grey shaded region corresponding to the 68% confidence interval. We adopt this M200 − c relation as a prior in our analysis in addition to a M200 prior based on Sonnenfeld et al. (2018).

In the text

Fig. 13. Radial mass profiles of the central deflector constrained by the LENSTRONOMY models. Top: circularly averaged convergence ⟨κ(< R)⟩ as a function of radius from LENSTRONOMY power-law (red) and composite (blue) models. The stellar (green) and dark matter (grey) distributions in the composite model are also individually illustrated. The shaded regions encompass the 16th and 84th percentile of the sampled profiles for the corresponding case. The grey points illustrate the surface brightness profile fitted with isophotes as described in Sect. 6.1.1. The amplitude of the isophotal fit is normalized to match with the triple Chameleon profile (green shaded region) at θE for the purpose of this illustration. The vertical black dashed lines mark the pixel size in the F160W band and the best fit Einstein radius. Bottom: ratio of the circularly averaged convergence profiles between the composite and the power-law models. At the Einstein radius the convergence slope deviates by 16−21% between the two model setups.

In the text

Fig. 14. Line-of-sight velocity dispersion profile (circularized) corresponding to the lens model posteriors in the power-law (red) and composite (blue) models from LENSTRONOMY. Isotropic stellar orbits are assumed in computing these velocity dispersion profiles. The solid lines correspond to the median and the shaded regions encompass the 16th and 84th percentiles. The composite profile predicts a decrease in the velocity dispersion towards the centre due to the flattened core – which is not observed in local massive ellipticals (e.g. Cappellari 2016; Ene et al. 2019) – thus pointing to the atypicality of the posterior mass profile in the composite model.

In the text

Fig. 15. Fitting the lens galaxy light profile from a much larger cutout than that used for lens modelling. Left: large cutout around the deflector galaxy in the F814W band, with the lensed arcs and quasar images subtracted using the best-fit of a power-law lens model. The galaxy extends much further beyond the Einstein radius with θeff/θE ≈ 1.7. Middle: model for the surface brightness profile of the deflector constructed using elliptical isophotes. Right: model residual showing that the model has captured the overall light distribution despite numerous overlapping objects. The model from middle panel was further approximated with an elliptical multi-Gaussian expansion (MGE; Cappellari 2002) to allow deprojection along the LOS for kinematic modelling.

In the text

Fig. 16. Comparison of the Fermat potential posteriors between two different model settings for the LENSTRONOMY composite model. The blue solid distributions correspond to our primary model settings. In the test model setup (black dashed distributions), we fix the triple Chameleon profile for the stellar mass and light distributions to best fit values from a separately optimized model with the power-law mass profile. We also mask the central region of the deflector galaxy in the test model setup. Additionally, the prior on the halo mass profile is different. Whereas we adopt a prior on the NFW scale radius rs in the primary setup, we adopt a combination of priors on M200 and M⋆ in the test setup. Despite multiple differences in the model settings, the Fermat potential differences are consistent within 1σ with each other.

In the text

Fig. 17. Checking for systematics in the kinematic computation. Our primary settings for kinematic computation adopts the F814W light profile with 3% uncertainty in θeff (red). The orange contours are for the case with 20% uncertainty in θeff, and the black contours are for the case with F160W light adopted in the kinematic computation. The vertical blue line marks the measured LOS velocity dispersion, which has an uncertainty of 19 km s−1. The difference in the computed velocity dispersion between these cases is negligible, and it impacts the Fermat potential differences by ≲0.1%.

In the text

Fig. 18. External convergence distribution from Buckley-Geer et al. (2020) with additional weighting applied based on the predicted external shears for the power-law (red lines) and composite (blue lines) models from LENSTRONOMY (solid lines) and GLEE (dashed lines). Each illustrated GLEE distribution is a BIC-weighted combination of multiple κext distributions corresponding to the external shear constraint from individual lens model setups. In contrast, each illustrated LENSTRONOMY distribution is a single κext distribution corresponding to the combined (through BIC weighting) external shear value from all the model setups within a model family (i.e. power law or composite). The κext distributions used by one team were not revealed to the other team before the un-blinding to maintain independence.

In the text

Fig. 19. Fermat potential differences Δϕ = Δϕmodel(1 − κext) from LENSTRONOMY for the power-law (red) and composite (blue) model with the kinematics information folded into these individual model posteriors. The folding in of the kinematics information is performed through importance sampling from the posterior weighted by the kinematics likelihood (Eq. (41)). Next, the two model posteriors are joined together with equal weights and then the kinematics information is folded in. The combined posterior (purple) mostly resemble the power-law posterior, as the composite posterior is heavily down-weighted by the kinematics likelihood.

In the text

Fig. 20. MST-adjusted Fermat potential differences Δϕ = Δϕmodelλint(1 − κext) from LENSTRONOMY for the power-law (red) and composite (blue) model families. The internal MST parameter λint is estimated by combining the lens models with the measured velocity dispersion and the external convergence estimate. After adjusting for the MST, the power-law and composite model predictions for the Fermat potential differences become consistent with each other within 1σ.

In the text

	Fig. 21. Comparison between the two modelling teams for the predicted time delays (un-blinded) as a function of H0 for the three image pairs, in flat ΛCDM cosmology with Ωm = 0.3. Each posterior is the final combined posterior from the two mass-model setups: power-law and composite, including external convergence and stellar kinematics, and assuming λint = 1.
In the text

Fig. 22. Comparison of predicted time delays between the LENSTRONOMY (green) and GLEE (purple) teams. The three rows show the comparison for combined (top), power-law-only (middle), and composite-only (bottom) cases. The combined posteriors are consistent within < 1.2σ, the power-law model’s posteriors are consistent within < 0.6σ, and the composite model’s posteriors are consistent within < 2.1σ. In percentage, the maximum deviation for the power-law model is between AB image pair by 3.4%, and for the composite model is between AC image pair by 15%. The GLEE team combined the power-law and composite model posteriors with equal weighting after applying the kinematics weighting within each model family, whereas the LENSTRONOMY team combined the two model families with kinematics weighting, leading to the final combined posterior being dominated by the power-law model.

In the text

	Fig. 23. Comparison of lens model parameter differences for the power-law model between the LENSTRONOMY (green) and the GLEE (purple) teams.
In the text

	Fig. 24. Comparison of Fermat potential differences for the power-law model between the LENSTRONOMY (green) and the GLEE (purple) teams, without including the external convergence (top row), and with including the external convergence (bottom row).
In the text

Fig. 25. Comparison of the reconstructed pixelized PSFs by the LENSTRONOMY (first column) and the GLEE (second column) modelling procedures. Third column: illustrates the difference between the LENSTRONOMY and GLEE PSFs. The three rows correspond to the F160W, F814W, and F475X filters from top to bottom. The F160W PSF is supersampled with a supersampling factor of 3. We note that the illustrated supersampled LENSTRONOMY PSF was not used in the pre-un-blinding models and was reconstructed after the un-blinding to perform further tests. The original reconstructed PSF in LENSTRONOMY in F160W was not supersampled. The GLEE PSF FWHMs are smaller by ∼2−7% than the LENSTRONOMY ones.

In the text

Fig. 26. Deviations in the logarithmic slope γ of the power-law model with different PSF settings. In all panels, the dashed distributions show the fiducial constraints – LENSTRONOMY in green and GLEE in purple. The LENSTRONOMY fiducial constraint is from the highest BIC value power-law model (Table 2), and the GLEE fiducial constraint is from the ‘power-law fiducial’ setup (Table 1). First panel: when only the reconstructed PSFs are interchanged to optimize the models, the constraint from one software program moves towards the other’s fiducial constraint to be consistent within 1.2 − 1.6σ. We note that GLEE model artificially creates supersampled version of the LENSTRONOMY PSF through interpolation to perform the tests here. Second panel: when both the PSFs and weighted uncertainty maps are interchanged, the resultant γ becomes slightly more consistent between the software programs within 1.0 − 1.5σ. Therefore, we conclude that the differences in the PSF itself significantly contributes to the discrepancy in the power-law γ constraint and not the particular method of weighting the uncertainty map to account for PSF uncertainty. Third panel: when models from both of the software programs are optimized only with UVIS data and interchanged PSFs, the GLEE constraint does not shift significantly; however, the LENSTRONOMY constraint shifts significantly towards the GLEE constraints. This result indicates that the difference between the GLEE and LENSTRONOMY fiducial models are largely created by the difference in the IR PSF. Fourth panel: when both software programs use a supersampled PSF with a supersampling factor of 3 reconstructed by LENSTRONOMY, the resultant γ constraint agrees very well, within 0.6σ.

In the text