IX. Systematic comparison between lens modelling software programs Time-delay prediction for WGD 2038-4008

The importance of alternative methods for measuring the Hubble constant, such as time-delay cosmography, is highlighted by the recent Hubble tension. It is paramount to thoroughly investigate and rule out systematic biases in all measurement methods before we can accept new physics as the source of this tension. In this study, we perform a check for systematic biases in the lens modelling procedure of time-delay cosmography by comparing independent and blind time-delay predictions of the system WGD2038 − 4008 from two teams using two di ﬀ erent software programs: glee and lenstronomy . The predicted time delays from the two teams incorporate the stellar kinematics of the deﬂector and the external convergence from line-of-sight structures. The un-blinded time-delay predictions from the two teams agree within 1 . 2 σ , implying that once the time delay is measured the inferred Hubble constant will also be mutually consistent. However, there is a ∼ 4 σ discrepancy between the power-law model slope and external shear, which is a signiﬁcant discrepancy at the level of lens models before the stellar kinematics and the external convergence are incorporated. We identify the di ﬀ erence in the reconstructed point spread function (PSF) to be the source of this discrepancy. When the same reconstructed PSF was used by both teams, we achieved excellent agreement, within ∼ 0.6 σ , indicating that potential systematics stemming from source reconstruction algorithms and investigator choices are well under control. We recommend that future studies supersample the PSF as needed and marginalize over multiple algorithms or realizations for the PSF reconstruction to mitigate the systematics associated with the PSF. A future study will measure the time delays of the system WGD2038 − 4008 and infer the Hubble constant based on our mass models.


Introduction
The Hubble constant, H 0 , 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 0 , at the current epoch is extrapolated using Λ cold dark matter (ΛCDM) cosmology. This early-Universe probe resulted in constraints of H 0 = 67.4 ± 0.5 km s −1 Mpc −1 (Planck Collaboration 2020) and H 0 = 67.6 ± 1.1 km s −1 Mpc −1 (Aiola et al. 2020). In the local Universe, H 0 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 0 for the Equation of State of dark energy (SH0ES) team used Cepheids and parallax distances to calibrate the cosmic distance ladder, measuring H 0 = 73.04 ± 1.04 km s −1 Mpc −1 (Riess et al. 2022), which is in 5σ tension with the Planck measurement. The Carnegie-Chicago Hubble Project NHFP Einstein Fellow Packard Fellow NHFP Einstein Fellow used the tip of the red giant branch to calibrate the distance ladder, measuring H 0 = 69.6 ± 1.9 km s −1 Mpc −1 (Freedman et al. 2019(Freedman et al. , 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 0 = 73.9 ± 3.0 km s −1 Mpc −1 (Pesce et al. 2020), the Tully-Fisher method calibrated with Cepheids measured H 0 = 75.1 ± 0.2 ± 3.0 km s −1 Mpc −1 (Kourkchi et al. 2020), and the surface brightness fluctuation method measured H 0 = 73.7±0.7±2.4 km s −1 Mpc −1 (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 'timedelay distance' . The time-delay distance is inversely proportional to H 0 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 as-pects have improved by a large margin over the past decade (for a review with a historical perspective, see Treu & Marshall 2016). Inferring H 0 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 0 Lenses In the COSMO-GRAIL's Wellspring (H0LiCOW) and the Strong-lensing High Angular Resolution Programme (SHARP) collaborations measured H 0 = 73.3 +1.7 −1.8 km s −1 Mpc −1 from a sample of six strongly lensed quasar systems Bonvin et al. 2017;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 H 0 = 74.2 +2.7 −3.0 km s −1 Mpc −1 . 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 COS-MOgraphy (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 0 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 H 0 = 74.5 +5.6 −6.1 km s −1 Mpc −1 , that is, relaxing the power-law assumption leads to an increase in H 0 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 1 (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 0 shifted to H 0 = 67.4 +4.1 −3.2 km s −1 Mpc −1 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 0 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 0 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 1 Advanced Camera for Surveys 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 2 and lenstronomy 3 . 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 collaborationfive 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 RXJ1131−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;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 0 measurement, we leave the H 0 inference from our models to be done in the future. However, we predict the time delays for this system as a function of H 0 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 0 measurement from this system. The baseline models for comparison have two different lens model setups: (i) a power-law mass model and (ii) a twocomponent 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 TD-COSMO 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 0 cannot be inferred. Measuring H 0 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. Sects. 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.

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.

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. Shajib et al. 2020).
The delay ∆t XY between arrival times of photons corresponding to images labelled as X and Y is given by 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 ψ = 2κ. The convergence is the surface mass density scaled by the critical density as κ ≡ Σ/Σ crit with The Fermat potential φ is defined by combining the geometric delay term with the deflection potential as The so-called time-delay distance is defined as Each distance term contains a factor of H −1 0 , which cancel out such that D ∆t ∝ H −1 0 . Eq. 1 can be written in short form as

Mass-sheet degeneracy
The imaging observables of the lensing phenomenon -the image positions and the flux ratios -remain invariant under the transformation 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 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 RXJ1131−1231, HE 0435−1223, and ES J0408−5354 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 lim θ→∞ κ true (θ) = κ ext .
Therefore, κ ext can be interpreted as lensing mass in the 3D space far from or 'external' to the central deflector. Let κ model be the model convergence that can reproduce the imaging observables. However, due to the MSD, κ model is not a unique solution and we cannot ascertain that κ true = κ model . If we impose the condition Article number, page 3 of 34 A&A proofs: manuscript no. 2038_model_final lim θ→∞ κ model = 0, then κ model is a mass-sheet transform of κ true with the rescaling factor λ = 1/(1 − κ ext ) as If the external convergence κ ext can be independently estimated by studying the lens environment, then the true lensing convergence κ true can be recovered from κ model through the corresponding inverse MST. However, the lens model κ model that we actually constrain can be an internal MST of κ model as Interestingly, both κ model and κ model 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, 9, and 10 to write the relation between the true mass distribution κ true and the modelled mass distribution κ model as 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 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 Yıldırım et al. 2021). The time delay rescales with the MST as ∆t → ∆t = λ∆t.
As a result, we need to correct the time delays ∆t model predicted by the model κ model as In the next section, we describe our framework for the kinematics analysis.

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 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 The observable quantity is the luminosity-weighted LOS velocity dispersion, which we can obtain by solving the Jeans equation as where G is the gravitational constant, I(R) is the surface brightness, and M(r) is the 3D enclosed mass within radius r (Eqs. A15-A16 of Mamon & Łokas 2005). The function K β (r/R) depends on the parameterization of β ani (r). We adopt the Osipkov-Merritt parameterization given by where r ani is a scaling radius (Osipkov 1979;Merritt 1985a,b). For this parameterization, the form of K β (r/R) is given by with u ani ≡ r ani /R (Mamon & Łokas 2005). The observed aperture-averaged velocity dispersion is where ap denotes integration over the aperture and * S denotes convolution with the seeing. Thus, the lens-model-predicted LOS velocity dispersion can be written in the form 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 σ 2 ap, true = (1 − κ ext )λ int σ 2 ap, model .
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).

Bayesian inference
We denote the set of all the observables as 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 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 Sthat 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 Then, we can first perform the following sub-integral within the right-hand side of Eq. 23: 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 where k is the number of free parameters, N data is the number of data points, andL is the maximum of the likelihood function L.
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)  Specific implementations of the Bayesian inference framework presented in this section through sampling by each team are described in Sects. 5 and 6.

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 . The deflector redshift is z d = 0.230±0.002 and the source redshift is z s = 0.777±0.001 . In this section we describe the imaging data and spectroscopic measurements used in our analysis.

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 MULTI-ACCUM mode. The total exposure times are 2196.9 s, 1428.0 s, and 1158.0 s, respectively, in the three filters. We show a falsecolour red-green-blue (RGB) image of the system created from the HST imaging in Fig. 1.
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).

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 −1 from a 0 . 75×1 rectangular aperture, which is in agreement with a more recent measurement from the XShooter 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 0 . 9, and the exponent parameter of the Moffat PSF is β = 1.74.

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;Rusu et al. 2020). 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. (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.

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. More-over, 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.

Mass profiles
The two baseline lens model families we adopt are the powerlaw mass profile and the composite mass profile.
Article number, page 6 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.: Systematic comparison between lens modelling software programs 4.1.1. Power-law mass profile We adopted the power-law elliptical mass distribution (PEMD; Barkana 1998) defined as where γ is the logarithmic slope, θ E is the Einstein radius, and q m is the axis ratio. The coordinates (θ 1 , θ 2 ) 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.

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 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: where a 0 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.

Sérsic profile
The Sérsic profile is defined as 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.

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 κ 0 with the flux amplitude I 0 .

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 0 . 004 to account for offsets due to substructure in the lens or LOS. This uncertainty is small enough to satisfy astrometric requirements for cosmography . 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 5 . 6 × 5 . 6, 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 construc-Article number, page 7 of 34 A&A proofs: manuscript no. 2038_model_final tion 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 A × σ b img,i . 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.

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 (θ 1 , θ 2 ) 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. Fig. 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.

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 powerlaw 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.
A Gaussian prior for the M/L of the baryonic component is employed, using the stellar mass constraint from Agnello et al. (2018) of log 10 (M /M ) = 11.40 +0.01 −0.08 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 r s = 22 . 6 ± 3 . 1 based on the results of Gavazzi et al. (2007) for a sample of lenses in the SLACS survey . These lenses span a redshift and velocity dispersion range that includes WGD 2038−4008, with a mean virial mass of M vir = 1.4 +0.6 −0.5 × 10 13 h −1 M . 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. Fig. 3 shows the data and the lens model results in all three bands for our fiducial composite model, as well as the source reconstructions.

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).

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 Article number, page 8 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.: Systematic comparison between lens modelling software programs 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). In the 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. 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. 22 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 0 (via Eqs. 5 and 14).

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).L (the maximum likelihood of the model from the MCMC sampling) is the product of the AGN position likelihood, the pixellated im-age 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, σ 2 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 σ 2 BIC . We find σ 2 BIC ∼ 36 for the power-law models and σ 2 BIC ∼ 34 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  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 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 σ 2 BIC . 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.

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. Fig. 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.
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 ∼ 0 . 003, and the model is able to fit the quasar positions to an rms of ∼ 0 . 005.
The composite model fits the quasar positions to an rms of ∼ 0 . 01, 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 (0 . 04), which corresponds to ∼ 2 − 3%. At the Einstein radius, the composite model slope closely matches  Fig. 7. The unblinded 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%.

Lenstronomy modelling
In this section we describe the lenstronomy model setups and modelling results. The software package lenstronomy (Birrer & Amara 2018;) is a publicly available lens modelling software. 4 . 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.

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 com- putation in Sect. 6.1.3. We summarize the set of all the lens models combining these different settings in Sect. 6.1.4.

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 22 . 6 and standard deviation 3 . 1 for the NFW scale radius r s based on the measurements of Gavazzi et al. (2007) for a sample of SLACS survey lens systems ). 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; Article number, page 12 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.: Systematic comparison between lens modelling software programs  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 200 -c relation from Diemer & Joyce (2019) with an intrinsic scatter of 0.11 dex.

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. 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. 5 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 5 using the 'PROFILE_SHEAR' solver of lenstronomy shapelet components. We adopt three choices for {n IR max , n UVIS max }: {7, 11}, {8, 12}, {9, 13}.

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 {r IR L , r UVIS L }: {2 . 2, 3 . 6}, {2 . 3, 3 . 7} with r L . 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). 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.

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  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 10000 steps. We check for convergence of the chain by manually inspecting that the median and standard deviations of the parameters within the Article number, page 14 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.:   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.

Bayesian model averaging
We have 24 models that make up our set of models {S , M}, with each lens model family from M ≡ {power law, 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 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 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 σ model ∆BIC 2 of ∆BIC between 'neighbouring' models that differ with each other by one step in only one setting . We furthermore accounted for the numeric uncertainty in estimation of ∆BIC by taking the variance σ numeric ∆BIC 2 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 . Thus, our total variance in ∆BIC is We compute that σ model ∆BIC = 304 and σ numeric ∆BIC = 69. To implement the ∆S n weighting through sampling, we first follow  to obtain the absolute weight W n,abs of the n th model by convolving the ∆BIC with the evidence ratio function f (x) as where the evidence ratio function f (x) is defined using the BIC difference as Then, we obtain the relative weight W n,rel by normalizing the absolute weights by the maximum absolute weight as W n,rel = W n,abs max W n,abs .  In Sect. 6.4, we compare the two mass model families after performing the above model averaging procedure within each model family.

Lensing-only constraints on the Fermat potential difference
We constrain the image positions with uncertainty 0 . 002 for the power-law model, and with uncertainty 0 . 004 for the composite model. Given the longest predicted time-delay for this system, these precisions are well below the astrometric requirement of ∼0 . 02 uncertainty so that the astrometric uncertainty is subdominant to achieve ≤5% uncertainty in H 0 from this single system ).
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σ.
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 ξ rad ∝ θ E α E /(1 − κ E ) 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.
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.  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. centration posterior is consistent with the concentration prior within 1σ. Our combined posterior from the composite model setup provides the total halo mass log 10 (M 200 /M ) = 13.04 +0.14 −0.13 and the total stellar mass is log 10 (M /M ) = 11.87 +0.01 −0.03 . The total stellar mass is obtained by doubling the enclosed mass within the half-light radius of 3 . 2 corresponding to the F160W band. The projected dark matter fraction within the Einstein radius is 0.22 +0.06 −0.02 . The total baryon-to-dark-matter fraction is 0.07 +0.03 −0.02 , which is consistent with the upper limit set by the cosmic baryonic fraction 0.19 (Planck Collaboration 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 ( 0 . 2) 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 0 = 70 km s −1 Mpc −1 and Ω m = 0.3. This stellar core radius is consistent with the core radii measured in local elliptical galaxies

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 p(log 10 M 200 | log 10 M Chab ) where M Chab 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 M Chab 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 6 to obtain M Chab , 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 log 10 M Chab = 11.57 +0.16 −0.13 . Thus, from Eq. 38, our Gaussian prior on the halo mass log 10 M 200 has mean 13.5 and standard deviation 0.3. We additionally adopt a prior for the total stellar mass log 10 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 2×0.13 2 , log 10 M / < 11.51, A, 11.51 ≤ log 10 M / ≤ 11.88, A exp − (log 10 M / −11.88) 2 2×0.16 2 , log 10 M / > 11.88, where M / ≡ M /M , A is the amplitude that normalizes the probability distribution to have p(log 10 M / ) d(log 10 M / ) = 1. 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 0 in the SPSbased 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 log 10 (M /M ) = 11.84 +0.02 −0.01 , the halo mass is log 10 (M 200 /M ) = 13.3 +0.1 −0.3 , the total baryonic fraction is 0.04 +0.03 −0.01 , and the dark matter fraction within the Einstein radius is 0.28 +0.02 −0.03 . 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 powerlaw models.
The mass difference of 3.8×10 10 M within 0 . 2 (assuming flat ΛCDM cosmology with H 0 = 70 km s −1 Mpc −1 ) 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.

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 θ eff = 2 . 4. The aperture size in this numeric computation can only grow by a size of a pixel, which is 0 . 08. 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. (2016Birrer et al. ( , 2020 find that the uniform prior for log a ani is a less informative choice than the uniform prior on a ani . 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 combine the stellar kinematics and external convergence information with the lens model posterior in two different ways: used by the two teams was revealed to each other only after the unblinding. While this creates an inconsistency between the two teams, Rusu et al. (2020) show that for large shear values, the κ ext distribution is dominated by the shear constraint, and therefore the imposition of the 45 aperture or the lack thereof is expected to have a negligible impact.
Article number, page 20 of 34 3RZHUODZ &RPSRVLWH 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 p blinded ≡ p/p pl − 1 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. 6.5.1. The case with λ int = 1 Assuming λ int = 1, the model-predicted velocity dispersion can be written as 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 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 L 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.
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   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.
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 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 0 , the jerk parameter j 0 , the snap parameter s 0 , 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 0 based on its highest model evidence. We use the relation D A = D L /(1+z) 2 to convert the luminosity distance to angular diameter distance and transform the 2D posterior distribution of (q 0 , j 0 ) 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 λ pl int = 1.02 ± 0.15 for the power-law model and λ comp int = 1.79 ± 0.53 for the composite model. This large median value of λ comp int 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).
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 0 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 λ int ∝ σ 2 ap , the uncertainty on the λ int is thus expected to be twice the uncertainty of the velocity dispersion. Our obtained uncertainties of λ pl int and λ comp int are consistent with this expectation. The uncertainties on the predicted time delays from the final combined posterior with free λ int is 19-22%.

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 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.

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 Article number, page 23 of 34 A&A proofs: manuscript no. 2038_model_final requirements in Sect. 7.4. Finally, we discuss our findings in Sect. 7.5.

Predicted time delays
We illustrate the final predicted time delay from both teams as a function of H 0 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 compositeonly 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 powerlaw-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σ).

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 0 . 004 at maximum. This astrometric precision satisfies the requirement for cosmographic measurements .
In Fig. 23 we compare the lens model parameter posteriors and the predicted Fermat potential posteriors from the powerlaw 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 Article number, page 24 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.: Systematic comparison between lens modelling software programs  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. a 2 4 + b 2 4 0.005), which implies that allowing boxy-ness or discy-ness in the description of the mass profile is not required (for definitions of a 4 and b 4 , and their impact on H 0 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).
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 effectthat 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.

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 (0 . 08) as the drizzled image. The PSFs for the UVIS filters have the same  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%.
pixel resolution as the drizzled image (0 . 04) for both glee and lenstronomy.
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.
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.

Requirements for computational resources
The entire lenstronomy modelling procedure required ∼ 4 × 10 5 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 O(10 5 ) CPU hours only to run the final models and MCMC chains to produce the final un-blinded result. However, Article number, page 26 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.: 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. the total usage of CPU hours can be O(10 6 ) CPU hours including modelling trials, robustness tests, and reruns of chains to recover lost progress due to numerical issues.

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 de-lays 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 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σ. 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 illus-trates 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 Article number, page 28 of 34 Shajib, Wong, Birrer, Suyu, Treu et al.: Systematic comparison between lens modelling software programs 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 0 . 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 re-construction algorithms to avoid any potential bias. Furthermore, a supersampled PSF in the IR band is also recommended, as the drizzled pixel scale of 0 . 08 in the IR does not optimally sample the PSF.

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: 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;. The particular method of weighting the uncertainty map to account for the PSF uncertainty is nondominant 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 et al. 2022). lenstronomy Glee Comparison of the reconstructed pixelized PSFs by the lenstronomy (first column) and the glee (second column) modelling procedures. The 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. γ Density lenstronomy supersampled IR PSF 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 constraintslenstronomy in green and glee in purple. The lenstronomy fiducial constraint is from the highest BIC value powerlaw model (Table 2), and the glee fiducial constraint is from the 'powerlaw 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σ.
Article number, page 33 of 34