$\gamma$ distribution becomes consistent within $\sim 0.6\sigma$. However, the
$\gamma$ distribution of the glee model with its own supersampled PSF still
differs by $2.6\sigma$ 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.

Figure 25:  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. Note
that the illustrated supersampled lenstronomy PSF was not used in the pre-
unblinding models and was reconstructed after the unblinding to perform
further tests. The original reconstructed PSF in lenstronomy in F160W was not
supersampled. The glee PSF FWHMs are smaller by $\sim$2–7% than the
lenstronomy ones.  Figure 26:  Deviations in the logarithmic slope $\gamma$ of
the power-law model with different PSF settings. In all panels, the dashed
distributions show the fiducial constraints – lenstronomy in green and glee in
purple. The lenstronomy fiducial constraint is from the highest BIC value
power-law model (Table 3), and the glee fiducial constraint is from the
“Power-law fiducial” setup (Table 2). First panel: When only the reconstructed
PSFs are interchanged to optimize the models, the constraint from one software
program moves toward the other’s fiducial constraint to be consistent within
$1.2-1.6\sigma$. 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 $\gamma$ becomes slightly more consistent between
the software programs within $1.0-1.5\sigma$. Therefore, we conclude that the
differences in the PSF itself significantly contributes to the discrepancy in
the power-law $\gamma$ 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 toward 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 supersampling
factor 3 reconstructed by lenstronomy, the resultant $\gamma$ constraint
agrees very well within $0.6\sigma$.

###  7.4 Requirements for computational resources

The entire lenstronomy modelling procedure required $\sim$ $4\times 10^{5}$
CPU hours including initial modelling trials, running full MCMC chains of the
adopted models, post-processing of posteriors, and post-unblinding tests. The
estimated usage for glee models are $\mathcal{O}(10^{5})$ CPU hours only to
run the final models and MCMC chains to produce the final unblinded result.
However, the total usage of CPU hours can be $\mathcal{O}(10^{6})$ CPU hours
including modelling trials, robustness tests, and reruns of chains to recover
lost progress due to numerical issues.

###  7.5 Discussion

The final unblinded time-delay predictions agree within $<1.2\sigma$ between
the two modelling teams. As a result, the inferred Hubble constants from the
two teams based on the observed time delays will be consistent within
$<1.2\sigma$. 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 toward the same direction for both
teams. We were already aware prior to the unblinding that the composite model
is atypical with respect to previous observations, e.g., the velocity
dispersion profile significantly decreases toward the center for the one from
the lenstronomy (Fig. 14), the one from the glee team has a very low inner
dark matter fraction. Although such discrepancy between composite and power-
law model predictions have been observed in the previously analyzed systems
(e.g., Suyu et al. 2014), this system WGD 2038$-$4008 demonstrates the largest
discrepancy to date out of the systems analyzed by H0LiCOW/TDCOSMO. However,
unlike the previously analyzed systems, the Einstein radius of this system
encompasses only the very central region of the very extended lens galaxy,
thus the imaging observables probe a different region of the elliptical
galaxy, i.e., at $\sim\theta_{\rm eff}/3$ instead of $\sim\theta_{\rm eff}$.
This discrepancy in the composite model illustrates that this model is not an
adequate description for the mass distribution at the central region of all
elliptical galaxies. Rather an appropriate combination of a mass-to-light
ratio gradient and a SMBH can be necessary to sufficiently describe the mass
distribution at the scales considered here. In such cases when the different
mass model families fit the imaging observables equally well, but lead to
different predictions in the Fermat potential and kinematics, the observed
kinematics is crucial to act as the differentiator between the mass model
families through appropriate weighting of the kinematics likelihood. In the
future, spatially resolved velocity dispersion from integral field spectra –
e.g., from Muse/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 a large room for
discrepancy 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 $\gamma$,
as this parameter is the most sensitive lens model parameter to the predicted
time delays and thus the inferred Hubble constant. The $\gamma$ distributions
between the modelling teams are discrepant at $4\sigma$. 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 $\gamma$ constraints with differences below 0.5%. The $\gamma$
distributions are not expected to be identical from the two modelling
softwares due to differences in the numeric implementation, e.g., the
likelihood computation region, 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 $\gamma$ corresponding to $\sim$1% on H0.
However, the PSF reconstruction method can lead to systematic differences as
large as $\sim 4$%. This difference is reconciled in the time-delay
predictions from the two teams after combining the lens model posteriors with
the kinematics data and the estimated external convergence for the particular
system WGD 2038$-$4008\. However, it is inconclusive if such differences can
be similarly reconciled for all lens systems in general, and we leave this
investigation with a larger lens sample for a future study. If such
differences cannot be reconciled for other lens systems and these differences
do not average out when a sample of lenses is considered, then these
differences would be non-negligible in the long run when a large sample of
lens systems are combined to infer the Hubble constant. As it is currently not
possible to evaluate the appropriateness of one reconstructed PSF over the
other, we recommend to marginalize over different realizations of the PSF
reconstruction and also over different reconstruction algorithms to avoid any
potential bias. Furthermore, a supersampled PSF in the IR band is also
recommended, as the drizzled pixel scale of 0$\aas@@fstack{\prime\prime}$08 in
the IR does not optimally sample the PSF.

##  8 Summary and Conclusion

In this study, two teams independently modelled the lens system WGD
2038$-$4008 from 3-band HST imaging using two different software programs:
lenstronomy and glee. Two families of models were specified as baseline models
before performing the modelling – one family with power-law ellipsoidal mass
distribution and the other family with a two-component mass distribution
separately accounting for the dark and baryonic components. The baseline
settings were pre-specified to allow fair comparison between the modelling
results. Individual teams were allowed to improve upon the baseline settings
as they deem appropriate, for example, choice of priors, numerical settings
specific to the software program being used. Both modelling procedures were
carried out blindly from each other’s results. The models were unblinded on
October 22, 2021 after internal review by scientists from the TDCOSMO
collaboration not directly involved with either modeling team, and no further
modifications in 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 investigate the observed systematic differences
between the model outputs from the two teams and identify that differences in
the reconstructed PSF is the dominant source of systematic differences. The
main results of this study are as follows:

  * •

The final predicted time delays from lenstronomy are: $\Delta t_{\rm
AB}=-5.0^{+0.2}_{-0.2}$ d, $\Delta t_{\rm AC}=-10.0^{+0.4}_{-0.3}$ d, and
$\Delta t_{\rm AD}=-24.2^{+1.0}_{-0.7}$ d; and the ones from glee are: $\Delta
t_{\rm AB}=-4.4^{+0.4}_{-0.5}$ d, $\Delta t_{\rm AC}=-9.4^{+0.7}_{-0.8}$ d,
and $\Delta t_{\rm AD}=-23.0^{+1.8}_{-2.4}$ d. These values assume the flat
$\Lambda$CDM cosmology with $H_{0}=70$ km s-1 Mpc-1, and $\Omega_{\rm m}=0.3$.
The negative value of $\Delta t_{\rm AB}$, for example, signifies that image B
lags image A. This system is currently being monitored under the COSMOGRAIL
program to measure the time delays (Eigenbrod et al. 2005). Once the time
delays are measured, the mutual agreement between the predicted values will
result in a mutually consistent inference of the Hubble constant.

  * •

The logarithmic slope $\gamma$ and the external shear $\gamma_{\rm ext}$ of
the power-law model deviate by 4$\sigma$ (3.9%) between the lenstronomy and
glee models. This discrepancy is predominantly created by the difference in
the reconstructed pixelized PSFs. If the same PSF is used by both modelling
softwares, then the resultant $\gamma$ and $\gamma_{\rm ext}$ distributions
agree within $\sim$0.6$\sigma$ ($\sim$0.5%), which is compatible with the
$\sim$1% precision goal in the Hubble constant measurement from a large sample
of $\sim$40 quad lenses (e.g., Shajib et al. 2018; Birrer & Treu 2021). The
particular method of weighting the uncertainty map to account for the PSF
uncertainty is non-dominant in this discrepancy.

  * •

Our composite model posteriors are not generally in good agreement with the
one from the power-law model, and the discrepancy is more prominent for
lenstronomy models due to adopting more stringent physical priors on the halo
mass. This inconsistency points to inadequacy of our composite model to
describe the mass distribution for this particular lens system WGD
2038$-$4008; WGD 2038$-$4008 is atypical compared to previously analyzed
TDCOSMO/H0LiCOW/STRIDES lenses in that the $\theta_{\rm E}/\theta_{\rm eff}$
is relatively smaller 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 7 systems
previously analyzed by the TDCOSMO/H0LiCOW/STRIDES 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. 2021). 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 supernovae will robustly measure the Hubble constant to $\lesssim 1$%
precision (e.g., Jee et al. 2016; Birrer & Treu 2021; Birrer et al. 2022).

###### Acknowledgements.

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

## References

  * Agnello et al. (2018) Agnello, A., Lin, H., Kuropatkin, N., et al. 2018, MNRAS, 479, 4345
  * Aiola et al. (2020) Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmology Astropart. Phys., 2020, 047–047
  * Astropy Collaboration (2013) Astropy Collaboration. 2013, A&A, 558, A33
  * Astropy Collaboration (2018) Astropy Collaboration. 2018, AJ, 156, 123
  * Barkana (1998) Barkana, R. 1998, ApJ, 502, 531
  * Barkana (1999) Barkana, R. 1999, ASCL, ascl:9910.003
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Bertin & Lombardi (2006) Bertin, G. & Lombardi, M. 2006, ApJ, 648, L17
  * Birrer (2021) Birrer, S. 2021, ApJ, 919, 38
  * Birrer & Amara (2018) Birrer, S. & Amara, A. 2018, Physics of the Dark Universe, 22, 189
  * Birrer et al. (2015) Birrer, S., Amara, A., & Refregier, A. 2015, ApJ, 813, 102
  * Birrer et al. (2016) Birrer, S., Amara, A., & Refregier, A. 2016, J. Cosmology Astropart. Phys., 8, 020
  * Birrer et al. (2022) Birrer, S., Dhawan, S., & Shajib, A. J. 2022, ApJ, 924, 2
  * Birrer et al. (2020) Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, A165
  * Birrer et al. (2021) Birrer, S., Shajib, A. J., Gilman, D., et al. 2021, JOSS, 6, 3283
  * Birrer & Treu (2019) Birrer, S. & Treu, T. 2019, MNRAS, 489, 2097–2103
  * Birrer & Treu (2021) Birrer, S. & Treu, T. 2021, A&A, 649, A61
  * Birrer et al. (2019) Birrer, S., Treu, T., Rusu, C. E., et al. 2019, MNRAS, 484, 4726
  * Blakeslee et al. (2021) Blakeslee, J. P., Jensen, J. B., Ma, C.-P., Milne, P. A., & Greene, J. E. 2021, ApJ, 911, 65
  * Blandford & Narayan (1992) Blandford, R. D. & Narayan, R. 1992, ARA&A, 30, 311
  * Bolton et al. (2006) Bolton, A. S., Burles, S., Koopmans, L. V. E., Treu, T., & Moustakas, L. A. 2006, ApJ, 638, 703
  * Bonfini & Graham (2016) Bonfini, P. & Graham, A. W. 2016, ApJ, 829, 81
  * Bonvin et al. (2017) Bonvin, V., Courbin, F., Suyu, S. H., et al. 2017, MNRAS, 465, 4914
  * Bradley et al. (2020) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, Zenodo
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Buckley-Geer et al. (2020) Buckley-Geer, E. J., Lin, H., Rusu, C. E., et al. 2020, MNRAS, 498, 3241
  * Cappellari (2002) Cappellari, M. 2002, MNRAS, 333, 400
  * Cappellari (2016) Cappellari, M. 2016, ARA&A, 54, 597
  * Chen et al. (2019) Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743
  * Chen et al. (2021) Chen, G. C.-F., Fassnacht, C. D., Suyu, S. H., et al. 2021, A&A, 652, A7
  * Chen et al. (2016) Chen, G. C.-F., Suyu, S. H., Wong, K. C., et al. 2016, MNRAS, 462, 3457
  * Diemer (2018) Diemer, B. 2018, ApJS, 239, 35
  * Diemer & Joyce (2019) Diemer, B. & Joyce, M. 2019, ApJ, 871, 168
  * Ding et al. (2021) Ding, X., Treu, T., Birrer, S., et al. 2021, MNRAS, 503, 1096
  * Dobler et al. (2015) Dobler, G., Fassnacht, C. D., Treu, T., et al. 2015, ApJ, 799, 168
  * Dullo (2019) Dullo, B. T. 2019, ApJ, 886, 80
  * Dutton et al. (2011) Dutton, A. A., Brewer, B. J., Marshall, P. J., et al. 2011, MNRAS, 417, 1621
  * Efstathiou (2021) Efstathiou, G. 2021, MNRAS, 505, 3866
  * Eigenbrod et al. (2005) Eigenbrod, A., Courbin, F., Vuissoz, C., et al. 2005, A&A, 436, 25
  * Ene et al. (2019) Ene, I., Ma, C.-P., McConnell, N. J., et al. 2019, 878, 57
  * Falco et al. (1985) Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Freedman (2021) Freedman, W. L. 2021, ApJ, 919, 16
  * Freedman et al. (2019) Freedman, W. L., Madore, B. F., Hatt, D., et al. 2019, ApJ, 882, 34
  * Freedman et al. (2020) Freedman, W. L., Madore, B. F., Hoyt, T., et al. 2020, ApJ, 891, 57
  * Gavazzi et al. (2007) Gavazzi, R., Treu, T., Rhodes, J. D., et al. 2007, ApJ, 667, 176
  * Gilman et al. (2020) Gilman, D., Birrer, S., & Treu, T. 2020, A&A, 642, A194
  * Golse & Kneib (2002) Golse, G. & Kneib, J.-P. 2002, A&A, 390, 821
  * Goodman & Weare (2010) Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65–80
  * Greene et al. (2013) Greene, Z. S., Suyu, S. H., Treu, T., et al. 2013, ApJ, 768, 39
  * Hilbert et al. (2009) Hilbert, S., Hartlap, J., White, S. D. M., & Schneider, P. 2009, A&A, 499, 31
  * Hoeting et al. (1999) Hoeting, J. A., Madigan, D., Raftery, A. E., & Volinsky, C. T. 1999, Statistical Science, 14, 382–417
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
  * Jee et al. (2016) Jee, I., Komatsu, E., Suyu, S. H., & Huterer, D. 2016, J. Cosmology Astropart. Phys., 4, 031
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., & Others. 2001, SciPy: Open source scientific tools for Python
  * Kennedy & Eberhart (1995) Kennedy, J. & Eberhart, R. 1995, in Proceedings of ICNN'95 - International Conference on Neural Networks (IEEE)
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt (IOS Press BV, Amsterdam, Netherlands), 87 – 90
  * Knox & Millea (2020) Knox, L. & Millea, M. 2020, Phys. Rev. D, 101, 043533
  * Kochanek (2020) Kochanek, C. S. 2020, MNRAS, 493, 1725–1735
  * Kourkchi et al. (2020) Kourkchi, E., Tully, R. B., Eftekharzadeh, S., et al. 2020, ApJ, 902, 145
  * Lejeune et al. (1998) Lejeune, T., Cuisinier, F., & Buser, R. 1998, A&AS, 130, 65
  * Lewis & Bridle (2002) Lewis, A. & Bridle, S. 2002, Phys. Rev. D, 66, 103511
  * Liao et al. (2015) Liao, K., Treu, T., Marshall, P., et al. 2015, ApJ, 800, 11
  * Madigan & Raftery (1994) Madigan, D. & Raftery, A. E. 1994, Journal of the American Statistical Association, 89, 1535
  * Mamon & Łokas (2005) Mamon, G. A. & Łokas, E. L. 2005, MNRAS, 363, 705
  * Mehrgan et al. (2019) Mehrgan, K., Thomas, J., Saglia, R., et al. 2019, ApJ, 887, 195
  * Melo et al. (2021) Melo, A., Motta, V., Godoy, N., et al. 2021, A&A, 656, A108
  * Merritt (1985a) Merritt, D. 1985a, MNRAS, 214, 25P
  * Merritt (1985b) Merritt, D. 1985b, AJ, 90, 1027
  * Millon et al. (2020) Millon, M., Galan, A., Courbin, F., et al. 2020, A&A, 639, A101
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Oliphant (2015) Oliphant, T. E. 2015, Guide to NumPy, 2nd edn. (USA: CreateSpace Independent Publishing Platform)
  * Osipkov (1979) Osipkov, L. P. 1979, Pisma v Astronomicheskii Zhurnal, 5, 77
  * Pesce et al. (2020) Pesce, D. W., Braatz, J. A., Reid, M. J., et al. 2020, ApJ, 891, L1
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A6
  * Poulin et al. (2019) Poulin, V., Smith, T. L., Karwal, T., & Kamionkowski, M. 2019, Physical Review Letters, 122, 221301
  * Refregier (2003) Refregier, A. 2003, MNRAS, 338, 35
  * Refsdal (1964) Refsdal, S. 1964, MNRAS, 128, 307
  * Riess et al. (2021) Riess, A. G., Yuan, W., Macri, L. M., et al. 2021, arXiv:2112.04510 [astro-ph] [arXiv:2112.04510]
  * Rusu et al. (2017) Rusu, C. E., Fassnacht, C. D., Sluse, D., et al. 2017, MNRAS, 467, 4220
  * Rusu et al. (2020) Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440
  * Schneider et al. (1992) Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses
  * Schneider & Sluse (2014) Schneider, P. & Sluse, D. 2014, A&A, 564, A103
  * Schwarz (1978) Schwarz, G. 1978, The Annals of Statistics, 6, 461–464
  * Scolnic et al. (2018) Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101
  * Shajib et al. (2020) Shajib, A. J., Birrer, S., Treu, T., et al. 2020, MNRAS, 494, 6072
  * Shajib et al. (2019) Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649
  * Shajib et al. (2018) Shajib, A. J., Treu, T., & Agnello, A. 2018, MNRAS, 473, 210
  * Shajib et al. (2021) Shajib, A. J., Treu, T., Birrer, S., & Sonnenfeld, A. 2021, MNRAS, 503, 2380
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, American Institute of Physics Conference Series, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405
  * Sonnenfeld et al. (2018) Sonnenfeld, A., Leauthaud, A., Auger, M. W., et al. 2018, MNRAS, 481, 164
  * Speagle (2019) Speagle, J. S. 2019, arXiv e-prints, arXiv:1904.02180
  * Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Nature, 435, 629
  * Suyu (2012) Suyu, S. H. 2012, MNRAS, 426, 868
  * Suyu et al. (2013) Suyu, S. H., Auger, M. W., Hilbert, S., et al. 2013, ApJ, 766, 70
  * Suyu et al. (2017) Suyu, S. H., Bonvin, V., Courbin, F., et al. 2017, MNRAS, 468, 2590
  * Suyu & Halkola (2010) Suyu, S. H. & Halkola, A. 2010, A&A, 524, A94
  * Suyu et al. (2012) Suyu, S. H., Hensel, S. W., McKean, J. P., et al. 2012, ApJ, 750, 10
  * Suyu et al. (2010) Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201
  * Suyu et al. (2006) Suyu, S. H., Marshall, P. J., Hobson, M. P., & Blandford, R. D. 2006, MNRAS, 371, 983
  * Suyu et al. (2014) Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35
  * Sérsic (1968) Sérsic, J. L. 1968, Atlas de Galaxias Australes
  * Treu et al. (2006) Treu, T., Koopmans, L. V., Bolton, A. S., Burles, S., & Moustakas, L. A. 2006, ApJ, 640, 662
  * Treu & Marshall (2016) Treu, T. & Marshall, P. J. 2016, A&A Rev., 24, 11
  * Van de Vyvere et al. (2021) Van de Vyvere, L., Gomer, M. R., Sluse, D., et al. 2021, arXiv:2112.03932 [astro-ph] [arXiv:2112.03932]
  * Waskom et al. (2014) Waskom, M., Botvinnik, O., Hobson, P., et al. 2014, seaborn: v0.5.0 (November 2014)
  * Wong et al. (2017) Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895
  * Wong et al. (2020) Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420
  * Yıldırım et al. (2021) Yıldırım, A., Suyu, S. H., Chen, G. C.-F., & Komatsu, E. 2021, arXiv:2109.14615 [astro-ph] [arXiv:2109.14615]
  * Yıldırım et al. (2020) Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, MNRAS, 493, 4783