© The Authors 2023

1 Introduction

Gravitational lensing describes the effect of the gravitational potential of a massive object on the light of a background source. In the case of strong lensing (SL), the light is deflected, such that multiple images of the source are observed.

Because SL is sensitive to both luminous and dark matter (DM), it is a powerful tool for probing the distribution of the total galaxy mass and thus for giving insight into the DM distribution of galaxies (e.g., Dye & Warren 2004; Barnabè et al. 2012; Sonnenfeld et al. 2015; Schuldt et al. 2019; Shajib et al. 2021). Mass clumps in the lensing galaxy affect the image magnification, and substructures can be detected by analyzing flux-ratio anomalies of point-like sources (e.g., Dalal & Kochanek 2001; Moustakas & Metcalf 2002; Nierenberg et al. 2014, 2017, 2020; Hsueh et al. 2020; Gilman et al. 2020) or distortions in the arcs formed by the lensing of spatially extended background galaxies (e.g., Koopmans 2005; Vegetti et al. 2010, 2012). Clumps with a mass lower than 10⁸ solar masses can be detected in this way. We can also analyze the structural parameters of early-type lens galaxies, such as the average slope of the total mass density profile or the DM fraction, with models (e.g., Gavazzi et al. 2007; Auger et al. 2010; Lagattuta et al. 2010; Barnabè et al. 2011; Shu et al. 2015).

In addition, the magnifying effect of SL can be used to observe the earliest galaxies and quasars in the Universe (e.g., Kneib et al. 2004; Bradley et al. 2008; Coe et al. 2013; Salmon et al. 2020). Modeling these high-redshift objects enables us to study their rotation curves and masses, for example (e.g., Jones et al. 2010; Chirivì et al. 2020; Rizzo et al. 2020). This helps us to understand galaxy evolution and structure formation in the early Universe. Despite major progress in the last years, this research field has many open questions.

Another important application of SL is the measurement of cosmological parameters such as the Hubble constant H₀ (e.g., Suyu et al. 2017; Wong et al. 2019; Shajib et al. 2019a; Chen et al. 2020; Millon et al. 2020c; Birrer et al. 2020). This parameter, which gives the current expansion rate of the Universe, is of great interest because early-universe measurements from the cosmic microwave background (CMB) give a value of 67.36 ± 0.54 km s⁻¹ Mpc⁻¹ (Planck Collaboration VI 2018), whereas local measurements from the SH0ES project (Riess et al. 2021) based on the cosmic distance ladder using Cepheids and type Ia supernovae (SNe) obtain a higher value of 73.04 ± 1.04 km s⁻¹ Mpc⁻¹. This 5σ disagreement is commonly known as the Hubble tension and might be an indication of physics beyond a flat ACDM. The flat ACDM currently is the standard cosmological model. SL can also be used to probe alternative cosmological models (e.g., Jullo et al. 2010; Oguri et al. 2012; Cao et al. 2015; Krishnan et al. 2021). To resolve the Hubble tension, independent methods are necessary, including megamasers (e.g., Pesce et al. 2020), standard sirens (e.g., Abbott et al. 2017), surface brightness fluctuations (e.g., Khetan etal. 2020), expanding photospheres of supernovae (e.g., Schmidt et al. 1994), and gravitational lensing.

To measure H₀ with SL, a strongly lensed variable background source such as a quasar or an SN is required. Because of the different path length and the gravitational potential differences on the way of the photons, the characteristic brightness fluctuations reach the observer at different times. We can measure this time delay in the light curves by monitoring the multiple images and comparing these flux variations with each other (e.g., Millon et al. 2020a,b).

In addition to an accurate measurement of the time delays, the lens potential is required to determine H₀, which is obtained by modeling the observed lensing system. Major progress has been made in this field in recent decades, in part because of the development of advanced modeling software such as GLEE (Suyu & Halkola 2010; Suyu et al. 2012) and Lenstronomy (Birrer et al. 2015; Birrer & Amara 2018). Modeling software like these allow us to constrain the model parameters of the lens and the source, such as their position and radial mass profile. A downside of this approach is that we need at least several weeks to model a lens system with spatially extended sources because it is time consuming to model the parameter space, for instance, with Markov chain Monte Carlo (MCMC) methods. Moreover, the modeling process itself is interactive. In addition, the input files required by GLEE (e.g., point spread function (PSF) and error map) have to be obtained manually from the available data in some cases.

With the increasing number of wide-field imaging surveys, including the Hyper Suprime-Cam (HSC; Aihara et al. 2018), the Dark Energy Survey (DES; Dark Energy Survey Collaboration 2016), the upcoming Euclid (Laureijs et al. 2011; Scaramella et al. 2022), and the Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezic et al. 2019), the number of known strong-lensing systems is growing rapidly. Even though the number of detected lenses per square degree lies between 0.5 and 10 (Li et al. 2020) depending on the image depth, the number of lens candidates will increase strongly as a result of the large survey areas. For example, the LSST will cover a total of ~20 000 deg² in multiple bandpasses over its planned ten-year run time (Tyson 2002; Lochner et al. 2022). We expect new images of billions of galaxies, about 100 000 of which are strong-lensing systems (Collett 2015). To be able to keep up with the increasing sample size, we need techniques that automate the lens-modeling process. By automating the creation of the input files, the calculation of the starting values of lens parameters and image positions, and the optimization and sampling of the parameter space, we can save many hours of user interaction during the modeling process. Several projects in recent years have already successfully demonstrated the automation of lens modeling, especially when initial estimates of the lens mass distribution were to be obtained (e.g., Shajib et al. 2019b; Nightingale et al. 2021). Techniques are also being developed based on machine learning (e.g., Hezaveh et al. 2017; Perreault Levasseur et al. 2017; Park et al. 2021; Schuldt et al. 2021, 2022), although the modeling of strongly lensed quasars via neural networks has yet to be tested on real instead of mock data. Another method is the massive parallelization of graphics processing units (Gu et al. 2022) to speed up the lens modeling. All these approaches cannot yet replace the detailed manual lens-by-lens analysis that is necessary to produce models that are accurate enough for time-delay cosmography (e.g., Suyu et al. 2010, 2013; Birrer et al. 2019; Wong et al. 2019; Shajib et al. 2019a; Rusu et al. 2020), but rather serve as a first step toward this and provide important information for follow-up observations.

In this paper, we present a new automated procedure for modeling of strongly lensed quasars that utilizes the lens-modeling software GLEE and efficiently employs the user time. The modeling is uniform (i.e., has the same assumptions) for all systems. We developed a pipeline written in Python that works with minimum user input. Providing only high-resolution imaging data and marking important regions in the field, the user can model in a nearly fully automated way with GLEE. The automated modeling can be performed with single-band and multiband data. We focus on lensed quasars to model the lens mass distribution with an isothermal or power-law profile and the source brightness distribution on a grid of pixels. To test our code, we uniformly modeled a sample of nine strongly lensed quasars observed with the Hubble Space Telescope (HST). Furthermore, we predict the time delays and Fermat potential differences between the multiple images of the quasar. We emphasize that our approach differs from that of papers aimed at deriving cosmological parameters from individual lenses, because we prioritize speed over accuracy. Owing to the adopted shortcuts (e.g., we do not iteratively correct the PSF and usually model only the images in filters with visible arcs from the host galaxy) and the uneven quality of the data, we do not in general expect our results to reach cosmography grade. The output of this pipeline is aimed as a first step toward cosmography-grade models, and it is meant to assist in prioritizing follow-up monitoring and further modeling.

Schmidt et al. (2022, hereafter S22) have also developed an automated modeling pipeline using Lenstronomy and applied it to a sample of 30 lensed quasars. Our sample of 9 quasars is a subset of the systems in S22, allowing us a direct comparison of the results between the two modeling pipelines as a blind test because the lead modelers of the two teams obtained their lens models independently and only compared the results after the models were completed. We note that our approach and that of S22 are similar, but not identical. Important differences include the modeling of the PSF, which is known to be crucially important (Shajib et al. 2022), and the strategy adopted to cope with insufficient data quality (S22 imposed informative priors, but we do not). Therefore, we do not expect the two efforts to agree within the formal uncertainties, which are only a representation of the statistical error and not of the systematics, which are notoriously dominant. This should therefore be kept in mind when interpreting the comparison results. Nevertheless, as we show below, important lessons can be gleaned from it.

The outline of the paper is as follows: in Sect. 2 we describe the steps in modeling a strong-lensing system with our automation code. We briefly describe each system and the observational data in Sect. 3. The modeling results are presented in Sect. 4. The performance of the automated modeling compared to manual modeling is discussed in Sect. 5, together with a comparison with the models obtained with Lenstronomy of the same systems (S22). In Sect. 6 we summarize our results.

Throughout the paper, we adopt a flat ACDM cosmology with H₀ = 70km s⁻¹ Mpc−¹ and Ω_M = 1 − Ω_A = 0.3. Parameter estimates are given by the median of its one-dimensional marginalized posterior probability density function, and the quoted uncertainties show the 16th and 84th percentiles (corresponding to a 68% credible interval).

2 Automated modeling process

In this section, we describe the steps of modeling a strongly lensed quasar with GLEE in detail and how they are conducted in our automation code. We adopt different mass and light profiles, and use the Bayesian framework of MCMC sampling and simulated annealing (Kirkpatrick et al. 1983) provided by GLEE. The goal is to model the lens mass distribution first using the source and image positions, and then by reconstructing the surface brightness distribution of the lens and background source galaxies. Figure 1 shows the strongly lensed quasar J2100-4452 as an example. The figure shows the lens galaxy in the image center, four multiple images of the background quasar, and the arc, which is the lensed light of the quasar host galaxy. The code models the mass distribution of the lens galaxy and the light components of the system (i.e., lens, satellite galaxies of the lens, lensed quasar, and arc light). The modeled constituents of the lensing system are labeled in Fig. 1.

Because the modeling is typically first performed in a single filter and further bands may be included after the model has stabilized, we present in Sect. 2.1 an automation code dedicated for the single-band modeling. We then introduce the automation code for the multiband modeling in Sect. 2.2 and finally predict the time delays between the multiple images of strongly lensed quasars in Sect. 2.3.

2.1 Single-band modeling

In this section, we discuss the automated modeling of a lensing system with single-band data. We describe the semi-automated preparation of the required input files and each modeling step and their substeps in the following subsections. An overview of the basic modeling procedure is shown as a flowchart in Fig. 2 and summarized in Table 1. In particular, Table 1 shows the varying and fixed parameters in each modeling step.

2.1.1 Preparation

To model a strong-lensing system with our newly developed automated modeling codes, the user needs to provide information about the position and approximate light structure of the lens, the quasar images, and the satellites. In addition, several input files (see Table 2) are needed. Each input file is created in separate codes that are described in the following.

Science image cutout We start by creating a cutout of the lensing system from the preprocessed reduced data such that it includes the immediate environment (≲10") that is to be taken into account during the modeling. In the automated pipeline, the user creates a DS9 (Joye & Mandel 2003) region file containing the cutout area in the field and the positions of the lens, satellite galaxies, and the quasar images. Optionally, the background flux is subtracted by calculating the average flux of an empty patch in the field. The modeling code then creates the cutout, which is used for the subsequent modeling process, and calculates starting values for the lens mass and the light parameters and image positions from the region information.

Error map The error map accounts for both background noise and Poisson noise. The background noise σ_background is approximated as a constant that is set to be the standard deviation from a small patch in the science image without astrophysi-cal sources. The Poisson noise σ_poisson is computed from the weighted image (exposure time map) and the science image. For images in units of counts per second, we calculate

$${\sigma }_{\text{poisson,}i}^2=\left|\frac{d_i}{t_i\times G}\right|,$$(1)

with d_i the observed intensity¹, t_i the exposure time for pixel i, and G the gain of the charge-coupled device (CCD) for the specific filter. In the case of data units of counts, we have

$${\sigma }_{\text{poisson,}i}^2=\left|\frac{d_i}{G}\right|.$$(2)

We add σ_poisson in quadrature to the background noise level only if σ_poisson > 2 × σ_background to avoid overestimating the noise in regions without flux from astrophysical objects. We check this condition for each pixel i to obtain the final error map:

$${\sigma }_{\text{tot,}i}^2=\{\begin{array}{ll}{\sigma }_{\text{background,}i}^2+{\sigma }_{\text{poisson}}^2\hfill & \text{if}{\sigma }_{\text{poisson}}>2\times {\sigma }_{\text{backgruund}};\hfill \\ {\sigma }_{\text{background,}i}^2\hfill & \text{otherwise.}\hfill \end{array}$$(3)

Point spread function If the PSF is not directly provided with the science image, the PSF can be created with our PSF code. For this, the user marks several nonsaturated stars in the field. The field spans ~2.5 × 2.5 arc min in the case of the HST WFC3 data we used. The PSF code then automatically creates cutouts of these stars that are centered at the brightest pixel. Since the star is not necessarily located at exactly the center of the brightest pixel, the PSF code interpolates and shifts the cutout within fractions of a pixel to center it in the cutout. In particular, the profile is shifted such that the brightness of the pixels to the left and right of the brightest pixel is within a margin of 10%, and similarly for the pixels above and below the brightest pixel. This causes the central 3 × 3 region to be roughly symmetric. To obtain the PSF, the recentered stars are normalized and then stacked by calculating the median per pixel among the cutouts, and they are renormalized such that the sum of all pixel values is 1. In addition, the automation code requires two PSFs, one for modeling the light of the lensed quasar, and one for the extended source. The first PSF includes the diffraction spikes as that PSF image is used directly as the spatial profile for the quasar images. The second PSF is cropped at the first diffraction minimum of the Airy disk because it is used for the convolution of the parameterized light distribution and arcs, and the computation time increases significantly with increasing size of the PSF. The cropped PSF is also used to model the lens light. Subsampling of these two PSFs is required when the pixel sizes are large relative to the full width at half maximum of the PSF in order to avoid numerical inaccuracies caused by an under-sampled PSF². The PSF code subsamples the PSF by a custom factor n such that its pixel scale is $\frac{1}{n}$ of the original pixel size. The two subsampled PSFs are renormalized by the PSF code.

Masks GLEE requires two masks for the modeling. The so-called arc mask contains the region of the lensed quasar and its host (red regions in Fig. 1). It is used during the modeling of the lens light to exclude quasar and host light from the modeling, and it later gives the pixel region that is used to reconstruct the distribution of the source surface brightness (SSB). The lens mask specifies the pixels of the light of surrounding luminous objects such as neighboring stars (green areas in Fig. 1). GLEE sequentially fits the pixels of the lens light, excluding regions of the arc mask and neighboring objects specified by the lens mask, as we describe in the following subsections. The masks are created by the user, who marks the corresponding regions, saves the DS9 region file of the marked areas, and runs a mask-generation code on these DS9 region files.

Fig. 1 Lensing system and regions for mask creation. Top: image cutout of J2100-4452 to illustrate its different components which are modeled with our automation code. Bottom: image cutout of J2100-4452 with regions provided by the user to create the masks. The region required for creating the arc mask, which covers the light of the lensed quasar and arcs, is plotted in red. For the lens mask, the user provides the region(s) of the surrounding luminous objects (green) that contain the pixels that are ignored by the model.

Fig. 2 Overview flowchart for single-band modeling of a strongly lensed quasar system with our new automated code. The modeling consists of two main steps: the lens mass distribution modeling with source or image positions, and the reconstruction of the source surface brightness (SSB; which is the unlensed light distribution of the quasar host galaxy). Each step is divided into smaller modeling steps. The details are described in flowcharts for each step and in the corresponding subsection, as referenced in the flowchart.

2.1.2 Modeling the lens light

The first modeling step is to reconstruct the light of the lens galaxy and the multiple quasar images in order to subtract their contribution and reveal the underlying arc light more prominently. We model the lens light with the Sérsic profile (Sérsic 1963), which is parameterized as

$$I_{\text{S}}\left(x,y\right)=A_{\text{S}}\exp \left[-k\left\{{\left(\frac{\sqrt{{\left(x-x_{\text{S}}\right)}^2+{\left(\frac{y-y\text{s}}{q\text{s}}\right)}^2}}{r_{\text{eff}}}\right)}^{\frac{1}{n_{\text{s}}}}-1\right\}\right].$$(4)

An overview of the parameters in the Sérsic profile and their priors that were used in the automated modeling is provided in Table 3. The constant k normalizes r_eff such that r_eff is the halflight radius in the direction of the semi-major axis.

The Gaussian prior on the centroid coordinates is used to stabilize the optimization during lens light modeling, and it is later changed to a flat prior when the code models the light of the arc (Sect. 2.1.5). The prior range on the axis ratio q_S and Sérsic index n_S are relaxed during multi-band modeling as more constraints are available when additional filters are added. As described in Sect. 2.1.1, we do not subsample the PSF for images with pixels smaller than 0.05" and use the cropped and subsampled PSF for images with a pixel scale greater than 0.05".

The parameters are fit to the observed intensity value $I_i^{\text{obs}}$ of pixel i by minimizing the ${\chi }_{\text{SB}}^2$ of the surface brightness (SB)

$${\chi }_{\text{SB}}^2={\displaystyle \sum _{i=1}^{N_{\text{p}}}\frac{{\left|I_i^{\text{obs}}-\text{PSF}\otimes I_{\text{S,}i}\right|}^2}{{\sigma }_{\text{tot,}i}^2}},$$(5)

where N_p is the number of pixels in the cutout that are not masked, ${\sigma }_{\text{tot},i}^2$ is the total noise of pixel i from Eq. (3), and the crossed circle represents the convolution of the PSF with the predicted Sérsic intensity I_S,i of pixel i. The modeling code runs cycles that consist of one simulated annealing in the beginning and then one MCMC chain using a sampling covariance matrix obtained from the previous chain, which is updated after every cycle. To obtain the uncertainties and optimized parameter values, we sample the parameter space with an MCMC method. Each chain samples 100000 points, and the code removes the burn-in phase, which we conservatively set to be the first 50% of the points in the chain. The modeling runs until two criteria are met. The first criterion is

$${\chi }_{\text{red}}^2=\frac{{\chi }_{\text{SB}}^2}{\text{DOF}}\le 1,$$(6)

with the degrees of freedom (DOF), which is the number of pixels in the cutout minus the number of masked pixels minus the number of free parameters. The number of masked pixels is the sum of all pixels inside the arc mask and of nearby objects masked with the lens mask. The second criterion is approximate convergence of the chain, that is,

$$\text{Δ}\log P\le 5.$$(7)

This is the difference in the log-likelihood between the median of the first 2000 points in the chain after burn-in is removed and the median of the last 2000 points. We visualize the lens light modeling steps with a flowchart in Fig. 3.

If the chain of the lens light model with one Sérsic profile (approximately) converges (i.e., fulfills Eq. (7)) without reaching the ${\chi }_{\text{red}}^2$ threshold in Eq. (6), the code checks whether a point source at the center of the lens improves the model given the possibility of an active galactic nucleus (AGN) within the lens galaxy. The light distribution of a point source is represented by the PSF. For this, the code checks how different point sources in an amplitude range between 0.01 and 100 change the ${\chi }_{\text{SB}}^2$ of the model. If ${\chi }_{\text{SB}}^2$ is lowered by more than 10, a point-source component with the corresponding amplitude is added to the model. This is motivated by the Bayesian information criterion (BIC; following Rusu et al. 2019): the point source adds one additional parameter (the amplitude). A ${\chi }_{\text{SB}}^2$ decrease by 10 ensures that the BIC does not increase. The model is subsequently optimized and sampled until the criterion in Eq. (7) is fulfilled.

When the criterion defined with Eq. (7) is met, the code checks the criterion of Eq. (6). If the latter is fulfilled, the code stops the lens light modeling. When satellites are present in the cutout field, the code continues to model the light of the satellites. If the criterion of Eq. (6) is not fulfilled, one additional Sérsic profile with the same centroid is added to the model for the primary lens light. The code runs optimizing and sampling cycles until the criterion defined by Eq. (7) is met. We model the light of the satellites with Sérsic profiles when the chain of the model for the primary lens light has converged according to Eq. (7). The user is queried to replace the masks such that the satellites are no longer masked. The model is then optimized and sampled until the criterion in Eq. (7) is fulfilled. If this converged model of the primary lens light and of the satellite galaxy light still does not fulfill Eq. (6), the user is asked to either add a third Sérsic profile, refine the masks, or directly continue with the next modeling step. In the first two cases, the code again runs an optimizing/sampling cycle until the AlogP criterion (Eq. (7)) is met.

2.1.3 Modeling the quasar light

After the lens light modeling procedure (see Sect. 2.1.2 and Fig. 3), that is, when we have obtained a good fit for the lens and the satellite light, the code continues by modeling the light of the lensed quasar. Now we mask only the luminous objects surrounding the lens system, that is, we remove the arc mask such that the model now includes the lensed quasars and arcs. For each quasar image, we add a point-like light component that is represented by the PSF and is initially centered on the positions given by the user. The PSF used here was subsam-pled for images with pixels >0.05" and extended ~4"×4" to include the diffraction spikes. We only allowed the quasar positions and amplitudes to vary; the Sérsic lens light parameters were fixed. To infer the best-fit parameters for the quasar images, we again used simulated annealing and an MCMC method with a sampling covariance matrix obtained from an earlier MCMC chain. We stopped the quasar light modeling when approximate convergence of the chain (Eq. (7)) was achieved.

2.1.4 Modeling the lens mass distribution with source or image positions

With the quasar image positions obtained from the quasar light modeling (see Sect. 2.1.3), we can model the mass distribution of the lens galaxy using the mapped source positions and subsequently the image positions. The modeling code queries the user for the choice of the lens mass profile. Our automated modeling code supports both a pseudo-isothermal elliptical mass distribution (PIEMD; Kassiola & Kovner 1993) and a singular power-law elliptical mass distribution (SPEMD; Barkana 1998). The steps conducted in the modeling of the lens mass distribution with source and image positions are visualized in Fig. 4.

The dimensionless surface mass density (or convergence) of the SPEMD profile is given by

$${\kappa }_{\text{SPEMD}}\left(x,y\right)=\left(1-{\tilde{\gamma }}_{\text{PL}}\right){\left[\frac{{\theta }_{\text{E}}}{\sqrt{q_{\text{m}}{\left(x-x_{\text{m}}\right)}^2+\frac{{\left(y-y_{\text{m}}\right)}^2}{q_{\text{m}}}}}\right]}^{2{\tilde{\gamma }}_{\text{PL}}},$$(8)

where (x_m, y_m) is the lens mass centroid, q_m is the axis ratio of the elliptical mass distribution, r_c is the core radius, which we set to 1 × 10⁻⁴, and ${\tilde{\gamma }}_{\text{PL}}=\left({\gamma }_{\text{PL}}-1\right)/2$ is the radial profile slope (where γ_PL is the power-law slope of the three-dimensional mass density $\rho \left(r\right)\propto r^{-{\gamma }_{\text{PL}}}$.

The convergence κ is related to the lens potential via ∇²ψ(θ) = 2κ. The scaled deflection angle for computing the deflection of light rays is the gradient of the lens potential: α(θ) = ∇ψ(θ).

An overview of the mass parameters and their priors is shown in Table 4. The isothermal PIEMD profile is a special case of the SPEMD profile, where ${\tilde{\gamma }}_{\text{PL}}=0.5$ = 0.5 holds (although with slightly different parameterizations for the κ) in Eq. (8). In case of the SPEMD profile, our uniform modeling pipeline first keeps the power-law index fixed to the isothermal case ${\tilde{\gamma }}_{\text{PL}}=0.5$. Only later, when the code models the light of the arc (Sect. 2.1.5), do we allow ${\tilde{\gamma }}_{\text{PL}}$ to vary between 0.3 and 0.7.

The lens mass parameters are first constrained by modeling a source position that reproduces the observed image positions of the quasars. First, the position of the source is obtained by using the lens equation

$${\beta }_i={\theta }_i-{\alpha }_i\left({\theta }_i\right),i=1,\dots ,N_{\text{im}},$$(9)

where the image positions θ_i = (x_QSO,i, y_QSOi,) are mapped back to the source plane with the scaled deflection angle α_i(θ_i). N_im is the image multiplicity of the system. The lens mass parameters η are then varied to minimize the quantity

$${\chi }_{\text{sr}}^2={\displaystyle \sum _{i=1}^{N_{\text{im}}}\frac{{\left|{\beta }_i\left(\eta ,{\theta }_i\right)-{\beta }^{\mathrm{mod}}\left(\eta ,{\theta }_i\right)\right|}^2}{{\left(\frac{{\sigma }_i}{\sqrt{{\mu }_i}}\right)}^2}},$$(10)

with β_i the mapped source position for image i, β^mod the modeled source position (i.e., the mean of the mapped source positions, weighted by the magnification), σ_i the positional uncertainty on the image plane, and μ_i the magnification of image i. We adopted 0.004"as the positional uncertainty, which accounts for both the astrometric uncertainties from PSF fitting in Sect. 2.1.3 and for astrometric perturbations of typically several milliarc-seconds (mas) due to mass substructures in the lens galaxy (e.g., Chen et al. 2006). For HST images where the PSF is well sampled, the centroid position of the PSF can be measured with 2-4 mas for galaxy-scale lenses (e.g., Lehar et al. 1999; Sluse et al. 2012), which is substantially more precise than the pixel size; in Sect. 5.3 we show that our astrometric uncertainties from PSF fitting reach an accuracy of 2 mas. In order to minimize${\chi }_{\text{sr}}^2$, we used simulated annealing to optimize the parameters until ${\chi }_{\text{sr}}^2<5$.

Based on the model obtained by minimizing the ${\chi }_{\text{sr}}^2$, which is based on the mapped source positions, the parameter values are further optimized by using the image positions as constraints. When a satellite galaxy is located close to the lensing system, we can model its mass distribution with a singular isothermal sphere (SIS) profile. The satellite mass centroid is fixed to the satellite light centroid obtained in Sect. 2.1.2. In addition, the user can choose to add an external shear profile to the model after the code has constrained the lens mass parameters with the modeled source position. This accounts for the tidal gravitational field of objects surrounding the lensing system. The shear magnitude is calculated by ${\gamma }_{\text{ext}}=\sqrt{{\gamma }_{\text{ext},1}^2+{\gamma }_{\text{ext},2}^2}$ with γ_ext,1 and γ_ext,2 the components of the shear matrix. The parameters η are varied to minimize

$${\chi }_{\text{im}}^2={\displaystyle \sum _{i=1}^{N_{\text{im}}}\frac{{\left|{\theta }_i^{\text{obs}}-{\theta }_i^{\text{pred}}\left(\eta ,{\beta }^{\mathrm{mod}}\right)\right|}^2}{{\sigma }_i^2}},$$(11)

with ${\theta }_i^{\text{obs}}$ and ${\theta }_i^{\text{pred}}$ the observed and the predicted position for image i, respectively. Again, the code optimizes with simulated annealing until ${\chi }_{\text{im}}^2<5$. It then samples with MCMC chains until approximate convergence (Δlog P ≤ 5).

Fig. 3 Flowchart showing the decision tree for the lens light modeling. The modeling code obtains the best configuration of Sérsic and point-source profiles to fit the primary lens light. The quality check is made with a ${\chi }_{\text{red}}^2$ approach and a criterion for stability in the posterior values. In addition, the code also has the option of modeling the light of satellite lens galaxies.

Fig. 4 Flowchart showing the decision tree for modeling the lens mass distribution with source or image positions. After the user specifies a lens mass profile (PIEMD or SPEMD), the code models the lens mass distribution first by using the modeled source position, and then by using the positions of the multiple quasar images with the option to include an external shear component to the model. Both steps use simulated annealing until the respective χ2 is below 5. The image position modeling is completed after an approximately converged MCMC chain is achieved. After this, the code continues the lens mass modeling with the arc light modeling and reconstruction of the SSB distribution.

2.1.5 Modeling the lens mass distribution with surface brightness distribution

When the lens mass distribution is modeled by using the multiple quasar image positions, we proceed to constrain the lens mass parameters by modeling the SB distribution of the unlensed source and the light of the arc, which is the lensed light of the quasar host galaxy. We illustrate the following modeling steps with a flowchart in Fig. 5.

Our automated modeling pipeline uses GLEE to reconstruct the most probable SSB distribution on a pixelated grid through a linear inversion of the pixel intensities in the arc mask (Suyu et al. 2006). The optimizing or sampling of the posterior probability distribution includes both the ${\chi }_{\text{SB}}^2$ term in Eq. (5) and a prior term that acts on the source pixels with a quadratic regularizing function. The SSB distribution is then mapped back to the image plane to reconstruct the arc light. We use the properties of the arc light to further constrain the lens mass parameters. The inferred quasar light intensity in the image cutout (see Sect. 2.1.3) is a superposition of the real quasar amplitude and a roughly constant contribution from the lensed host light. To reconstruct the light of the host galaxy, we lower the parameter value of the quasar amplitude. The new amplitudes are determined such that the condition

$$I_{\text{circ,avrg}}=3\times I_{\text{arc,avrg}}$$(12)

is fulfilled, where I_cir,avrg is the average intensity in a circle with a radius of 4 pixels around the center of the quasar after the quasar and lens light is removed, and I_arc,avrg is the average intensity of a region provided by the user that contains a small patch of arc light. The factor of 3 is chosen empirically such that the flux is continuous across the arc. Our automation code then adjusts each quasar amplitude such that the criterion in Eq. (12) is met. For the optimizing or sampling process, we first allow the lens mass parameters, now also including the power-law index ${\tilde{\gamma }}_{\text{PL}}$ of the SPEMD profile, and the quasar amplitudes to vary while fixing all the parameters of the lens or satellite light profiles and the quasar centroids. From this point on, we remove the Gaussian prior of the lens galaxy centroids in the lens mass and Sérsic profiles.

To sample the parameter space efficiently, we use EMCEE, an ensemble sampler based on MCMC chains that is highly parallelizable (Foreman-Mackey et al. 2013). With EMCEE, we obtain a sampling covariance matrix for the parameter space, which makes the subsequent sampling with Metropolis-Hastings MCMC more efficient. We consider this modeling part successful when the Metropolis-Hastings chain has approximately converged, now with a stricter requirement:

$$\text{Δ}\log P\le 2.$$(13)

The quasar images are much brighter than the light of the arc, and residuals of quasar light fitting with imperfect PSF can often overwhelm the signal in the arcs. Therefore, we allow for larger uncertainties in the quasar image areas by boosting the error map in these regions. The user marks the areas dominated by quasar light, and the code will boost all the pixels of the error map within this region such that their normalized residuals are within ±1σ. The error map can also be boosted in the satellite galaxy region if the user chooses to do so. This can be necessary if the satellite galaxy is very close to the arc of the lensing system. In the next modeling step, now including the boosted error map, we fix all quasar light parameters and allow the lens mass and the lens light parameters to vary. If the uncertainties in the satellite region are not boosted, the satellite light parameters are varied as well. We again use EMCEE to obtain a sampling covari-ance matrix and optimize or sample with simulated annealing or an MCMC until approximate convergence of the chain (Eq. (13)) is achieved. In a last step, the code checks the model for full convergence by analyzing the discrete power spectrum of the chain (following Dunkley et al. 2005). If full convergence is not yet achieved, the code runs optimizing and sampling cycles with an increased chain length by now sampling 200 000 points. When full convergence is achieved, the code stops and the modeling is completed.

2.2 Multiband modeling

In addition to the single-band automation code described in Sect. 2.1, we developed a code to simultaneously model systems with strongly lensed quasars using multiband data. The modeling procedure in most parts is analogous to the single-band case and is described in Table 5. We present an overview flowchart in Fig. 6.

When we obtain a single-band lens mass model in one band (the reference band) and a lens and quasar light model in the other bands with the pipeline described in Sect. 2.1, we proceed to jointly model the lens and quasar light in all the considered bands. We chose the F 160W band as our reference band unless otherwise specified because the arc light is most prominent in this wavelength. We use the modeled quasar positions of the individual bands to align the multiple filters. The lens mass distribution model in the reference filter serves as a starting model for the multiband modeling.

In the first step, we pair the reference filter with one of the additional filters and model the lens and quasar light together with MCMC chains. We fix the amplitudes of the lens light components in the reference filter, and force the lens light amplitudes in the other filter to follow the same ratio of light amplitudes as in the reference filter. In addition, we link the remaining lens light parameters (centroid, axis ratio, position angle, effective radius, and Sérsic index) between the filters. The quasar light centroids are linked as well, while the quasar amplitudes vary independently. After a MCMC chain is approximately converged with ΔlogP ≤ 5, we fix the lens light amplitudes of this filter, and add the next filter, which is modeled together with the reference filter in the same way as the filter before. We also allow variation in the quasar light parameters of the previous filters. When approximate convergence (ΔlogP ≤ 5) is achieved for all pairs, we model the lens and quasar light parameters of all filters simultaneously and now allow all lens light amplitudes to vary independently while keeping the remaining lens light parameters and quasar light centroids linked.

Like in the single-band modeling, the modeled quasar positions are then used for image position modeling (see Sect. 2.1.4) to constrain the lens mass parameters. In the last step, the light of the arc and the SSB distribution are reconstructed analogously to the single-band modeling case described in Sect. 2.1.5. Finally, we run additional MCMC chains to achieve full convergence of our lens model parameters.

Fig. 5 Flowchart for arc light modeling and SSB reconstruction. In the first stage, the quasar amplitudes are lowered, and in the second stage, the error map is boosted in the bright quasar regions in order to model the light of the arc. The lens mass parameters are allowed to vary in this modeling step. In a final step, the code runs optimizing and sampling cycles until full convergence of the chain is achieved.

2.3 Time-delay estimates

Strong-lensing time-delay cosmography requires time-delay measurements and thus observational programs that monitor the lensed quasars. To facilitate the scheduling of the monitoring, estimations of the time delays of new quasar systems are useful. Furthermore, a comparison of the (blind) prediction of the time delays to the subsequently measured time delays provides a crash test of mass modeling procedures (Treu et al. 2016).

When we obtain a good model for the lens potential ψ(θ) of the system from the previous subsections, we can use it to compute the relative time delays between the different images. Measured redshifts of the lens (z_d) and of the source (z_s) allow us to convert the dimensionless surface mass density into a physical quantity. We can then calculate the time-delay distance

$$D_{\text{Δ}t}=\left(1+z_{\text{d}}\right)\frac{D_{\text{d}}{D}_{\text{s}}}{D_{\text{ds}}}$$(14)

by assuming a cosmological model. The distances on the right-hand side D_d, D_s and D_ds are the angular diameter distance to the deflector or lens galaxy, to the source galaxy, and between the deflector and source galaxy, respectively. With the time-delay distance, GLEE can calculate the time delays between the quasar images i and j,

$$\text{Δ}{t}_{ij}=\frac{1}{c}{D}_{\text{Δ}t}\text{Δ}{\tau }_{ij},$$(15)

where

$$\text{Δ}{\tau }_{ij}=\left[\frac{1}{2}{\left({\theta }_i-\beta \right)}^2-\psi \left({\theta }_i\right)\right]-\left[\frac{1}{2}{\left({\theta }_j-\beta \right)}^2-\psi \left({\theta }_j\right)\right]$$(16)

is the Fermat potential difference between quasar images i and j, which is obtained from the final lens mass model. We have a separate code that calculates the time-delay distance from Eq. (14) and the relative time delays from Eq. (15), including 1σ uncertainties by using the final MCMC chain of the model parameters from the modeling code and the cosmological parameters and redshifts provided by the user.

3 Observations

Our sample consisted of the nine strongly lensed quasars DES J0029-3814 (Schechter et al., in prep.), DES J0214-2105 (Agnello & Spiniello 2019; Lee et al. 2019), DES J0420-4037

(Ostrovski et al., in prep.), PS J0659+1629 (Delchambre et al. 2019), 2M1134-2103 (Lucey et al. 2018; Rusu et al. 2019), J1537-3010 (Lemon et al. 2019; Delchambre et al. 2019; Stern et al. 2021), PS J1606-2333 (Lemon et al. 2018), PS J1721+8842 (Lemon et al. 2018), and DES J2100-4452 (Agnello & Spiniello 2019). A detailed description of the observations and the peculiarities of each lensing system are presented by S22. In Fig. 7 we show a color image that was created with the three HST bands F160W, F475X, and F814W, which were used as the red, blue, and green channels, respectively.

4 Modeling results

In this section, we describe the modeling results of our sample of lensing systems obtained with the autonomous modeling pipeline that we presented in Sect. 2. The input files for each lensing system, that is, the PSF, error map, and masks, were obtained by following Sect. 2.1.1. Each system was modeled with a SPEMD profile and external shear. The lens light was modeled with two Sérsic profiles unless specified otherwise in the following subsections. For two systems (J0659+1629 and J1606-2333), we additionally modeled the light and mass distribution of a close-by satellite galaxy. We decided which bands to include by visual inspection of the residuals after modeling lens and quasar light, that is, the remaining light from the arcs, which are the lensed host galaxy of the quasar. If there was significant arc light in one wavelength band, then we included that band in the final multiband model.

In Table A.1 we present figures of the modeled light and reconstructed source for each of the nine lensing system in our sample. We show the observed image (third column), the modeled light and critical curves (fourth column), and the normalized residuals in a range between −5σ and 5σ (fifth column). For each pixel, the normalized residuals show the difference between data and model, normalized by the estimated standard deviation. We show cropped images instead of the full cutout that we used for modeling for better visibility and indicate 1" with a white line. The panels are oriented such that north is up and east is left. The sixth column shows the reconstructed SSB distribution of the quasar host galaxy on the image plane. The seventh column shows the SSB distribution on the source plane with plotted caustic curves in red, and the mean weighted source position of the quasar as a blue star. We show rulers with length of 0.5" in x- and y-direction because the pixel size can be different in these directions.

In Table B.1 we present the median parameter values of the lens mass distribution together with their 1σ uncertainties from the final chain of each multiband model. We also show the best-fit lens and satellite light parameters in the F 160W band (the reference filter). The centroid coordinates are given with respect to the bottom left corner of the image cutout. The light parameters in the F475X and F814W bands are presented in Tables B.2 and B.3, respectively.

In Table 6 we summarize the single-band ${\chi }_{\text{red}}^2$ of each modeled band and a total ${\chi }_{\text{red,tot}}^2$, which is the sum of the independent single-band χ2 divided by the DOF of the multiband model.

We excluded pixels masked in the lens mask and the boosted quasar light regions in computing the DOF (see Sect. 2.1.2). For two systems (2M1134-2103 and J1606-2333) the ${\chi }_{\text{red}}^2$ is > ${\chi }_{\text{red}}^2$ of 1.84 and 1.36, respectively). For these systems, the uncertainties could be underestimated. In Table C.1 we present the convergence κ, total shear strength γ_tot, and magnification $\mu =1/\left({\left(1-\kappa \right)}^2-{\gamma }_{\text{tot}}^2\right)$ at the (modeled) image positions.

In addition, we calculated the Fermat potential differences at the multiple modeled image positions and predict the relative time delays for each system. We used the redshifts mentioned in Sect. 3. If the lens redshift was not yet available, we assumed z_d = 0.5, which is an approximate median value based on the sample of currently known lensed quasars with redshift measurements. The results are shown in Tables C.2 and C.3.

Because the lens light is described by multiple Sérsic components, we calculated the second brightness moments of the primary lens light model to compare the median centroid, axis ratio, and position angle to that of the mass distribution. The result is shown in Fig. 8. The mass and light centroid align well within one pixel (0.08") in most cases, with the exception of PS J0659+1629, which shows a large offset in the x-centroid. The position angle of mass and light agree mostly well within ~10° for five of the nine systems. The largest offset is also reached for J0659+1629 with ~50°. For the axis ratio, the models do not follow the 1:1 relation so closely. Five models tend to be more elliptical in mass compared to the light, where four of them either have a satellite or are highly sheared. This may explain the low mass-axis ratio.

Individual modeling details for each lensing system are discussed further in the following subsections.

Fig. 6 Overview flowchart for multiband modeling of a lensed quasar system with the multiband automated code. Details are described in flowcharts for each step and in the corresponding subsections, as referenced in the flowchart.

Fig. 7 Color images for each of the nine lenses in our sample created with the three HST bands F160W (red), F 475X (blue), and F814W (green).

4.1 DES J0029-3814

This lensing system was modeled solely in the IR F160W band. In the remaining bands, sufficient arc light from the host galaxy could not be identified by eye. We used two Sérsic profiles and a central point-source component for the light from the primary lens. The PSF was estimated by selecting two stars in the field. The best-fit parameters in Table B.1 show that the centers of the mass and light distribution are offset by ~0.04" in the x-direction and by ~0.02" in the y-direction. We obtain a high shear magnitude of ~0.2 with a low position angle (14°), which may be due to the overdense environment north of the lensing system.

4.2 DES J0214-2105

We modeled this system in all three considered bands (F160W, F475X, and F814W) and note that the quasar amplitudes in the F475X band are strongly underfit after step 5(a) (see Table 5). The reason might be that the faint arc light and an imperfect PSF model result in higher residuals in regions of multiple quasar images. The underfitting is not directly obvious from the residuals because the missing quasar amplitude is compensated for by high source intensity values in a single pixel. To avoid negative impact of this band on the lens mass model, we present the F160W single-band model. The PSF was approximated with one star from the field in the F160W band and two stars in each of the F475X and F814 bands. The centers of the mass and light distribution are very well aligned, with an offset within ~0.01".

4.3 DES J0420-4037

This lensing system was modeled in all three bands. As in the case of J0214-2105, the quasar amplitudes in the F475X band are strongly underfit. We present the F160W single-band model results. We approximated the PSF in the F160W band by choosing two stars in the field; for the F475X and F814 bands, we used one star. Mass and light show a small offset of <0.01 in the x-direction and ~0.01" in the y-directions. The modeling residuals indicate several light clumps near or in the arc. We confirmed two of them to be counter images, which indicates a second lensed background source.

4.4 PS J0659+1629

This lensing system was modeled in the IR F160W band because the arc light in the other bands was insufficient. We also included an SIS profile for the satellite mass and a Sérsic profile for the satellite light. The satellite is very close to the northern-most quasar image. We therefore boosted the error map in the satellite region in the same way as the quasar images in the last modeling step. The PSF was estimated with four stars from the field. We found two very faint light clumps in the residuals, which are at positions that suggest them to be counter images of a second, lensed background source at a similar redshift as the first source. The dark region in the source light reconstruction at the satellite position does not have an effect on the mass model because its counter image lies outside the reconstruction region (arc mask). For this system, the quasar amplitudes are also underfit after step 5(a) (see Table 1), which might have an effect on the lens mass model. This is difficult to overcome with automated modeling procedures and require individual tweaks that are deferred to future work for cosmographic analysis. This system shows the largest difference between mass and light in our sample, with an offset of ~0.1" in the x-direction and ~0.05" in y-direction. Furthermore, the modeled lens mass distribution is much more elliptical, with Δq = q_S − q_m ~ −0.4, and the position angle is offset by −45°. The reason might be the satellite galaxy close to image D.

4.5 2M1134-2103

This system has a faint lens and very bright quasar images, which makes an exact modeling of the lens light and spectroscopic redshift of the lens difficult. We performed the modeling in the F160W and F814W bands. The PSF for the F160W band was estimated by using five stars. We chose three stars to build the PSF in the F814W band. The median mass and light centroid align very well with an offset smaller than ~0.01" in the x-direction and ~0.01" in the y-direction. This system has the highest external shear magnitude γ_ext of the sample. This is consistent with the elongated shape and the overdense environment in the field of view of this system (Rusu et al. 2019).

Fig. 8 Comparison between the position and structural parameters of the lens mass and light for eacht lensing system. Top: difference between the median mass and light centroid of our models, with Δx = xs − xm and Ay = ys − ym. Bottom left: comparison between the median mass and light position angle. Bottom right: comparison between the median mass and light-axis ratio. The dashed lines show lines without a centroid offset (top) and the lines of equal position angle and axis ratio (bottom left and bottom right).

4.6 J1537-3010

This system was modeled in all three bands. In each band, we chose five stars from the field to approximate the PSF. The centers of the mass and light distribution align very well within ~0.01". In the model residuals shown in Table A.1, we find a small offset between the quasar images and the arc in bands F475X and F814W. A possible explanation for this offset is that these two bands are in the rest-frame UV of the quasar host galaxy, and thus are likely dominated by light from specific star-forming regions in the host galaxy. The power-law slope ${\tilde{\gamma }}_{\text{PL}}~0.3$ hits the lower bound of the prior, which might be due to our imperfect PSF models in the F475X and F814W models. As we show in Sect. 5.2, the PSF has a significant effect on ${\tilde{\gamma }}_{\text{PL}}$. This system has one of the most prominent arcs in all three HST filters in our sample, which is ideal for future cosmographic analysis.

4.7 PS J1606-2333

This lensing system was modeled in the F160W band with two Sérsic profiles and a central point source component for the primary lens light. Although we initially included the other bands because they show significant arc light, the source reconstruction shows no visible source (above the noise level) in these bands. We excluded these bands and optimized only the single-band model. For this lensing system, the PSF was stacked from five neighboring stars. The bright light clump in the south below image B that is embedded in a faint extended structure. We assumed this clump to be a satellite at the same redshift as the lens galaxy, and include an SIS profile for the satellite mass and a Sérsic profile for the satellite light. The error map was boosted at the satellite position in the last modeling step because of its proximity to the arc. Mass and light align very well within ~0.01". The power-law slope of this system is close to the lower bound of the prior with ${\tilde{\gamma }}_{\text{PL}}~0.31$. Again, this might be due to an imperfect PSF model.

4.8 PSJ1721+8842

This is an unusual and interesting lens system with six quasar images. It is composed of two quasar sources at the same red-shift, with one quasar lensed into a quad (images A, B, C, and D) and the other one lensed into a double (images E and F; Lemon et al. 2022). Most of the arc light in this system comes from the double system, therefore we only included these regions for the arc light modeling and SSB reconstruction. We modeled this system in all three considered bands with two Sérsic profiles and a central point-source component for the primary lens light. We picked three stars from the field to approximate the PSF in the F160W band, and one star in each of the F475X and F814W bands. The second source component was added manually because the current version of the automation code does not support an automated detection of multiple source components. It is represented as a green star in the source reconstruction plot in Table A.1. Quasar image F is offset relative to the arc peak SB, which may be due to the quasar being kicked out of its host galaxy (Mangat et al. 2021). We used F475X as the reference band because image F is not clearly distinct from the arc in the F160W band where the arc is the brightest. The offset introduces a challenge to the quasar and arc light modeling. At the position of the second source component, the source intensity is negative. The mass and light centroids are offset by -0.04" in the x-direction, and align very well within -0.01" in the y-direction.

4.9 DES J2100-4452

The final model of this system includes only the F160W band. The F475X and F814W bands show small indications of arc light, but including these bands in the model led to no visible source reconstruction. The PSF was created by using five stars. A dust lane appears to cross the lens galaxy from north to south. Its influence on the model needs to be investigated in future work. The dust lane can be seen best in the bluest band (in our case, the F475X band), and given its rest frame wavelength of 1.3 microns in the F160W band, the effect on the model is probably very weak. The mass is offset relative to the light by ~0.04" in the x-direction and by ~0.06" in the y-direction.

5 Discussion

5.1 Comparison of automated modeling to manual modeling with GLEE

In this section, we compare the automated modeling with our automation code and the manual modeling of a lensing system with GLEE based on the total time actively spent by the user on modeling a lensing system. We report a significant reduction of preparation time for the input files, which are created in separate automated codes. The creation of the image cutout is sped up from ~15 to ~5 min in the automated modeling. The time for creating the error map stays roughly the same because this was already automated prior to this work. The mask regions in both cases were obtained by eye and were marked manually, but the masks are now created fully automatically from the region files. This automation halves the user input time from roughly one hour in the manual case to -30 min when the automation code is used. The time for creating the PSF is significantly reduced from ~2 h to ~15 min in the case of the automated modeling. When the PSF code is used, the only step the user has to conduct in order to obtain the PSF is to find suitable stars in the field and save their positions. The time-consuming technical steps such as saving the region file, stacking the files, subsampling if requested, and normalizing to create the PSF, are fully automated.

For the actual lens modeling process, we can also report a significant time reduction. The initial setup of the GLEE configuration file, which also includes estimations for the starting values of the parameters, is about one hour for manual modeling. In the automated case, the user only has to provide a region file, which can be created within 10 min. The setup of the configuration file and the estimations for the starting values are fully automated, so that it is not necessary for the user to learn how to use GLEE (e.g., how to set up the GLEE configuration file). Quality control and modifications to the configuration file between different modeling steps are also fully automated, which again saves several hours over the full modeling process. Moreover, there is no waiting time between the multiple MCMC chains, and the code obtains the optimal sampling parameters automatically. Together with short interaction between the user and the code and the inspection of the output, we estimate an average of one hour of total user time per lensing system for the automated modeling process. Table 7 gives a comparison between user input time in the automated and in the manual modeling case.

The computational time for the arc light modeling and SSB distribution reconstruction mainly depends on the arc mask and PSF size. The optimizing and sampling cycles with simulated annealing and MCMC chains use one core at a time. We used EMCEE, which can be highly parallelized, to obtain an initial sampling covariance matrix at the beginning of the last modeling steps (5(a) and 5(b) in single-band, and 4(a) and 4(b) in multi-band modeling) with a chain length of 400 000 on a cluster with Intel Xeon E5-2640 v4 CPU using 30 cores. The average computation time on the cluster for each EMCEE chain is -5 h. Because of typical overhead in parallel computing, we estimate that the single core computation time for each EMCEE chain is ~100 h. EMCEE sampling is used twice in the automated modeling procedure, and thus we have a total EMCEE computation time of ~200 h per core. With this setup and modeling specifications, we estimated the average computation time per core per lens system for each modeling step described in Sect. 2, and summarize them in Table 8. The computationally expensive step of arc light modeling and SSB distribution reconstruction takes 15–20 days with a single core. Because of the parallelization of EMCEE, the effective computation time ranges between 7 and 14 days. When only approximate (and not full) convergence of a chain is sufficient, models can already be expected after 5–7 days. These models do not differ much from the final model that is fully converged and can be used for tentative analysis. We expect the computation time to be roughly the same in both the manual and the automated modeling.

5.2 Comparison of lens modeling results between GLEE and Lenstronomy

5.2.1 Blind comparison

All nine systems of our sample are part of a more extensive sample of 30 quads that were modeled independently by S22 with a similar, automated pipeline based on the lens modeling software Lenstronomy. Both modeling frameworks are dedicated to high-resolution images of lensed quasars, and thus use similar modeling assumptions for the lenses: the lens mass distribution is modeled with a power-law profile with external shear, and the lens light distribution is modeled with Sérsic profiles.

There are important differences between the two procedures that need to be taken into account. First, the two softwares have different implementations for the SSB reconstruction: GLEE uses a pixelated grid of the source intensity distribution, whereas Lenstronomy uses circular Sérsic profiles with additional optional shapelet components. Second, S22 always modeled the three available bands of HST images, while in this work, we focused only on those in which the arcs are clearly visible by eye, typically F160W. Even though these are the most informative bands in our nine-system sample, they do not capture the full information, especially because the resolution is better in the UV/visible (UViS) data. Third, different priors have been chosen for some of the key parameters. S22 adopted informative priors from the SLACS sample (Bolton et al. 2006, 2008; Auger et al. 2010) on the radial slope of the mass density profile and the ellipticity of the mass distribution in order to avoid unphysical solutions when the parameters are poorly constrained by the data. In this study, we opted to adopt uniform priors within some bounds to constrain the parameters solely by the data at hand. In the best cases, the priors should not matter because the likelihood should dominate the posterior, but in practice, as we show below, many of the systems are poorly constrained by the data, and therefore the priors do matter. Fourth, S22 adopted an iterative scheme to improve the PSF estimate based on the multiple images of the quasar, sampled at the scale of the reduced data. In contrast, in this work, we used an subsampled PSF estimated from images of stars in the field, without additional corrections. As shown by Shajib et al. (2022), the PSF is a crucial ingredient for cosmography-grade inference.

Keeping the caveats in mind, we compared the lens mass parameters, external shear parameters, convergence and total shear at the quasar image positions, image magnifications, Fermat potential differences, and predicted time delays. The adopted cosmology is the same as in S22. We plot in Fig. 9 the results of the two modeling codes together with a line showing the identity.

In addition, we present the difference of the median values of both modeling codes against the GLEE parameter values to better illustrate the absolute differences between the final results.

We find excellent agreement in the Einstein radius of the primary lens galaxy, where the median values agree within 4% for eight of the nine systems. This is expected, and it is reassuring that SL provides robust masses enclosed within the Einstein radius. An apparent outlier is J0659+1629 with a difference of ~0.2 , which is due to the presence of a satellite galaxy whose mass is degenerate with that of the primary lens. We calculated a value for the "effective" Einstein radius (the radius of the circle for which the enclosed mean κ is 1) that included the satellite mass and compared it to the radius that was obtained by S22 for this purpose. Our value of ${\theta }_{\text{E,eff,GLEE}}={2.368}_{-0.02}^{+0.02}{}^{″}$ matches very well ${\theta }_{\text{E,eff,Lenstronomy}}={2.368}_{-0}^{+0.01}{}^{″}$ for J0659+1629 (the distribution for θ_{E,eff,Lenstronomy} is heavily skewed such that the median coincides with the 16th percentile). We thus conclude that for all systems, we recover the same Einstein radius to within a root-mean-square (rms) scatter of 1.6% when the proper accounting for satellites is considered.

Flattening, position angle of the mass distribution, and external shear magnitude and direction are strongly degenerate with each other (and to some extent with the slope), and should therefore be looked at holistically. To illustrate this, we show the posterior distributions of the mass parameters and external shear of J0659+1629 in Fig. D.1, where we highlight the strong degeneracy between the Einstein radius of the main deflector and the satellite.

For six out of the nine systems, the position angles of the lens mass distribution match within the errors (1 or 2 σ at most). The outliers are J0659+1629, for which GLEE converges to a mass distribution that is highly flattened, more than the light distribution, and therefore likely not fully correct, while S22 - driven by the prior - converges to a much rounder distribution, and J0214-2105 and J0420-4037, for which the discrepancy is due to a combination of differences in flattening, PA, and external shear. We note that for this system, the slope of the mass density profile is poorly constrained without a prior, suggesting that the information content of the data is not sufficient.

The mass-flattening parameter agrees within ~0.1 (larger than the formal uncertainty, which is therefore underestimated) except for two cases: (i) J0659+1629 discussed above, and (ii) J0029+3814, for which S22 reported a stronger flattening than GLEE, perhaps driven by the prior on the axis ratio of the light distribution. The inferred external shear magnitudes agree within 0.05 (again larger than the formal uncertainty); the two most discrepant systems are J1537-3010 and J2100-4452, for which the slope of the mass density profile is relatively poorly constrained without a prior, and which are therefore likely to have a poor information content of the data. Reassuringly, we find a good match of the shear position angle for the two systems with the highest shear magnitude.

As mentioned above, a main source of difference is the prior of the power-law slope ${\tilde{\gamma }}_{\text{PL}}$, which reverberates through the other parameters. Not surprisingly, the inferred power-law slope ${\tilde{\gamma }}_{\text{PL}}$ of the Lenstronomy models obtained by S22 follows their imposed Gaussian prior, which is centered on a close-to-isothermal value of ${\tilde{\gamma }}_{\text{PL}}=0.54$. In contrast, the slope values of the GLEE models span the uniform prior range of [0.3,0.7], with values closer to the bounds than to an isothermal value. As we discussed before, the difference arises from a combination of two factors: differences in the PSF, and information in the data that is insufficient to constrain the slope.

Not surprisingly, given the differences highlighted above, the inferred quantities κ, γ_tot, and μ at the image positions show significant deviations (see Fig. E.1). We calculated the average relative scatter of these quantities for each lensing system, and find that it ranges from ~3 to >100%, and the statistical uncertainties are underestimated in nearly all cases.

The relative Fermat potential differences and predicted time delays (also shown in Fig. E.1) follow the 1:1 line overall, although they often disagree within the statistical uncertainties. For these two quantities, we removed the outlier image C of J0659+1629 for better visibility in the plot, because its time delay is one order of magnitude longer than all other images. Table 9 lists the relative deviation in Fermat potential differences δ(Δτ)/Δτ_GLEE| with δ(Δτ) = Δτ_GLEE − Δτ_Lenstronomy. We estimated the statistical uncertainties on these relative deviations by symmetrizing the uncertainties on the GLEE and Lenstron-omy parameter values (e.g., on Δτ_GLEE and Δτ_Lenstronomy) and assuming that they follow Gaussian distributions.

This comparison shows us that for only one system (J1721+8842) are the relative Fermat potential differences and predicted time delays within 10%. For seven of the nine systems, the deviation of the relative Fermat potential differences of the multiple images is 30% or more on average. Therefore, for most systems, a more detailed modeling and/or better data are needed to bring all these models to cosmography-grade level. In some cases, the UVIS exposures may provide further constraints on the mass and light profile parameters, even though the source contribution is fainter than in the IR, especially because the resolution is twice as high in the UVIS.

Fig. 9 Direct comparison of lens mass parameter and external shear values between our models and Lenstronomy models from S22.

5.2.2 Post-blind analysis

After completing the blind comparison between the results, we consider the efforts that have been made after unblinding to understand the origin of the discrepancies. To test the influence of the power-law slope ${\tilde{\gamma }}_{\text{PL}}$ on the physical quantities κ and Δτ, we rescaled the GLEE and Lenstronomy values of Δτ with $\Delta {\tau }_{\text{iso}}\simeq \Delta \tau /2{\tilde{\gamma }}_{\text{PL}}$ (e.g., Suyu 2012), and calculated new values with Eq. (8) and ${\tilde{\gamma }}_{\text{PL}}=0.5$, thus assuming isothermal mass profiles. In this case, we obtained a much better agreement between the two modeling results overall, as shown in Fig. E.1. With these assumptions, all κ values now agree within 20%, and five of the nine systems have a scatter ≤10%. The trend is similar for the relative Fermat potential differences at the multiple image positions, where six of the nine systems agree within ~20% (see the last column δ(Δτ_iso)/|Δτ_GLEE,iso| in Table 9). Two systems, J0420−4037 and J0659+1629, now have a higher relative deviation in the multiple images, but their absolute values of Fermat potential differences are small. Therefore, most of the discrepancy in κ and Δτ between the two independent models are due to the different power-law slope values. The two modeling pipelines yield different slope values, as discussed above and as expected, because the priors and PSFs differ.

To distinguish the two effects (priors and PSF), we used the final PSF of the Lenstronomy modeling by S22 in our automated modeling procedure. The creation of the PSF of S22 is different than in our pipeline because it includes iterative PSF correction, while our PSF is simply constructed from stars in the field (see Sect. 2.1.1). Moreover, the GLEE PSF is subsampled by a factor of 3 compared to the pixel scale of the data, while the Lenstronomy PSF uses the pixel scale of the original data.

To ensure that the information content is sufficient (and thus the prior is less important), the two teams focused for this comparison on two systems with high signal-to-noise ratio (S/N) data and visible arcs: J1606-2333 and J2100-4452. We remodeled these two systems in the F160W (IR) band with our GLEE-based pipeline, but used the not subsampled Lenstronomy PSF. We compared the results with those based on the GLEE PSF. We show the comparison in Figs. 10 and E.2. In both figures, we present the GLEE blind results as presented in this work and the new post-blind results using the Lenstronomy PSF. In most cases, the values obtained with the Lenstronomy PSF are closer to the modeling results of S22. The Einstein radius of J1606-2333 is lower with the new PSF because it is degenerate with the Einstein radius of the satellite, which is now larger. The effective Einstein radius of this system remains unchanged. The influence on the power-law slope is evident, as the slope values move much closer to the values of S22, and it is closer to isothermal, without imposing a prior. This also leads to a better agreement of other quantities such as κ or the Fermat potential differences to the Lenstronomy modeling results, as shown in Fig. E.1. We conclude from this that for data with a sufficiently high signal-to-noise ratio, the PSF is a crucial ingredient for accurately inferring the power-law slope, and thus the Fermat potential. Crucially, the PSF can be directly reconstructed from the data, thus improving the accuracy of cosmographic analysis. This result reinforces the necessity of PSF reconstruction, which has become best practice in recent years for cosmographic analyses using GLEE and Lenstronomy by members of our team and collaborators (e.g., Wong et al. 2019; Birrer et al. 2019; Shajib et al. 2019b; Chen et al. 2019). The counterpoint is that data of insufficient quality should not be used for cosmographic analysis, unless an accurate and informative prior is available.

We conducted the same test with a subsampled version of the Lenstronomy PSF. It is subsampled in the same way as the GLEE PSF. The modeling results do not agree as well with the Lenstronomy values, and are closer to the GLEE values. This comparison shows that although subsampling in the IR band is important (Suyu et al. 2012, Shajib et al. 2022), the modeling results depend on the way in which the PSF is subsampled. It should be done during the PSF reconstruction process, and not as a simple interpolation after the corrections of the PSF.

To conclude, we recall that the goal of this work was not to produce cosmography-ready models, but to automate the modeling procedure for a wide range of lens morphologies that can subsequently enable the construction of cosmography-grade models; thus differences in the results of the two modeling teams are expected. Our analysis shows that the origin of the differences can generally be understood, and that the two most significant factors are data quality and accuracy of PSF modeling.

Fig. 10 Change in lens mass parameter and external shear values when using the Lenstronomy PSF with our GLEE automation code for J1606−2333 and J2100-4452.

5.3 Comparison of astrometry

To assess the precision of our astrometry, we compared the relative positions of the multiple quasar images from our model with the modeled positions from S22 and Luhtaru et al. (2021) and the measured positions from the Gaia satellite data release 3 (Brown et al. 2021). Like in our pipeline, S22 obtained positions by modeling the surface brightness distribution, while Luhtaru et al. (2021) used geometric properties in the image plane. The comparison with S22 and Luhtaru et al. (2021) was performed for all nine systems in our sample and the comparison with Gaia for seven of the nine systems. For the remaining two systems, no Gaia data are available.

In all comparisons, J1721+8842 has the strongest offsets for all four images of the quad, which is expected because of its complexity (see Sect. 4.8). When this system is excluded, we obtain the following results with respect to S22: We have an r.m.s. scatter of ~6 mas in the x-direction (RA) and ~5 mas in the y-direction (Dec). The comparison with Luhtaru et al. (2021) reveals an r.m.s. scatter of ~1.7 mas in the x- and y-directions. The r.m.s. scatter for the Gaia comparison is 1.7 mas in the x- and 2.3 mas in the y-direction. The top panel of Fig. 11 shows the difference in quasar image positions for the comparison with S22, Luhtaru et al. (2021), and Gaia for all systems that are in all three samples (seven systems). The images of J1721+8842 are marked in orange as the mentioned outlier. The bottom panel shows the difference in quasar image positions for all systems that are in common with Luhtaru et al. (2021) and Gaia. To give an idea of how small the differences are, we include a box centered around zero offset with dimensions ±5mas. The vast majority of points lie within this box. Although we did not perform PSF reconstruction in our automated modeling pipeline, we have recovered the astrometry of the quasar images within ~2 mas of the Gaia measurements (which are most precise or accurate because Gaia was designed to measure astrometry).

Birrer & Treu (2019) provided an estimate of how the astrometric uncertainties translate into the uncertainty on the H₀ inference. The ~2 mas rms offset to Gaia in our analysis translates into a ~0.2 mas offset in the source plane, which leads to uncertainties on H₀ well below 5%, which was chosen by Birrer & Treu (2019) as an estimate of the total uncertainties (i.e., from modeling the lens mass potential and the contribution from the mass along the line of sight). This means that the astrometric uncertainties in this work do not contribute significantly to the cosmographic error budget, which shows that our automated pipeline can meet the astrometric requirements.

Fig. 11 Relative difference in quasar image positions between our results and S22 (squares), Luhtaru et al. (2021, circles), or Gaia (triangles) for systems that overlap in the multiple samples. The image positions of the outlier J1721+8842 are marked in orange (top) or are excluded (bottom).

5.4 Possible improvements of the automated lens modeling pipeline

To further improve the automated uniform modeling in combination with detailed cosmographic analyses of lenses, we plan to implement several additions to the code. For the uniform modeling approach, we reconstruct the PSF from several stars in the field. This approach introduces inaccuracies that are negligible for uniform modeling, but need to be minimized for a detailed analysis, as discussed in Sect. 5.2. We will automate a method to iteratively update the initial PSF using the multiple lensed quasar images. This PSF correction might also resolve the tendency of the power-law slope parameter ${\tilde{\gamma }}_{\text{PL}}$ to move to the boundaries of the prior.

To speed up the modeling procedure and make it more convenient, we can automate the detection of extended structures in the cutout, for instance, the positions of multiple quasar images and objects that need to be masked. Other automated modeling codes already show that an accurate detection in the scope of uniform modeling is often possible , for example in Savary et al. (2022) and Rojas et al. (2021), who also showed that objects are occasionally misidentified, and the automated modeling has to be stopped. Our current approach of manually identifying objects is thus a compromise to ensure very high accuracy. This works best for a medium-size sample (e.g., a few dozen), but is not applicable for samples of some hundred lenses or more.

The initial background subtraction of the data can be made more accurate by allowing for the selection of multiple sky regions. This is especially important for systems that show a gradient in background flux. Another improvement can be using a p-value threshold for the lens light modeling (Sects. 2.1.2 and 2.2). This would decrease the probability that the ${\chi }_{\text{red}}^2$ criterion (Eq. (6)) is met only due to specific noise realizations.

Other modeling teams have shown that optimizing methods that are based on gradient descent and can be heavily parallelized are superior to conventional methods in terms of computation time, for instance, GIGA-Lens (Gu et al. 2022). These directions are worth exploring for future developments of automated lens modeling.

6 Conclusions

We presented a new automated modeling code for the uniform modeling of strongly lensed quasars in galaxy-scale systems with high-resolution data based on the well-tested software GLEE. We additionally developed several codes that create necessary input files (i.e., PSF, error map, and masks) that are used for the modeling. The lens mass distribution of the system is modeled in two steps. In the first step, the modeling code uses the predicted source position and the observed image positions to constrain the distribution of the lens mass parameters. In the second step, the light distributions of the lens, the multiple lensed quasar images, and their host galaxy are modeled. The latter is used to reconstruct the surface brightness distribution of the source. With the SSB distribution, the code models the light of the arc to further constrain the lens mass parameters. The modeling code and the creation of input files are nearly fully automated, requiring only minimum user input. Quality control and technical steps in the modeling process are fully automated.

We tested the modeling code on a sample of nine strongly lensed quasars and obtained a good light fit and a good alignment of the mass and light distributions. The current automated pipeline is versatile enough to model a variety of lens systems in a uniform manner to obtain lens mass models. The models provide robust estimates of Einstein radii and astrometry of the quasar images. The accuracy of parameters such as the power-law slope depends crucially on the data quality and on the details of the modeling. This is highlighted by a blind comparison of our models with those of an automated modeling framework from S22, based on the lens modeling code Lenstronomy. The two approaches are similar in the choice of parameterization and philosophy, but have a few crucial differences: the description of the lensed galaxy light (pixellated vs. Sérsic+shapelets), the number of modeled bands (user choice vs. all available bands), priors (uniform vs. informative), and the PSF (initial guess based on stars vs. corrections based on the QSO images for S22). Despite the differences, Einstein radii and mass flattening agree well between the two studies, with a few exceptions that can be traced to inadequacy of the data or the PSF modeling. Other quantities such as convergence and image magnifications show differences that are generally larger than the estimated formal statistical uncertainties, but they are still encouraging in amplitude, considering the automated approach.

The differences primarily arise from two effects, as shown by our post-blind analysis. First, when the data quality is insufficient to constrain the slope of the power-law mass density profile, the S22 results are driven by the prior, while our current results tend to be more uncertain and constrained by the bounds of the uniform prior. We showed that if the same power-law slope is imposed, the agreement between the codes improves significantly. Second, the reconstruction of the PSF is a key factor in determining the slope and other parameters for data of sufficient quality. We illustrated this point by modeling two high-quality systems with the same PSF as S22, and we found a much better agreement. We conclude that for systems with high S/N, the PSF can be directly reconstructed from the data and the power-law slope can be stably inferred (this confirms the results obtained via detailed modeling by Wong et al. 2019; Birrer et al. 2019; Shajib et al. 2019b, 2022; Chen et al. 2019).

In terms of our automated modeling pipeline, we conclude that the models are a good starting point, but more work, particularly PSF corrections, and in some cases, better data, are needed to reach cosmography grade. Nonetheless, these automated modeling results provide important information about the lens systems, such as the approximate time delays, to help schedule the monitoring of the system. The lensing systems modeled with our automated pipeline are being followed up observation-ally to acquire the redshift and velocity dispersion of the lens, time delays, and environment properties, for example. We can use the modeling results that are presented here and models from lensing systems that are obtained with the automation code in the future, as input models for a more detailed modeling, in a fraction of the user time compared to conventional modeling.

Acknowledgements

We thank P. Schechter for data on the astrometric comparison and for helpful discussions. We further thank the anonymous referee for helpful comments. This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with programs HST-GO-15320 and HST-GO-15652. Support for the two programs was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. S.E., S.S., and S.H.S. thank the Max Planck Society for support through the Max Planck Research Group for SHS. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (LENSNOVA: grant agreement No 771776; COSMICLENS: grant agreement No 787886). This research is supported in part by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC-2094 - 390783311. T.S. and T.T. acknowledge support by the National Science Foundation through grant NSF-AST-1906976 and NSF-AST-1907396 "Collaborative Research: Toward a 1% measurement of the Hubble Constant with gravitational time delays". T.T. acknowledges support by the Packard Foundation through a Packard Research Fellowship. This work made use of NumPy (Oliphant 2015), SciPy (Jones et al. 2001), astropy (Astropy Collaboration 2018), matplotlib (Hunter 2007), and draw.io at https://www.draw.io.

Appendix A Final modeling results

Appendix Β Final model parameters

Appendix C Results and predictions for other lensing quantities

Appendix D Degeneracies between the mass parameters

Fig. D.1 Posterior distributions for the main mass parameters and external shear of J0659+1629. The correlation between the Einstein radius of the lens and the satellite is highlighted in red. The three shaded areas show the 1, 2, and 3σ credible regions. The one-dimensional histograms show the marginalized posterior distribution for the selected mass parameters, and the vertical lines mark the 1σ confidence intervals.

Appendix E Plots for the comparison with Lenstronomy models

Fig. E.1 Comparison of κ, γtot, μ, Fermat potential difference Δτ, and predicted time delay Δt between the GLEE and Lenstronomy models. Additionally, we compare κiso and Δτiso, where we assume isothermal instead of power-law mass profiles.

Fig. E.2 Change in κ, γtot, μ, Fermat potential difference Δτ, and predicted time delay Δt when using the Lenstronomy PSF with our GLEE automation code for J1606−2333 and J2100−4452 in the F160W (IR) band.

References

All Tables

All Figures

Fig. 1 Lensing system and regions for mask creation. Top: image cutout of J2100-4452 to illustrate its different components which are modeled with our automation code. Bottom: image cutout of J2100-4452 with regions provided by the user to create the masks. The region required for creating the arc mask, which covers the light of the lensed quasar and arcs, is plotted in red. For the lens mask, the user provides the region(s) of the surrounding luminous objects (green) that contain the pixels that are ignored by the model.

In the text

Fig. 2 Overview flowchart for single-band modeling of a strongly lensed quasar system with our new automated code. The modeling consists of two main steps: the lens mass distribution modeling with source or image positions, and the reconstruction of the source surface brightness (SSB; which is the unlensed light distribution of the quasar host galaxy). Each step is divided into smaller modeling steps. The details are described in flowcharts for each step and in the corresponding subsection, as referenced in the flowchart.

In the text

Fig. 3 Flowchart showing the decision tree for the lens light modeling. The modeling code obtains the best configuration of Sérsic and point-source profiles to fit the primary lens light. The quality check is made with a ${\chi }_{\text{red}}^2$ approach and a criterion for stability in the posterior values. In addition, the code also has the option of modeling the light of satellite lens galaxies.

In the text

Fig. 4 Flowchart showing the decision tree for modeling the lens mass distribution with source or image positions. After the user specifies a lens mass profile (PIEMD or SPEMD), the code models the lens mass distribution first by using the modeled source position, and then by using the positions of the multiple quasar images with the option to include an external shear component to the model. Both steps use simulated annealing until the respective χ2 is below 5. The image position modeling is completed after an approximately converged MCMC chain is achieved. After this, the code continues the lens mass modeling with the arc light modeling and reconstruction of the SSB distribution.

In the text

Fig. 5 Flowchart for arc light modeling and SSB reconstruction. In the first stage, the quasar amplitudes are lowered, and in the second stage, the error map is boosted in the bright quasar regions in order to model the light of the arc. The lens mass parameters are allowed to vary in this modeling step. In a final step, the code runs optimizing and sampling cycles until full convergence of the chain is achieved.

In the text

	Fig. 6 Overview flowchart for multiband modeling of a lensed quasar system with the multiband automated code. Details are described in flowcharts for each step and in the corresponding subsections, as referenced in the flowchart.
In the text

	Fig. 7 Color images for each of the nine lenses in our sample created with the three HST bands F160W (red), F 475X (blue), and F814W (green).
In the text

Fig. 8 Comparison between the position and structural parameters of the lens mass and light for eacht lensing system. Top: difference between the median mass and light centroid of our models, with Δx = xs − xm and Ay = ys − ym. Bottom left: comparison between the median mass and light position angle. Bottom right: comparison between the median mass and light-axis ratio. The dashed lines show lines without a centroid offset (top) and the lines of equal position angle and axis ratio (bottom left and bottom right).

In the text

	Fig. 9 Direct comparison of lens mass parameter and external shear values between our models and Lenstronomy models from S22.
In the text

	Fig. 10 Change in lens mass parameter and external shear values when using the Lenstronomy PSF with our GLEE automation code for J1606−2333 and J2100-4452.
In the text

	Fig. 11 Relative difference in quasar image positions between our results and S22 (squares), Luhtaru et al. (2021, circles), or Gaia (triangles) for systems that overlap in the multiple samples. The image positions of the outlier J1721+8842 are marked in orange (top) or are excluded (bottom).
In the text

Fig. D.1 Posterior distributions for the main mass parameters and external shear of J0659+1629. The correlation between the Einstein radius of the lens and the satellite is highlighted in red. The three shaded areas show the 1, 2, and 3σ credible regions. The one-dimensional histograms show the marginalized posterior distribution for the selected mass parameters, and the vertical lines mark the 1σ confidence intervals.

In the text

	Fig. E.1 Comparison of κ, γtot, μ, Fermat potential difference Δτ, and predicted time delay Δt between the GLEE and Lenstronomy models. Additionally, we compare κiso and Δτiso, where we assume isothermal instead of power-law mass profiles.
In the text

	Fig. E.2 Change in κ, γtot, μ, Fermat potential difference Δτ, and predicted time delay Δt when using the Lenstronomy PSF with our GLEE automation code for J1606−2333 and J2100−4452 in the F160W (IR) band.
In the text