© The Authors 2026

1. Introduction

In recent years, the Λ cold dark matter (ΛCDM) cosmological paradigm has been challenged by observations, as tensions on several cosmological parameters now reach a statistically significant level. The most significant of these tensions is about the Hubble constant H₀, which gives the current expansion rate of the Universe. While the most recent H₀ measurements using the cosmic microwave background (CMB) obtain H₀ = 67.4 ± 0.5 km s⁻¹Mpc⁻¹ (Planck Collaboration VI 2020), local measurements find higher values, (e.g. Verde et al. 2019; Di Valentino et al. 2025). In particular, the distance ladder using parallaxes, water megamasers, Cepheids, and type-Ia supernovae yields H₀ = 73.17 ± 0.86 km s⁻¹Mpc⁻¹ (Breuval et al. 2024), which is in a 5σ tension with the CMB value. Given this discrepancy, it is essential to develop alternative cosmological probes that are fully independent of those listed above and that are competitive in terms of precision and accuracy.

Strong gravitational lensing occurs when a massive and compact object (the lens) intersects the line-of-sight (LOS) to a more distant galaxy, forming multiple images of the same background source. The lensed images appear at positions on the sky that minimise the photons’ travel time to the observer, following Fermat’s principle. These optical paths differ in length, resulting in a time delay between the arrival times of the photon in the different multiple images. These time delays are proportional to the so-called time delay distance D_Δt (e.g. Refsdal 1964; Suyu et al. 2014), a quantity that is inversely proportional to H₀. The method, first proposed by Refsdal (1964), requires the time-delay measurement and modelling of the mass along the LOS up to the source redshift, including the main deflector while breaking the inherent mass-sheet degeneracy (MSD, e.g. Falco et al. 1985; Gorenstein et al. 1988; Kochanek 2002; Schneider & Sluse 2013, 2014; Blum & Teodori 2021).

The TDCOSMO collaboration (e.g. Millon et al. 2020c; Birrer et al. 2020; Shajib et al. 2023) is working towards a 1% precision measurement of H₀ using time-delay cosmography (TDC) with strongly lensed quasars, one of the few cosmological probes that are truly independent of the local distance ladder. The initial set of H₀-measurements by our TDCOSMO collaboration (e.g. Suyu et al. 2010, 2017; Bonvin et al. 2017; Millon et al. 2020c; Wong et al. 2020) used empirically motivated parametric models, which is an ‘assertive’ (Treu et al. 2022) model assumption consistent with several independent observations from stellar dynamics (Cappellari 2016) or galaxy simulations (Navarro et al. 1997), and also from joint lensing-dynamics studies (Shajib et al. 2021). However, this assertive model assumption implicitly breaks the MSD, and thus risks underestimating the uncertainty beyond the sole constraining power of the data. Therefore, TDCOSMO-IV (Birrer et al. 2020) incorporates one additional degree of freedom in the mass model that is maximally degenerate with the MSD and constrains this additional degree of freedom using stellar kinematics, which breaks the MSD. As expected, the precision of the measurement decreases, but the value remains compatible with previous ones. Recently, spatially resolved dynamics of the lensing galaxy yielded from integral field unit (IFU) data have been shown to significantly improve the constraints on the mass profile and the H₀-measurement on individual lenses (Shajib et al. 2018; Yıldırım et al. 2020; Shajib et al. 2023).

The latest milestone measurement of the collaboration (TDCOSMO Collaboration 2025, hereafter TDCOSMO25) used new IFU data from the JWST, Keck, and VLT telescopes with the improved methodology of Knabel et al. (2025). The hierarchical inference, jointly modelling all TDCOSMO lenses together with a subset of the Sloan Lenses ACS and Strong Lenses in the Legacy Survey samples, in combination with information on Ω_m from independent datasets such as the Pantheon Supernovae catalog (Brout et al. 2022) or DESI’s baryonic acoustic oscillation measurement (Abdul Karim et al. 2025) on the TDCOSMO lenses, yields respectively $H_0=74.3_{-3.7}^{+3.1}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ and $H_0=74.8_{-3.4}^{+3.5}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$. This milestone marks significant improvements in both accuracy and precision. One of the main ways to further improve precision is to increase the number of time-delay lenses in the sample (Birrer & Treu 2021).

As they are more magnified, quadruply lensed quasars are easier to find serendipitously than doubles. Hence, lensed quasar datasets were initially dominated by them since the image multiplicity and the three independent delays also give more constraints on the mass modelling of the main deflector. However, these systems are rare in the sky, and one of the current limiting factors of the measurements is the sample size. As wide-field surveys emerged with deeper imaging, the discovery selection bias was attenuated (Lemon et al. 2020, 2023; Dux et al. 2024; Lemon et al. 2024) and doubles now represent 3/4 of known lensed quasars. They are generally easier to monitor from the ground because of their wider separation, and therefore have the potential to greatly reduce statistical noise. The extension of the TDCOSMO sample to doubles will accelerate progress towards the ultimate goal of 1% precision (Treu et al. 2022).

Furthermore, quads generally require more exquisite alignment of deflector and source, and higher ellipticity/shear than doubles. Therefore, restricting the TDCOSMO sample to quads may expose the H₀ measurement to unknown selection biases (Sonnenfeld et al. 2023). Analysing a sample of doubly lensed quasars offers a way to improve statistical precision and unveil and mitigate any potential biases that may be present.

The present paper is the first TDCOSMO paper to analyse a doubly imaged quasar¹, namely HE 1104−1805 discovered by Wisotzki et al. (1993) with a source at redshift z_s = 2.32 and deflector at redshift z_d = 0.729 (Wisotzki et al. 1993; Courbin et al. 1998; Lopez et al. 1999; Lidman et al. 2000). Following the collaboration protocol, the measurement was performed blindly to prevent experimenter bias. The mean value of the posterior distribution function of the Hubble Constant was kept hidden from the investigators and unblinded on July 8th 2025, during a collaboration teleconference. The paper was mostly written before unblinding, with placeholder figures and tables, and completed after unblinding. The unblinded results were published without modification.

The paper outline is as follows. Sect. 2 introduces the basics of TDC. The time-delay measurement is presented in Sect. 3 before turning to the spatially resolved kinematics of the lens with the ESO-VLT MUSE instrument in Sect. 4.1. Sect. 5 describes the mass contribution of galaxies along the LOS, and Sect. 6 shows the different models we developed to measure the time delay distance to HE 1104−1805. Sect. 7 gives our cosmological inference and H₀ value after unblinding of the whole analysis chain. Finally, Sect. 8 provides the limitations of our work, and we conclude in Sect. 8.4.

2. Time-delay cosmography formalism

The determination of H₀ using the time delay observable in strong lens systems was first suggested by Refsdal (1964) and further developed by Schneider et al. (1992). In this section, we review the main points of this formalism, highlighting how the redshifts of the source (z_s) and the lens (z_d), the time delay between each image (Δt_AB), the resolved kinematics, the LOS, the mass model of the lens and the position of the source are used to measure H₀².

2.1. Lensing formalism

The time delay Δt_AB between two images A and B of a given strong lens system is given byinfer

$$\begin{array}{cc}\hfill \mathrm{\Delta }{t}_{\mathrm{AB}}& =\frac{D_{\mathrm{\Delta }t}}{c}(\tau ({\mathit{\theta }}_{\mathbf{A}},\mathit{\beta })-\tau ({\mathit{\theta }}_{\mathbf{B}},\mathit{\beta })),\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \text{with}{D}_{\mathrm{\Delta }t}& \equiv (1+z_{\mathrm{d}})\frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{d}\mathrm{s}}},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \text{and}\tau (\mathit{\theta },\mathit{\beta })& =\frac{1}{2}{(\mathit{\theta }-\mathit{\beta })}^2-\varphi \left(\mathit{\theta }\right),\hfill \end{array}$$(3)

where θ is the apparent angular position of an image and β the position of the source, D_Δt is the so-called ‘time-delay distance’ defined as a function of the lens redshift z_d, and the angular distances D_d, D_s and D_ds, respectively, between the observer and the lens, the observer and the source, and the lens and the source. Finally, the Fermat potential τ(θ, β) depends on the lens potential ϕ(θ), which can be determined by modelling the mass distribution of the lens.

We can generalise Eq. (3) in the case where multiple lenses lie at different redshifts, creating P different lens planes. The time delay then becomes (e.g. Schneider et al. 1992)

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{AB}}={\displaystyle \sum _{i=1}^{\mathrm{P}}}\frac{D_{\mathrm{\Delta }t,i,i+1}}{c}[\frac{{({\mathit{\theta }}_{\mathrm{A},i}-{\mathit{\theta }}_{\mathrm{A},i+1})}^2}{2}-\frac{{\mathbf{(}\mathit{\theta }}_{\mathrm{B},i}-{\mathit{\theta }}_{\mathrm{B},i+1}{)}^2}{2}\\ \hfill -{\zeta }_{i,i+1}({\varphi }_i({\mathit{\theta }}_{\mathrm{A},\mathrm{i}})+{\varphi }_i({\mathit{\theta }}_{\mathrm{B},\mathrm{i}}))],\end{array}$$(4)

$$\begin{array}{c}\hfill D_{\mathrm{\Delta }t,i,i+1}\equiv (1+z_i)\frac{D_i{D}_{i+1}}{D_{ii+1}},\end{array}$$(5)

$$\begin{array}{c}\hfill {\zeta }_{i,i+1}\equiv \frac{D_{i,i+1}{D}_{\mathrm{s}}}{D_{i+1}{D}_{\mathrm{is}}},\end{array}$$(6)

By defining D_Δt, eff ≡ D_{Δt, l, s} as the time-delay distance between the main lens plane and the source plane, the time delay can be rewritten

$$\begin{array}{c}\hfill \mathrm{\Delta }{t}_{\mathrm{AB}}=\frac{{D_{\mathrm{\Delta }t}}^{\mathrm{eff}}}{c}\mathrm{\Delta }{\varphi }_{\mathrm{AB}}^{\mathrm{eff}}.\end{array}$$(7)

With the effective Fermat potential

$$\begin{array}{c}\hfill \mathrm{\Delta }{\varphi }^{\mathrm{eff}}(\theta )={\displaystyle \sum _{i=1}^{\mathrm{P}}}\frac{1+z_i}{1+z_{\mathrm{d}}}\frac{D_i{D}_{i+1}{D}_{\mathrm{ds}}}{D_{\mathrm{d}}{D}_{\mathrm{s}}{D}_{i,i+1}}[\frac{{({\theta }_i-{\theta }_{i+1})}^2}{2}-{\zeta }_{i,i+1}{\psi }_i({\theta }_i)].\end{array}$$(8)

As pointed out by Falco et al. (1985), the observables of a strong lensing system, such as the image position and the flux ratios, are invariant with respect to the mass-sheet transform (MST). This transformation of the convergence κ(θ) and the source plane coordinates β can be written as

$$\begin{array}{cc}\hfill {\kappa }_{\lambda }\left(\theta \right)& =\lambda \kappa \left(\theta \right)+(1-\lambda ),\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\mathit{\beta }}_{\mathit{\lambda }}& =\lambda \mathit{\beta }.\hfill \end{array}$$(10)

This degeneracy between the mass of the lens and the size of the source is called the mass-sheet degeneracy (MSD) and implies a degeneracy between H₀ and the time delay.

Although a pure mass sheet extending up to infinity in the lens plane is unphysical, the MSD is the manifestation of two physical effects: (1) the lensing effect of all the galaxies along the LOS other than the main lens. We refer to this effect as the external MSD, and it can be treated as equivalent to a mass sheet with a surface mass density κ_ext in units of the critical density for lensing; (2) the internal mass sheet degeneracy, which results from a small mismatch between the true mass profile of the lens and its model.

The total internal and external MSD parameters defined in Eq. (9) can therefore be decomposed into two terms

$$\begin{array}{c}\hfill \lambda =(1-{\kappa }_{\mathrm{ext}}){\lambda }_{\mathrm{int}},\end{array}$$(11)

where λ_int is the internal MSD parameter.

The external part of the MSD can be statistically derived if κ_ext is estimated from the mass of the galaxies along the LOS. This can be done either by counting the galaxies in a large aperture and comparing the density of the LOS to numerical simulations (see e.g. Wells et al. 2023, for detail) or with weak lensing (Tihhonova et al. 2018, 2020). The time-delay distance D_Δt can then be corrected D_Δt  = D_Δt′/(1−κ_ext), where D_Δt′ is the time-delay distance from mass modelling without accounting for κ_ext. The detailed computation of κ_ext is presented in Section 5.4.

Following previous TDCOSMO single-lens analyses (Wong et al. 2017; Rusu et al. 2020; Shajib et al. 2020; Wong et al. 2024), we assumed no internal MST (i.e. λ_int = 1). The next collaboration milestone that will hierarchically model all TDCOSMO lenses will release this parameter.

2.2. Kinematic modelling

The kinematic observable is the luminosity-weighted line-of-sight projected stellar velocity dispersion, denoted as σ^LOS. It is measured by targeting stellar absorption lines and quantifying their width with high-resolution spectra.

As in TDCOSMO25, we used the spherically aligned Jeans anisotropic modelling (JAM) method. This method builds on the orbital distribution f(x, v) of position and velocity of the stars in 3D motion in the galactic potential Φ. This potential is described by the steady-state collisionless Boltzmann equation (Binney & Tremaine 1987, hereafter, BT, Eqs. 4–13b)

$$\begin{array}{c}\hfill {\displaystyle \sum _{i=1}^3}(v_i\frac{\partial f}{\partial x_i}-\frac{\partial \mathrm{\Phi }}{\partial x_i}\frac{\partial f}{\partial v_i})=0.\end{array}$$(12)

We then obtain a single spherical Jeans equation (BT Eqs. 4–54)

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left[{\rho }_{\ast }(r){\sigma }_r^2(r)\right]}{\mathrm{d}r}+\frac{2{\beta }_{\mathrm{ani}}(r){\rho }_{\ast }(r){\sigma }_r^2(r)}{r}=-{\rho }_{\ast }\frac{\mathrm{d}\mathrm{\Phi }(r)}{\mathrm{d}r},\end{array}$$(13)

with ρ_*(r) the stellar density distribution and β_ani the orbital anisotropy, defined as a function of the ratio of tangential over radial velocity dispersion components, σ_θ and σ_r

$$\begin{array}{c}\hfill {\beta }_{\mathrm{ani}}\equiv 1-\frac{{\sigma }_{\theta }^2}{{\sigma }_r^2}.\end{array}$$(14)

Given σ^LOS, the values of σ_r and σ_θ are degenerate, leaving β_ani nearly unconstrained.

This introduces a degeneracy between σ^LOS and the 3D mass profile known as the mass-anisotropy degeneracy (MAD) (Binney & Mamon 1982; Merritt 1985). As demonstrated in Appendix A.2, the lens galaxy of HE 1104−1805 is a slow rotator. We hence followed TDCOSMO Collaboration (2025), and use the measurements of Cappellari et al. (2007) on slow rotating galaxies to apply the uniform prior $0.87\le \frac{{\sigma }_{\theta }^2}{{\sigma }_r^2}\le 1.12$.

A solution of Eq. (13) is given by

$$\begin{array}{c}\hfill {\sigma }_r^2(r)=\frac{G}{{\rho }_{\ast }(r)}{\displaystyle {\int }_r^{\infty }}\frac{M(s){\rho }_{\ast }(s)}{s^2}{\mathcal{J}}_{\beta }(r,s)\mathrm{d}s,\end{array}$$(15)

with M(r) the mass enclosed within the radius r and

$$\begin{array}{c}\hfill {\mathcal{J}}_{\beta }(r,s)=exp\left[{\displaystyle {\int }_r^s}2{\beta }_{\mathrm{ani}}(r\prime )\frac{\mathrm{d}r\prime }{r\prime }\right].\end{array}$$(16)

We then get the modelled velocity dispersion along the LOS, ${{\sigma }^{\mathrm{LOS}}}_{\mathrm{model}}$ with (BT, Eqs. 4–60)

$$\begin{array}{c}\hfill {\left({{\sigma }^{\mathrm{LOS}}}_{\mathrm{model}}\right)}^2=\frac{2}{{\mathrm{\Sigma }}_{\ast }(R)}{\displaystyle {\int }_R^{\infty }}[1-{\beta }_{\mathrm{ani}}(r)\frac{R^2}{r^2}]\frac{{\rho }_{\ast }(r){\sigma }_r^2(r)}{\sqrt{r^2-R^2}}r\mathrm{d}r,\end{array}$$(17)

with R the projected radius and Σ_* the enclosed surface stellar density, which can be constrained from the luminosity profile of the lens I(R) assuming, for instance, a constant mass-to-light ratio ΥΣ_*(R) = ΥI(R).

To compare this with the observed velocity dispersion along the LOS, σ^LOS, we weigh ${{\sigma }^{\mathrm{LOS}}}_{\mathrm{model}}$ with the point spread function (PSF, 𝒫) convolved light profile of the lens

$$\begin{array}{c}\hfill {\left({{\sigma }^{\mathrm{LOS}}}_{\mathrm{sph}}\right)}^2=\frac{{\int }_{\mathrm{A}}I(R){\left({{\sigma }^{\mathrm{LOS}}}_{\mathrm{model}}\right)}^2\ast \mathcal{P}\mathrm{d}A}{{\int }_{A}I(R)\ast \mathcal{P}\mathrm{d}A}.\end{array}$$(18)

Following the methodology of TDCOSMO25, we did not account for the inclination angle of the lens galaxy during the JAM modelling. To prevent bias raised by Huang et al. (2025), we adopted the axisymmetric correction factor${\sigma }_{\mathrm{axi}}^{\mathrm{ap}}/{\sigma }_{\mathrm{sph}}^{\mathrm{ap}}$ proposed by these authors. This factor was computed using the shape distribution from Li et al. (2018) and assuming random inclination. We then bin this factor in the same fashion as the IFU data (see Section 4.2 for more details) to obtain a correction factor for each bin, respectively: 0.992 ± 0.004, 0.977 ± 0.001, and 0.993 ± 0.004.

The final σ^LOS to be compared with data is then given by

$$\begin{array}{c}\hfill {\sigma }^{\mathrm{LOS}}={{\sigma }^{\mathrm{LOS}}}_{\mathrm{sph}}·\frac{{\sigma }_{\mathrm{axi}}^{\mathrm{ap}}}{{\sigma }_{\mathrm{sph}}^{\mathrm{ap}}}\end{array}$$(19)

The prediction of the observed LOS velocity dispersion from any model, irrespective of the approach, can be decomposed into a cosmology-dependent ratio of angular diameter distances and a cosmology-independent part following Birrer et al. (2016, 2019)

$$\begin{array}{cc}\hfill {\left({\sigma }^{\mathrm{LOS}}\right)}^2& =(1-{\kappa }_{\mathrm{ext}})\frac{D_{\mathrm{s}}}{D_{\mathrm{d}\mathrm{s}}}\frac{c^2}{J({\xi }_{\mathrm{lens}},{\xi }_{\mathrm{light}},{\beta }_{\mathrm{ani}})},\hfill \end{array}$$(20)

$$\begin{array}{cc}& =\frac{(1-{\kappa }_{\mathrm{ext}})}{1+z_{\mathrm{d}}}\frac{D_{\mathrm{\Delta }\mathrm{t}}}{D_{\mathrm{d}}}\frac{c^2}{J({\xi }_{\mathrm{lens}},{\xi }_{\mathrm{light}},{\beta }_{\mathrm{ani}})},\hfill \end{array}$$(21)

where the dimensionless quantity J relies on the deflector model parameters (ξ_lens and ξ_light) and β_ani. Moreover, J considers observational conditions and luminosity-weighting within the dispersion measurement aperture, as previously demonstrated in studies such as Binney & Mamon (1982) and Treu & Koopmans (2004).

As the angular diameter distances are sensitive to Ω_m only to the second order, they are primarily sensitive to H₀. Large samples of lenses and combination with probes independent from Time-delay Cosmography are required to constrain other parameters (e.g. Linder 2011; Shajib et al. 2025).

3. Time-delay measurement

3.1. Optical photometric monitoring

As one of the first discovered strongly lensed quasars, multiple attempts at measuring the time delay of HE 1104−1805 with different observation campaigns have been published in the literature. First, Wisotzki et al. (1998) published 18 spectroscopic observations with the 3.6 m ESO telescope spanning over 5 years and estimated that the time delay between images A and B, Δt_AB, should be between 109 and 292 days. Reanalysis of these data by Gil-Merino et al. (2002) and Pelt et al. (2002) showed that this sampling is insufficient for a precise delay. With a better sampling of the R-band photometry between 1997 and 2006 with the OGLE and SMARTS programs³, Poindexter et al. (2007) estimated Δt_AB = 152 ± 3 days using a Legendre-polynomial fitting technique. Later, Morgan et al. (2008) re-estimated to Δt_AB = 162.2 ± 6.1 days by jointly fitting the time delay and accretion disk size of the quasar by exploiting the microlensing signal contained in the light curves. However, the method employed depends on physical assumptions about the quasar accretion disk profile and population of microlenses that could bias the measurement.

In this work, we used three new unpublished datasets and applied the established PyCS spline-fitting method to re-estimate the time delay of HE 1104−1805 in a data-driven way without relying on any prior assumption about physical properties of the system. The first dataset is provided by the SMARTS telescope, which monitored it from 2003 to 2016 with a ∼ weekly cadence. Next, the COSMOGRAIL program (Courbin et al. 2005) observed it every ∼4 days with ECAM and C2 at the Leonhard Euler 1.2 m Swiss Telescope from 2013 to 2018 and during the 2017 season with WFI at the MPG/ESO 2.2 m Telescope every ∼2 days. These three datasets are called SMARTS, ECAM, and WFI, respectively.

For each dataset, the quasar images’ photometry was measured with the same methodology as Millon et al. (2020a) that we briefly summarise here. After bias, flat, and sky corrections, each exposure’s PSF is modelled as a Moffat profile fit combined with regularised pixel adjustments of five stars with similar luminosity as the lens in the field. The MCS deconvolution algorithm (Magain et al. 1998) is then applied to the lensed images in order to deblend the flux of an image from its counterpart, the lensing galaxy, and the arc.

Finally, a deconvolution of reference stars is carried out using the previously built PSF. Stars displaying the most photometric stability are chosen to calculate a median photometric normalisation coefficient for each exposure. This step helps to correct image-to-image systematics resulting from PSF variations across epochs. Fig. 1 shows a stack of the exposures taken with the WFI, with the stars used to compute the PSF and the normalisation of each epoch.

Fig. 1. Left panel: Stack of 280 WFI exposures of the field of view used for light curve extraction. Stars nominated PSF1 to PSF5 are used to model the PSF of each exposure, and stars designated as N1 to N5 are used to normalise the flux between each exposure. Right panel: HE 1104−1805 imaging in the filter F160W band using HST WFC3. The main lensing galaxy is denoted as G, whereas the main perturbers within 5″ of the lens are numbered from P1 to P9.

The resulting light curves are shown on the top panel of Fig. 2. The three monitoring datasets presented are merged into a single light curve by fitting for possible magnitude and flux shifts between instruments on the overlapping parts of the light curves. The resulting combined dataset is referred to as ‘ECAM+SMARTS+WFI’.

Fig. 2. Top panel: HE 1104−1805 R-band light curve obtained by joining three datasets: SMARTS (yellow and blue squares), ECAM (purple and green circles), and WFI (brown and pink triangles). For clarity, the light curve of B was shifted by –1.05 mag. The inset zooms on the WFI dataset span, showcasing its superior sampling to the ECAM dataset. Bottom panel: Example of a simultaneous fit of an intrinsic light curve with η = 45 days and nml = 15. The time shift obtained by this fit gives a point estimate of the system’s time delay. By repeating such measurements with randomised starting points 800 times, we obtained the time-delay measurement of this (η, nml) configuration.

3.2. Curve shifting procedure

The time delay was computed using the method described extensively by Millon et al. (2020b), relying on free knot spline fitting implemented in the PyCS python package⁴ (Millon et al. 2020a). Free-knot splines are piece-wise polynomials where the mean number of days between two consecutive knots is controlled by the parameter η, which filters out variation with shorter time scales. A single free-knot spline is used to model the intrinsic variability of the quasar, which is identical in both light curves; meanwhile, an additional spline is fitted on each light curve to model the extrinsic variability caused mainly by microlensing. The number of degrees of freedom is controlled by the hyperparameter η for the intrinsic spline. For the extrinsic splines n_ml sets the number of knots throughout the whole curve. We fitted the intrinsic, image A and B extrinsic components, and the time delay between the light curves. A visual first guess of the time delay is used as a starting point of the fit ⁵. The bottom panel of Fig. 2 shows an example of such a fit, giving us a point estimate of the time delay.

We repeated the fit 800 times with starting points randomised uniformly in a range of 10 days around the first guess value, we took the median value as the time delay measurement for a given (η, n_ml).

To assess the uncertainty of this measurement in a data-driven way, a set of mock light curves is generated with built-in time delays uniformly drawn in an interval of ±10 days around the median value obtained previously. we used the generative noise model described by Millon et al. (2020b) to ensure that the simulated light curves have the same noise properties as the real data. The distribution of the difference between the measurement and the truth for each of the mocks gives an estimation of the systematic uncertainty (mean of the distribution) and the statistical uncertainty (standard deviation) of the time delay estimation. The final uncertainty of the time delay estimation is the sum in quadrature of these two sources of uncertainty.

The choice of values for η and n_ml is guided by the fact that quasars vary significantly faster than microlensing. Moreover, whether a given feature is mathematically attributed to the intrinsic or extrinsic splines depends on the flexibility balance between the splines. We therefore reproduced the same measurement for different values of these two hyperparameters. Following Millon et al. (2020b), we took η ∈ [25, 35, 45, 55] days and n_ml ∈ [9, 11, 13, 15] days when fitting the SMARTS, ECAM, and the merged datasets. In contrast, the microlensing in the much shorter WFI dataset is fitted with a spline having either one or two knots (n_ml ∈ [1, 2]) since the microlensing variability is expected to alter the light curve on time scales longer than a year.

The left panel of Fig. 3 shows the time-delay measurements obtained for each hyper-parameter configuration (η, n_ml) considered. The measurements obtained are consistent with each other with variable precision. We note that for a given n_ml, fits with η = 55 are slightly lower, which shows that the intrinsic fit is marginally too rigid. We are therefore conservative by keeping this value and should not use higher values for η. We gradually marginalised the most precise measurement with the others in order of increasing precision until the tension with the rest of the sample is lower than the threshold τ. Following the methodology of Millon et al. (2020b), we used τ = 0.5 (i.e. no individual measurement is further than 0.5σ away from the final measurement) to have a conservative uncertainty estimation without being biased by outlier measurements.

Fig. 3. Left panel: Measurements of the time delay for different configurations of (η, nml) with the merged ECAM+SMART+WFI dataset. Right panel: Measurement of ΔtAB for the different datasets. Since the datasets overlap, they cannot be considered fully independent. We therefore use our estimate on the joint dataset ECAM+SMARTS+WFI (green data point) as our final estimate. In both panels, the ‘combined τ = 0.5’ value was obtained by marginalising over the measurements listed above it.

Since the three variability components are degenerate and not physically informed, the resulting splines cannot be interpreted as the true source and microlensing variations. However, a visual inspection shows that the magnitude scales of the extrinsic variations are below 0.3 mag over several years, consistent with the expected microlensing behaviour (e.g. Mortonson et al. 2005; Schneider et al. 2006; Mosquera & Kochanek 2011). We, however, note that the number of features in the extrinsic splines is higher than in most systems, which hints towards a more complex microlensing behaviour, as observed by Schechter et al. (2003) in previously published light curves, making the extrinsic fit sensitive to the hyperparameter changes.

The same experiment is carried out separately for each dataset, and the right panel of Fig. 3 shows the time delay measured in each case, the combined value, and the measurement on the merged ECAM+SMARTS+WFI dataset. The marginalisation over the three datasets would be relevant if each dataset were subject to different sources of systematic errors in the photometric measurement. However, the photometry was computed with the same method for each dataset, and no discrepancies were observed between time delays measured with ECAM and WFI data on other lens systems (e.g. Millon et al. 2020b,a). As shown by the lower panel of Fig. 2, the extrinsic variations are complex during the 2017 season. The particular case of WFI is challenging because it covers only 253 days: given the long time delay, the light curves overlap on less than 100 days. Even though distinct features are visible and well-sampled in the WFI data, the complexity of the microlensing noticed previously cannot be adequately modelled with a single season. The time delay with the WFI dataset is, thus, not robust enough for the marginalised time delay value to be trusted. Therefore, we preferred to use the ECAM+SMARTS+WFI dataset, which allowed us to use the entire period of the light curve to constrain the microlensing behaviour as well as to take advantage of the high cadence and high signal-to-noise (S/N) display of the WFI dataset.

We did not take into account a possible microlensing time delay (Tie & Kochanek 2017) because it impacts the time delay only by a day in the worst cases (Bonvin et al. 2018), which is negligible given the length of the time delay. Therefore, we use the time-delay $\mathrm{\Delta }{t}_{\mathrm{AB}}=176.3_{-10.3}^{+11.4}$ days, i.e. a 5.8% precision for our cosmography analysis.

4. Spatially resolved kinematics of the lens

4.1. Integral field unit spectroscopy with MUSE

A 660-second IFU observation of the system was carried out with the multi-unit spectroscopic explorer (MUSE) instrument (Bacon et al. 2010) on March 13th 2019, in wide-field mode with adaptive optics (PI: A. Agnello, program 0102.A-0600(A) and 106.215F.001). The observation wavelength ranges from 4700.03 Å to 9351.28 Å with a spectral resolution of R ∼ 2500 corresponding to an instrumental dispersion of ${\sigma }_{\mathrm{inst}}=120\mathrm{km}{\mathrm{s}}^{-1}$. The spatial resolution is 0.2^″ pixel⁻¹ and the mean seeing throughout the wavelength after the adaptive optic correction is ∼0.57^″. The data reduction has been carried out as in Sluse et al. (2019), using the MUSE reduction pipeline.

This range includes the lens redshifted lines CaH [MYAMP] K lines at 3933 Å and 3968 Å along with the G-band at 4304 Å, while 5806 Å to 5965 Å wavelengths were cut out due to the notch filter. The data cube integrated across wavelength is shown in the top panel of Fig. 4. To ease spectral deblending of the quasar images from the lens galaxy, we subtracted the point source light at each wavelength of the data cube. To do so, at each wavelength, the point sources and lens galaxy light are fitted simultaneously with a Moffat and Sérsic profile with relative positions fixed to the ones measured on the HST data. For each wavelength frame, the Moffat part of the fit is then subtracted from the data to obtain the data cube displayed in the inset of Fig. 4.

Fig. 4. MUSE data cube summed over all wavelengths. The inset displays the point-source subtracted cube, revealing the lens galaxy light. The apertures used to extract individual spectra of the quasar images, lens centre, and perturbers shown in Fig. B.1 are represented by orange circles.

The spectrum of the centre of the lens galaxy obtained from this subtracted cube is shown in Fig. B.1 and displays observation with an S/N of 30 Å⁻¹, which is used to measure σ^LOS. We also retrieved new high-quality quasar spectra thanks to a small aperture integration on images A and B, and the spectra of some perturbers annotated on Fig. 1 (See Appendix B.1). Finally, we extracted the spectra of the two brightest perturbers (P5 and P6) with a sufficiently high S/N (∼ 11 Å⁻¹ and ∼ 7 Å⁻¹ respectively) to measure their redshifts.

4.2. Stellar template fitting

As shown by the left panel of Fig. 5, we masked the pixels from the lens light cube further than 1.8^″ from the centre of the galaxy, as well as an area of 0.8^″ around the position of images A and B. In order to maximise the number of constraints on the radial velocity profile, the frame is then divided into three concentric rings of 0.53^″ width with S/N from 14 Å⁻¹ to 30 Å⁻¹. The integrated spectra of each bin are then fitted using the penalised pixel fitting (pPXF) method⁶ (Cappellari 2017, 2023). This method models the galaxy spectrum with a weighted linear combination of stellar spectra, broadened by a convolution with the galaxy line-of-sight velocity distribution. Additive Legendre polynomials are included in the model to improve the robustness to template mismatch and dust reddening. Furthermore, a multiplicative polynomial accounts for spectral flux calibration mismatches between data and templates. Any contaminating signal from residuals of the quasar images blended with the galaxy is accounted for by including a scaled quasar spectrum presented in Fig. B.1.

Fig. 5. Velocity dispersion point-estimate in the three radial bins of HE 1104−1805’s lens galaxy. Left panel: Mean of the MUSE data cube after PSF subtraction and masking with overlaid bin numbers and contours. Right panel: Example of a pPXF fit of the integrated spectra of each bin, where data masked for the fit are marked in grey.

As shown in Fig. B.1, the most constraining features of the lens galaxy spectrum for the σ^LOS measurement are the CaH, CaK, and G-band absorption lines identified respectively at 6800 Å, 6828 Å and 7395 Å in the observed frame. To focus the fit on these lines and avoid regions of the spectrum that could be dominated by the C IV and C III emission of the quasars, we restricted the wavelength range to [6220, 7780] Å (in the rest-frame [3600, 4400] Å).

Following Knabel et al. (2025), we used three stellar template libraries: Indo-US (Valdes et al. 2004), MILES (Sánchez-Blázquez et al. 2006; Falcón-Barroso et al. 2011), and X-shooter spectral library (Verro et al. 2022), cleaned from incomplete or low-quality spectra and convolved them down to the same FWHM as the observation’s Line Spread Function in the galaxy’s rest frame, estimated at 1.46 Å with the empirical relation given by Bacon et al. (2023).

An exploratory set of fits allowed us to determine that an additive polynomial of degree 6 combined with a multiplicative polynomial of degree 0 yields the most robust measurements across the different template libraries.

An example of the σ^LOS estimate in each bin given by a fit with this fiducial setup is shown in Fig. 5.

We then combined the fits from different template libraries using the bayesian information criterion (BIC, Schwarz 1978). It is defined as:

$$\begin{array}{c}\hfill BIC=k·ln(n_{\mathrm{d}})-2ln(\widehat{L}),\end{array}$$(22)

where $\widehat{L}$ is the maximised value of the likelihood function of the model M, i.e. ${\displaystyle \widehat{L}=p(d|\widehat{{\xi }_{\max }})}$, where ${\displaystyle \widehat{{\xi }_{\max }}}$ are the parameter values that maximise the likelihood function, n_d is the number of data points (i.e. the number of unmasked spectral pixels), k is the number of parameters in the model.

The central value and statistical uncertainties are given by the BIC-weighted mean of the measurements obtained with the different libraries, while the systematic uncertainty is computed as the weighted-mean difference between the central value and each library measurement.

To test the robustness of the measurement against potential contamination of the CaK and G-band by sky emission lines, we repeated the measurement while masking the first ([3954;4000] Å) or the latter ([4275;4325] Å). We obtained measurements in excellent agreement with the original one and therefore use this one for the cosmographic analysis.

The final measurement, along with the corresponding covariance matrix, is shown in Table 1.

5. Line-of-sight analysis

5.1. Wide-field imaging

The HE 1104−1805 field was targeted by a multi-band imaging campaign to assess the contribution of line-of-sight structures to the lensing signal. The data acquired included u-band with a total exposure time of t_exp = 20 706 s with Megacam on the Canada France Hawaii Telescope (Program ID: 14AT01; PI Suyu); Subaru SuprimeCam imaging in the g (t_exp = 600 s), r (t_exp = 4800 s), and i (t_exp = 600 s) bands (Program ID: o14220; PI Fassnacht, the latter is shown in Fig. 6); and Subaru MOIRCS near-infrared imaging in the J (t_exp = 3600 s), H (t_exp = 5760 s), and Ks (t_exp = 2940 s) bands (Program ID: o15202; PI Fassnacht). For the current analysis of the line-of-sight, we only considered the SuprimeCam imaging.

Fig. 6. Wide-field imaging of HE 1104−1805 in the i-band. Red circles denote regions within 8″, 45″ and 120″ of the centre of the lens.

The SuprimeCam data were reduced in a two-phase process. The initial steps utilised the SDFRED reduction package to do the overscan and bias correction, flat-field correction, distortion and atmospheric distortion corrections, and masking of some bad regions. The second phase used custom Python-based scripts to mask bad pixels and to combine the individual exposures.

The SuprimeCam imager provides a field of view of approximately 30′ on a side, which is substantially larger than what is needed for the analysis presented in Sect. 5.4. We therefore used only a 4′ × 4′ section centred on the lens system to generate the galaxy catalogues. The catalogues were produced by running SExtractor (Bertin & Arnouts 1996) in dual-image mode, using the deep r-band data as the detection image and the i-band image as the measurement image. The photometric zeropoints were determined by comparing the instrumental magnitudes of stars in the field to the magnitudes measured by SDSS for those stars.

5.2. Perturbers and groups explicitly included in the lens model

As shown in Fig. 1, we identified nine perturbers within 5^″ of the lens in the HST imaging. Our initial step, hence, involved identifying the perturbers that we need to explicitly include to the lens model.

The ‘flexion shift’ was introduced by McCully et al. (2017) to determine which line-of-sight galaxies should be explicitly included in the lens modelling and which ones should be implicitly accounted for through the computation of κ_ext.

The specific influence of line-of-sight objects depends on whether the dominant-lens approximation holds, indicating that the critical density of these objects is much smaller than that of the primary deflector (Fleury et al. 2021).

The difference in lens image position induced by a given nearby galaxy is given by:

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{3}x=f\left(\beta \right)\frac{{({\theta }_{\mathrm{E},\mathrm{lens}}·{\theta }_{\mathrm{E},\mathrm{pert}})}^2}{{\theta }^3},\end{array}$$(23)

with θ_E, lens and θ_E, pert as the Einstein radii of the lens and perturber, and θ as the angular separation of the lens and perturber. The function f(β) equals 1 if the perturber is between the lens and the observer, while f(β) = (1−β)² if the perturber is behind the main lens. β is the pre-factor of the lens deflection in the multi-plane lens equation: β = (D_dpD_s)/(D_pD_ds). McCully et al. (2017) found that perturbers with Δ₃x ≤ 10⁻⁴ are accurately accounted for with a global external convergence term, whereas perturbers with a higher flexion shift should be individually modelled with a multi-plane treatment.

We used the M_* − σ^LOS relation obtained Zahid et al. (2016) on galaxies up to redshift 0.7 to estimate the mass-to-light ratio, M/L, of the P5 perturber using its velocity dispersion σ^LOS_P5 measured with the MUSE data (see Appendix B.2 for more details). The flux determined with a circular aperture around P5 in the F160W band converted into luminosity using a fiducial cosmology with H₀ = 70 km s⁻¹Mpc⁻¹ and Ω_m = 0.3. By applying the same M/L ratio to the other perturbers, we estimated their mass and σ^LOS. Given the uncertainty of the ${{\sigma }^{\mathrm{LOS}}}_{\mathrm{P}5}$ measurement, the uncertainty on the fiducial cosmological parameters and the fact that some perturbers are assumed to be further than the range studied by Zahid et al. (2016), can be neglected. By approximating the perturbers’ mass distributions with an SIS, we computed their Einstein radius, θ_E, SIS, with :

$$\begin{array}{c}\hfill {\theta }_{\mathrm{E},\mathrm{SIS}}=4\pi {\left(\frac{{\sigma }^{\mathrm{LOS}}}{c}\right)}^2\frac{D_{\mathrm{ps}}}{D_{\mathrm{s}}}\end{array}$$(24)

with D_ps and D_s the angular diameter between the source and perturber, and between the observer and the source, respectively.

We assumed that the perturbers without a redshift measurement lie between the closest perturber, P5 and the source. We therefore repeated the computation with redshifts uniformly drawn between z₅ and z_s, and took the mean, the 16th, and 84th percentile (1 − σ in a Gaussian distribution) of the resulting distribution as the final estimate.

We could then estimate the flexion shift caused by each perturber using Eq. (23). The results presented in Table 2 show that even when marginalising on their redshift, the flexion-shifts of these perturbers are all potentially larger than 10⁻⁴. They can therefore significantly influence the image’s position and should be included in the mass model.

In addition, the subtraction of the PSF, lens, and source light with a preliminary SIE model revealed the presence of a luminous component ∼0.6^″ from image B. Without a counter-image indicating that this feature belongs in the source, we treated this component as a perturber, hereafter P9. Because of its faintness, the spectra of P9 could not be extracted from the MUSE data cube, and its mass could not be estimated with the same technique employed for all other perturbers.

5.3. Identification of galaxy groups in the vicinity of the lens

Our search for galaxy groups along the line-of-sight relies on a catalogue of spectroscopic redshifts of field-of-view galaxies. To compile that catalogue, we have obtained multi-object spectroscopy of galaxies up to ∼2 arcmin using the FORS2 instrument mounted on the Antu VLT telescope (PID: 092.0515, PI: D. Sluse), and using GMOS at the Gemini-South Telescope (PID: GS-2013B-Q-28; PI: T. Treu). The observational setup, data reduction, and target selection strategy are similar to the one described in Sluse et al. (2017). Individual spectra typically cover the range 4500–9200 Å for FORS2 data and 4400–8200 Å for GMOS. Within a 30″ radius around the lens, we also extracted spectra of the sources present in the MUSE data used for kinematics (Sect. 4.2). We used 1D cross-correlation methods to measure redshifts: the xcsao task implemented in IRAF was used for the FORS data, and MARZ (Hinton et al. 2016) was used for the other data sets. This yields an ensemble of 131 reliable unique redshifts (96 with FORS2, 18 with GMOS, 17 with MUSE) distributed typically in the range z ∈ [0.04, 0.8]. We complemented the catalogue with redshift measurements from Momcheva et al. (2015) who targeted galaxies as faint as I ∼ 21 mag within a radius of 15 arcmin of HE 1104−1805. This sample adds 387 unique redshifts to the catalogue, out of which 47 are within a radius of 6 arcmin of the lens, i.e. the region where we expect groups to have first-order contribution to the main lens gravitational potential.

The identification of galaxy groups is done using the algorithm described in Sluse et al. (2017), Sluse et al. (2019), and Buckley-Geer et al. (2020), and we refer the reader to these publications for details. This algorithm is based on the group-finding algorithms of Wilman et al. (2005) and Ammons et al. (2014). For this measurement, we have computed the group centroid from the luminosity-weighted positions of the group members. The groups likely to perturb the modelling of the main lens are summarised in Table 2. Groups 4, 5, and 6 highlight the group finder’s sensitivity to initial conditions and interlopers, with some galaxies appearing in multiple groups. We therefore included them separately in the lens models (see Sect. 6). In Fig. D.1, we show, for each identified group, the right ascension and declination of the accepted and rejected member galaxies as well as the distances and velocities relative to the final group centroid.

5.4. External convergence (κ_ext) measurement

Similar to the previous time-delay cosmography studies (e.g. Rusu et al. 2017; Birrer & Amara 2018; Rusu et al. 2020), we used the number count technique to determine κ_ext. This quantity determines the contribution of objects not included in the lens modelling. Number counts involve measuring the galaxy number density near the lens as a summary statistic and comparing it to reference fields. This comparison helps determine whether the LOS is over- or under-dense compared to the average background (e.g. Fassnacht et al. 2011; Greene et al. 2013; Wells et al. 2023).

It can be summarised in four main steps:

• Lines of sight (typically several arcminutes wide squares) are drawn from a large comparison field and compared to the LOS of the lens. Modern survey datasets covering hundreds to thousands of square degrees ensure a sufficiently large comparison field to avoid sampling bias. The photometric catalogues of the fields are cleared from objects more distant than the lensed source and fainter than the used magnitude cut on the i < 24 (e.g. Fassnacht et al. 2011). This cut is designed to consider only objects with reliable photometry.

• The weight value of a given LOS W_i is then computed as the ratio of the weighted number counts for galaxies in the lens field to the same statistic in the reference field.

  $$\begin{array}{c}\hfill W_i=\frac{{\mathrm{\Sigma }}_j{w}_{j,\mathrm{lens}}}{{\mathrm{\Sigma }}_j{w}_{j,i}}\end{array}$$(25)

  where j indexes the galaxies in a given LOS and w is the weight of a galaxy. Different types of weights can be considered (e.g. Greene et al. 2013; Rusu et al. 2017; Wells et al. 2023) such as the galaxy’s distance to the center (1/r_j), its potential (m_j/r_j), its redshift (z_j(z_s − z_j)) or simply be the same for each galaxy (w_j = 1). The weight value obtained, therefore, acts as a tracer of the external convergence of the lens LOS.

• Similar weight values are computed in simulated fields with known κ_ext (Hilbert et al. 2009). Previous works (e.g. Rusu et al. 2017; Wells et al. 2023), have shown that using the Millennium Simulation yields a reliable estimate of κ_ext because it provides catalogues of galaxies for several values of the external convergence at sufficient points to represent the Universe accurately. Since the weight value relies on ratios, many dependencies on the simulation’s underlying cosmological parameters are expected to cancel out.

• The posterior probability of κ_ext can hence be computed as follows

  $$\begin{array}{c}\hfill p({\kappa }_{\mathrm{ext}}|d)={\displaystyle \int p_{\mathrm{sim}}}({\kappa }_{\mathrm{ext}}|W_i)p(W_i|d){\mathrm{\Pi }}_{i}d{W}_i.\end{array}$$(26)

  Here, p(W|d) represents the probability distribution of the weighting scheme W given the data, and p_sim(κ_ext|W) denotes the probability distribution of κ_ext in the simulated dataset, conditioned on a specific value of the weight W. Greene et al. (2013) showed that the combination of multiple evaluations relying on different weights (while adequately taking into account covariances between estimates) significantly improves the precision and robustness of the measurement.

Fig. 8. κext measurement showing that the line-of-sight of HE 1104−1805 is along a slightly underdense region of the Universe. Therefore, the fact that multiple perturbers are near the main lens is coincidental and the approximation κext ≈ κs holds.

As detailed previously, the treatment of perturbers within 5^″ from the lens is part of the explicit lens model. We, therefore, masked this region from the number counting and removed the galaxies forming the perturber groups presented in Table 2, as in Wells et al. (2023), we limited the field to 120^″ around the lens. The weights used are the number counts w₁ and the inverse distance w_1/r of objects. Their computation is displayed in Fig. 7 and Fig. 8 displays the final measurement ${\kappa }_{\mathrm{ext}}=-0.{019}_{-0.020}^{+0.058}$. Our procedure, in fact, determines κ_s, the convergence up to the source which we used as a proxy for κ_ext. The external convergence is technically computed using the observer-deflector and deflector-source convergences (κ_d and κ_ds) through κ_ext = (1−κ_d)/(1−κ_s) ⋅ (1−κ_ds). However, Tang et al. (2025) showed that ignoring this correction affects the H₀ measurement only when the LOS up to the deflector is significantly over- or under-dense, which is not the case here. At the population level, the expected bias on H₀ due to this approximation is about 0.1%, which is negligible compared to the present 4.6% precision of TDCOSMO25.

Fig. 7. Position and weight (size of the dot) of galaxies between us and the source. Black circles indicate the 45″ and 120″ apertures. Blue dots correspond to objects with i-magnitudes less than 24.

Therefore, this measurement is used for the rest of the cosmographic analysis, without further correction, in compliance with the perturber inclusion strategy.

6. Lens modelling

6.1. HST imaging

In the context of H0LiCOW (Suyu et al. 2017), Program 12889 (P.I: S.H. Suyu) took 12 exposures of 36 seconds along with 24 exposures of 599 seconds of HE 1104−1805 with wide-field camera 3 (WFC3) in the F160W⁷ band. A 4-point dither pattern with parallelogram shapes having half-pixel shifts was applied in order to reach a resolution of 0.08^″ per pixel. The imaging displayed in the left panel of Fig. 1 showcases the main challenges for modelling this system. First, the point source images are particularly bright with the PSF wings of both images aligned with the centre of the galaxy. Accurate modelling of the PSF’s wings is, therefore, crucial to model the light of the lens. Secondly, the quasar host’s light is very dim and distinguishable only after proper subtraction of the lens light (see Fig. C.2). This makes the lens light fitting a primordial and sensitive step to properly reconstruct the shape of the arc, a key constraint to mass modelling. Thirdly, the system’s environment is particularly crowded with eight luminous galaxies within 15″ of the main lens. The system is also characterised by a large and asymmetric separation of images A and B positions from the lens centre of 1.17^″ and 2.05^″, respectively.

6.2. Setup and workflow

6.2.1. Main deflector mass and light model

Following the methodology of the most recent TDCOSMO lens models (e.g. Shajib et al. 2020, 2022), we modelled the lens galaxy with two families of models: power law and a composite alternatively.

In the first case, the mass is parametrised as a power-law elliptical mass distribution (PEMD) for which the convergence at the position (θ₁,θ₂) in the frame aligned with the major and minor axis of the deflector is given by

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{PL}}({\theta }_1,{\theta }_2)=\frac{3-{\gamma }_{\mathrm{pl}}}{2}{\left(\frac{{\theta }_{\mathrm{E}}}{\sqrt{q_m{\theta }_1^2+{\theta }_2^2/q_m}}\right)}^{{\gamma }_{\mathrm{pl}}-1},\end{array}$$(27)

where θ_E is the Einstein radius, γ_PL the slope of the profile, and q_m its axis ratio. This convergence map is then rotated by the position angle ϕ_m to belong to the on-sky coordinate frame.

The lens light is modelled by a superposition of two Sérsic profile which gives the light intensity at the position (θ₁,θ₂) following

$$\begin{array}{c}\hfill I_{\mathrm{Sersic}}({\theta }_1,{\theta }_2)=I_{\mathrm{eff}}exp\{-b_n[{\left(\frac{\sqrt{{\theta }_1^2+{\theta }_2^2/q^2}}{{\theta }_{\mathrm{eff}}}\right)}^{1/n}-1]\}\end{array}$$(28)

with θ_eff the effective radius (i.e. the product average of semi-major and semi-minor axis), also defined as the half-light radius, thanks to the normalising factor b_n, I_eff the intensity at θ_eff. q is its axis ratio, and n the Sérsic index.

The positions of the two Sersic profiles are tied together throughout the fit, and the position of the mass component’s centre is fixed on the light component.

The chameleon profile can be parametrised as the difference between two non-singular isothermal ellipsoid (NIE) profiles (e.g. Suyu et al. 2014)

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{Chm}}({\theta }_1,{\theta }_2)=\frac{{\kappa }_{0,\mathrm{Chm}}}{(1+q)}[\frac{1}{\sqrt{{\theta }_1^2+{\theta }_2^2/q^2+4w_{\mathrm{c}}^2/{(1+q)}^2}}\\ \hfill -\frac{1}{\sqrt{{\theta }_1^2+{\theta }_2^2/q^2+4w_t^2/{(1+q)}^2}}]\end{array}$$(29)

where q is the axis ratio and w_c and w_t are sizes of the cores of the two NIE with w_t >w_c.

For the composite model, the baryonic mass is described by a double chameleon profile (i.e. two superposed chameleon profiles presented in Eq. 29).

The dark matter halo is modelled with an elliptical Navarro-Frenk-White (NFW) profile, whose volumic density is given by

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_0}{\left(\frac{r}{R_{\mathrm{s}}}\right){(1+\frac{r}{R_{\mathrm{s}}})}^2}\end{array}$$(30)

where ρ₀ is a normalisation and R_s is the scale radius.

Following the measurement of Gavazzi et al. (2007) on galaxies that have a similar σ^LOS as the one measured in this system, we expect the NFW’s scale radius to be within R_s = 58 ± 8 kpc. To prevent any prior bias, we took this value as a starting point within a wide uniform prior R_s ∼ 𝒰(25 100) kpc.

We used a double chameleon profile to model the light, and the fitted parameters are fixed to the parameters of the baryonic mass model.

To avoid the lens light modelling from being biased by background light while still predicting accurate flux at the position of the arc, we masked pixels further than 4.3^″ from the lens. We adjusted the mask locally to mask the light of all the perturbers and masked the lens’s central pixels because the light is not fitted perfectly in this region. Masking the few pixels close to the centre of the lens prevents the model from predicting source flux at this position, which can possibly bias the inferred mass profile.

6.2.2. Source light model

The visible source light is smooth, and we modelled it with a single elliptical Sérsic profile in both families of models. Contrary to previous TDCOSMO studies, no further structure is detected in the lensed arc. Thus, it is not required to add shapelets to describe extra complexity.

6.2.3. Perturbers’ model

To include the perturbers in the mass model, we added an SIS component to the mass model for each perturber. Even though we did not fit their light in the final model, a preliminary Sérsic fit gives us their position. Since P5 and P6 have different redshifts than the lens, we used the multi-lens-plane formalism introduced in Eqs. (4) and (8). For the other perturbers considered, we assumed that they lie in the same plane as the lens. We considered the perturbers P2, P5, and P6 as they have the highest flexion shifts, and P9 because of its proximity to image B. Furthermore, the perturber groups listed in Table 2 are parametrised with NFW profiles. We estimated their scale radii using their virial mass (based on measured R₂₀₀ and v_disp) and estimation of the NFW concentration parameter, R₂₀₀/R_S from Duffy et al. (2008). The groups 4, 5, and 6 are alternatively included in the model.

6.2.4. PSF modelling

To complete the light model, each component is convolved by the PSF. As a first guess of the PSF, five stars of the field were stacked together using the PSFr software (Birrer et al. 2022).

Since the spectrum of the stars in the field used for the PSF initial guess does not match the quasar one, we finalised the PSF model with four successive repetitions of the following sequence with the Lenstronomy software⁸ (Birrer & Amara 2018):

1. mass and light model parameter optimisation using Lenstronomy’s built-in particle swarm optimiser routine (e.g. Eberhart & Kennedy 1995; Birrer & Amara 2018)

2. iterative reconstruction of the PSF on lens and source subtracted residuals. We also adjusted the PSF noise map in a radius of 0.5^″ around its centre to adapt the constraining power of these pixels.

6.2.5. Posterior sampling and combination of different models

After obtaining a satisfactory PSF, shown in Fig. C.1, we used Lenstronomy’s built-in MCMC routine relying on EMCEE (Foreman-Mackey et al. 2013) to sample the light and mass parameters (ξ_light and ξ_mass) posterior by maximizing the likelihood expressed as

$$\begin{array}{c}\hfill \mathcal{P}({\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}}|\mathrm{d})\propto \mathcal{L}(\mathrm{d}|{\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}})p({\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}}),\end{array}$$(31)

where ℒ(d|ξ_light, ξ_mass) is the likelihood of the parameters given the data and p(ξ_light, ξ_mass) is the prior on the parameters. During this step, the data used to constrain the posterior are the unmasked pixels of the HST observation 𝒟_Img and the time delay, Δt_AB. The likelihood maximised by the sampling is therefore computed with

$$\begin{array}{cc}\hfill \mathcal{L}(\mathrm{d}|{\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}})& =\mathcal{L}(\mathrm{Img}|{\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}})\hfill \\ \hfill & +\mathcal{L}(\mathrm{\Delta }{t}_{\mathrm{AB}}|{\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}})\hfill \end{array}$$(32)

which is, up to an additive constant,

$$\begin{array}{c}\hfill =log[-\frac{1}{2}{\displaystyle \sum _i^{N_{\mathrm{pix}}}}\frac{{({\mathcal{D}}_{\mathrm{Img},\mathrm{i}}-M_{\mathrm{Img},\mathit{i}})}^2}{{\sigma }_i^2}]+log[-\frac{1}{2}\frac{{(\mathrm{\Delta }{t}_{\mathrm{AB}}-M_{\mathrm{\Delta }{t}_{\mathrm{AB}}})}^2}{{\sigma }_{\mathrm{\Delta }{t}_{\mathrm{AB}}}}],\end{array}$$(33)

where M_Img and M_ΔtAB are the model-predicted image and time delay of the system; σ_i is the pixel’s noise, and σ_ΔtAB the uncertainty on the time delay. We used ten times as many walkers as there are parameters in the model, and we used the last 1000 out of 20 000 iterations to construct the posterior of a given modelling setup. The best reconstructions of the HST imaging with power-law and composite models shown in Fig. 9 demonstrate that all relevant features in the light are correctly predicted. The convergence profile obtained by each model family in every perturber-inclusion scenario is displayed in Fig. 10. We notice that, with both model families, the convergence of the perturbers P5, P6, and G1 varies substantially between the different scenarios. This can be attributed to the fact that these components are aligned in the same direction relative to the main lens, resulting in a degenerate effect on the system. Conversely, the convergence of the perturbers P2, P9, and groups G2, G3 are more stable as they are more azimuthally spread. We also note that in every scenario, the composite model is cuspier than the power-law one.

Fig. 9. Reconstructed image of HE 1104−1805 using the best power-law and composite models. From left to right: The imaging data, reconstructed image, and normalised residuals are displayed for the power law model (PEMD + Shear + perturbers 2, 5, 6, 9 + groups 1, 2, 3, and 4) and composite models (NFW + Double chameleon + Shear + perturbers 2, 5, 6, 9 + groups 1, 2, 3, and 4).

Fig. 10. HE 1104−1805 modelled convergence κ with the two families of models and different perturbers and groups inclusion strategies. The top row shows a wide-field including all groups and the bottom row zooms in on the region closer to the main lens; white lines indicate the area enclosing 25, 50, and 75% of the total mass.

The BIC values computed with the maximum likelihood of each model are presented in Table 3. The variance of the BIC is σ_ΔBIC = 73 within the power-law models and σ_ΔBIC = 48 within the composite models.

Similarly to Birrer et al. (2019) and Shajib et al. (2022), we determined the weight of each posterior distribution following

$$\begin{array}{c}\hfill \mathcal{W}=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathrm{\Delta }BIC}}{\displaystyle {\int }_{-\infty }^{\infty }}f(x)\exp (-\frac{{(BIC-x)}^2}{2{\sigma }_{\mathrm{\Delta }BIC}})\mathrm{d}x,\end{array}$$(34)

where f(x) is the evidence ratio function given by

$$\begin{array}{c}\hfill f\left(x\right)\equiv \{\begin{array}{cc}1\hfill & \mathrm{if}x\le {\mathrm{BIC}}_{\min },\hfill \\ \exp ({\mathrm{BIC}}_{\min }-x)\hfill & \mathrm{if}x>{\mathrm{BIC}}_{\min },\hfill \end{array}\end{array}$$(35)

with BIC_min the BIC of the reference model with the minimum BIC value.

To create a single posterior out of N multiple models, we summed all the posteriors using their normalised weights

$$\begin{array}{c}\hfill \mathcal{W}=\frac{{\sum }_n{\mathcal{W}}_{n}p({\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}}|d)}{{\sum }_{n}max({\mathcal{W}}_n)}.\end{array}$$(36)

This procedure is done separately for the composite and power law models, and we compare the results in Fig. 11. For both models, the Einstein radius, the slope, and the time-delay distance are robust to the different perturber configurations, whereas the external shear shows multi-modality because the G6 group centre is azimuthally significantly different from groups G4 and G5 (see Fig. 10). A blind comparison of the two modelling families allowed us to see that a 2.1-σ disparity between the two families on γ induces a 1.3-σ difference in the Fermat potential and 1-σ for the time-delay distance between the two families. The external shear we inferred for HE 1104−1805 is higher than in other TDCOSMO lenses. As shown by Etherington et al. (2023), model-predicted γ_ext in power-law + shear configurations is often overestimated even in perturber-free systems because it is degenerate with potential azimuthal structures, which have a negligible impact on H₀ at the population level (Van de Vyvere et al. 2022). In Section 5.2, we considered only nine perturbers, out of which five are explicitly modelled, but many more are visible in the right panel of Fig.1. These perturbers are expected to induce external shear on the main lens and induce a large γ_ext. More marginally, the positions of the centre of mass of the groups of pertubers, G1-G6, are poorly constrained and might also inflate the value of γ_ext.

Fig. 11. BIC weighted comparison of the power-law and composite models. For each family, we marginalised over the three posteriors given by the three perturber inclusion scenarios, which results in the multi-modality of some of the parameters such as γ. We obtained ${\theta }_{\mathrm{E}}^{\mathrm{PL}}=1.356\pm 0.115$, γpl = 1.88 ± 0.04, ${\gamma }_{\mathrm{ext}}^{\mathrm{PL}}=0.14\pm 0.01$, ${\theta }_{\mathrm{E}}^{\mathrm{Comp}}=1.37\pm 0.057$, ${\gamma }_{\mathrm{ext}}^{\mathrm{Comp}}=0.12\pm 0.02$, and γComp = 2.03 ± 0.08.

7. Cosmographic inference

We combined the mass and light posteriors obtained with the lensing and time delay constraints with the resolved kinematics measured in Section 4.2 to measure D_d, D_Δt, and H₀.

Following the methodology of Shajib et al. (2023), we importance sampled the posterior distribution of the lens model obtained with Eq. (36), ξ_light, ξ_mass together with the cosmological distances D_d and D_Δt given the three radial bins of the velocity dispersion σ^LOS and the external convergence κ_ext measured

$$\begin{array}{c}\hfill \mathcal{P}({\xi }_{\mathrm{light}},{\xi }_{\mathrm{mass}},D_{\mathrm{d}},D_{\mathrm{\Delta }t}|{\sigma }^{\mathrm{LOS}},{\kappa }_{\mathrm{ext}})\\ \hfill \propto \mathcal{L}({\sigma }^{\mathrm{LOS}}|{{\sigma }^{\mathrm{LOS}}}_{\mathrm{modelled}})p({\kappa }_{\mathrm{ext}})\end{array}$$(37)

with p(κ_ext) the prior on κ_ext obtained in Section 5.4.

The likelihood of the modelled 3-bin velocity dispersion is computed with the following

$$\begin{array}{cc}& \mathcal{L}({\sigma }^{\mathrm{LOS}}|{{\sigma }^{\mathrm{LOS}}}_{\mathrm{modelled}})\hfill \\ \hfill & \propto exp[-\frac{1}{2}{({\sigma }^{\mathrm{LOS}}-{{\sigma }^{\mathrm{LOS}}}_{\mathrm{modelled}})}^{\mathrm{T}}{\mathcal{C}}^{-1}({\sigma }^{\mathrm{LOS}}-{{\sigma }^{\mathrm{LOS}}}_{\mathrm{modelled}})],\hfill \end{array}$$(38)

where ${{\sigma }^{\mathrm{LOS}}}_{\mathit{modelled}}$ is the 3-bins velocity-dispersion computed with the JamPy package⁹ (Cappellari 2008, 2020) through Eqs. (17) and (19). 𝒞 is the covariance matrix between σ^LOS measured in each bin determined in Section 4.2. In practice, the software uses multi-gaussian-expansions (Emsellem et al. 1994; Cappellari 2002, MGE) fits of the mass and light profiles to be able to deproject these surface densities into 3D ones in a straightforward way (Monnet et al. 1992). The light MGE fit is based on the lens light profile corresponding to the ξ_light with maximum likelihood. The PSF is modelled as a Moffat with β = 1.96 and FWHM = 0.54^″ based on the mean value of these parameters along the wavelength range used for the σ^LOS measurement (i.e. [6000-9000] Å) that is shown in the bottom panel of Fig. 4.

As explained in Section 2.2, we adopted a uniform prior 𝒰(−0.12, 0.13) on β_ani. The resulting kinematic fit for each modelling family is shown on the left panel of Fig. 12 while the associated posterior distribution of D_d and D_Δt is shown in the right panel. The kinematic modelling performed without allowing for the internal mass-sheet shows that the steepness of the composite model fails to correctly fit the observed velocity dispersion profile, especially in the innermost bin, where the best fit predicts 301 km s⁻¹, i.e 1.6 σ away from the observed value. This results in a 2.1 σ tension between the inferred D_d. The likelihood-weighted combination assigns 95 % of the final posterior to the power-law model, effectively discarding the composite model measurement.

Fig. 12. Left panel: Kinematic fit using power-law and composite models. Red and blue central dots show the best-fit values, while error bars are the 16th and 84th percentiles of 100 predicted velocity dispersions. Right panel: Posterior of Dd and DΔt using kinematics based on the power-law and composite mass models. Using likelihood weighting we measured $D_{\mathrm{\Delta }t}={5083}_{487}^{+504}$ Mpc and $D_{\mathrm{d}}={1921}_{-236}^{+254}$ Mpc.

Using the posterior distribution of D_d and D_Δt we then inferred H₀ in a flat-ΛCDM cosmology with the purposedly conservative uniform priors 𝒰(0, 150) km s⁻¹Mpc⁻¹ on H₀ and 𝒰(0.05, 0.5) on Ω_m. The resulting posterior distributions obtained for both models and the combination are displayed in Fig. 13. As with every TDCOSMO H₀ measurement, the entirety of the analysis was performed blindly to avoid confirmation bias. The final value was revealed on July 8th of 2025, $H_0=64.2_{-5.0}^{+5.8}$ km s⁻¹Mpc⁻¹, i.e. an 8.5% precision from a single doubly imaged lens.

Fig. 13. H0 measurement with HE 1104−1805 using the joint constraint of DΔt and Dd measured with the power-law, composite models, and both models combined. Our final measurement using the likelihood weighted combination of both models is $H_0=64.2_{-5.0}^{+5.8}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ for λint = 1.

8. Discussion and conclusion

Table 4 summarises all the cosmography-related measurements achieved for HE 1104−1805.

8.1. Time-delay

The measurement of the time delay is marginally higher than previous estimates (e.g. Poindexter et al. 2007; Morgan et al. 2008, measured Δt_AB = 162 ± 6.1 days). It relies on the longest and best-sampled light curve of that system, with the two most recent (ECAM and WFI) datasets having at least twice the cadence as the previous SMARTS dataset. It was estimated with a data-driven algorithm whose accuracy was assessed on multiple simulated and observed light curves, which is crucial, as we saw that the extrinsic variability is particularly challenging to model in this system. We therefore have the most robust estimate for this object.

8.2. Resolved kinematics

The main limitation of the resolved velocity dispersion measurement comes from the proximity and brightness of image A to the lens galaxy. Nevertheless, by subtracting the quasar images, we were able to extract spectra from three concentric bins with S/N > 14 Å⁻¹. Using state-of-the-art methodology to constrain the effect of stellar template fitting hyperparameters (polynomial degrees and stellar template libraries), we provided the velocity dispersion measurement in each bin. We ensured the robustness of these measurements against different spectral and spatial masks and proved that the lens is a slow rotator (see Appendix A.2). The latter justifies the use of spherically aligned JAM for the kinematic modelling.

8.3. Image and kinematic modelling

The image reconstruction was performed with the two main parameterisation families: power-law and composite. Three configurations of nearby perturbers were accounted for individually, with the addition of SIS and NFW components, and the resulting measurements were marginalised upon these configurations.

The main difference between the two models resides in the slope of the radial mass profile γ. Such a difference was already observed in modelling the quadruply lensed quasar WGD 2038−4008 (Shajib et al. 2022), where the system’s compactness only allowed to probe the inner regions of the mass profile. In the present case, the resolved velocity dispersion heavily favours the power-law model over the composite for this system. This may be because HE 1104−1805’s configuration offers too few constraints for a composite model, thus leading to potential biases or simply that the power law model is a better description over the range probed. The discrepancy between the external convergence, κ_ext, and the external shear, γ_ext, inferred independently might be counter-intuitive in light of simulations’ properties (e.g. Tang et al. 2025). However, in this work γ_ext is primarily encompassing the effect of perturbers within 5^″ of the main deflector, whereas κ_ext is inferred from galaxies up to 120^″ away, excluding the ones within 5^″. As discussed in Section 6.2.5, the large value of γ_ext may arise from genuine external shear or simply be inflated by inherent lens-modelling degeneracies. We could not disentangle the contribution of each effect with the available data, nor rule out unmodelled convergence. However, since the brightest perturbers we considered are already low-mass (see upper part of Table 2), any additional convergence from fainter perturbers should be small.

8.4. Implications for time delay cosmography

This work shows that the high precision of the time-delay measurement and the use of resolved kinematics allowed us to reach a precision of 8.5% on H₀ with a double lens quasar, which is comparable to previous TDCOSMO single lens measurements ranging from 4% to 28% (Millon et al. 2020c; Wong et al. 2024). The main limitation of this work is the fact that we fixed the internal mass-sheet degeneracy parameter, λ_int = 1. This parameter will be constrained during the next joint hierarchical modelling of the TDCOSMO lens sample. The H₀ value and precision of this system are therefore bound to change as H₀ → H₀λ_int, (δH₀ / H₀)² → (δH₀ / H₀)² + δλ_int².

Overall, this work is a stepping stone towards building a sample of double-lens quasar-based H₀ measurement, which is vital to investigate potential biases of the TDCOSMO method and to accelerate progress towards a 1% measurement of H₀.

Acknowledgments

E.P. is supported by JSPS KAKENHI Grant Number JP24H00221. F.C. is 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 (COSMICLENS: grant agreement No 787886). C.D.F. acknowledges support for this work from the National Science Foundation under Grant No. AST-2407278. A.G. acknowledges support by the SNSF (Swiss National Science Foundation) through a mobility grant P500PT_211034. M.M. acknowledges support by the SNSF (Swiss National Science Foundation) through mobility grant P500PT_203114 and return CH grant P5R5PT_225598. D.S. acknowledges the support of the Fonds de la Recherche Scientifique-FNRS, Belgium, under grant No. 4.4503.1. S.B. acknowledges support by the Department of Physics and Astronomy, Stony Brook University. C.L. acknowledges funding from the European Union’s Horizon Europe research and innovation programme under the Marie Sklodovska-Curie grant agreement No. 101105725. A.J.S. received support from NASA through STScI grants HST-GO-16773 and JWST-GO-2974. S.H.S. thanks the Max Planck Society for support through the Max Planck Fellowship. This work is supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. T.T. acknowledges support by NSF through grants NSF-AST-1906976, NSF-AST-1836016 and NSF-AST-2407277, and from the Moore Foundation through grant 8548. K.C.W. acknowledges support by JSPS KAKENHI Grant Numbers JP24K07089, JP24H00221. V.M. acknowledges support from ANID FONDECYT Regular grant number 1231418, Centro de Astrofísica de Valparaíso CIDI 21, and ANID, Millennium Science Initiative, AIM23-0001. This work used observations collected at the European Southern Observatory under ESO programmes 092.0515(B) and 0102.A-0600(A). This work used Python packages Numpy (Harris et al. 2020), Astropy (Astropy Collaboration 2022), and Scipy (Virtanen et al. 2020).

References

Appendix A: Velocity dispersion measurement additional material

A.1. PSF subtraction in the MUSE cube

As shown in Fig. A.1, we first performed the fit on ten stacked wavelengths to constrain the parameters as a function of wavelength when performing the frame-by-frame fit. The wavelength change of the Moffat index β and the Full-Width-Half-Maximum (FWHM) are consistent with expected MUSE PSF modelling (e.g. Fusco et al. 2020; Weilbacher et al. 2020). The relative instability of the ellipticity of the PSF can be attributed to the abrupt changes in the quasar and lens light shape due to emission and absorption lines.

Fig. A.1. PSF parameters fitted across wavelength. The values obtained when using stacks of 12.5 Åare used as constraining prior for the wavelength-by-wavelength fit. The shaded area indicates the range of wavelength used for the σLOS measurement, and the dotted lines highlight the mean value of the parameters in this range.

A.2. Is the lens rotating?

To investigate bias in the σ^LOS measurement due to possible rotation of the lens galaxy, we designed three spatial masks cutting half of the lens light area as shown in the left panel of A.2. The three masks are chosen to ensure an equivalent S/N for the three configurations while minimising the spatial overlap. With the same methodology as described in Section 4.2, we then measured the line-of-sight velocity of each bin $v_{1,2,3}^{\mathrm{LOS}}$. The bottom panel of Fig.A.2 shows that there is no sign of ordered rotation. These observations are, hence, consistent with the assumption that the lens galaxy is a slow rotator.

Fig. A.2. Top panel: three different masked MUSE cube cutouts. Bottom panel: measured velocity of the three bins in each case.

Appendix B: MUSE by-products

B.1. Spectra extracted from the MUSE data

We present here the spectra of the center of the lens galaxy (S/N ∼ 30Å⁻¹), perturbers P5 (∼ 11Å⁻¹), P6 (∼ 7Å⁻¹), and images A (∼ 529Å⁻¹) and B (∼ 266Å⁻¹) of the quasar extracted at the positions displayed in Fig.4.

Fig. B.1. Integrated spectra from the MUSE data cube with apertures shown in Fig. 4. The binned flux is computed by imposing a minimal S/N of 5 per wavelength bin. From top to bottom, Lens centre: from left to right, the stellar absorption lines CaH, CaK, H[[INLINE557]], and G-band (orange dashed lines) are consistent with the previously determined zd = 0.729. P5: from left to right, the galactic emission lines Ne III, He I, H[[INLINE559]], and H[[INLINE560]] are highlighted with green dashed lines and allowed us to measure z5 = 0.3575. We also identify the stellar absorption lines H[[INLINE562]], HEI, CaH, and CaK allowing a velocity dispersion measurement. P6: from left to right the lines Mg I, Ca II O II, H[[INLINE563]], and O III doublet allowed us to measure z6 = 0.505. Quasar images: from left to right, we recognise the C IV, C III, C II, and Mg II emission lines.

Fig. B.2. Measurement of the P5 stellar kinematics. Point-estimate of σLOSP5.

B.2. Velocity dispersion of the perturber P5

We applied the same stellar template fitting process as performed with the lens galaxy spectra to the spectrum of P5 shown in Fig. B.1. In this case, we used the observed wavelength range [(4600,4700):(5600,5700)] Å. As shown by Fig. B.2, the absorption lines are not altered by the sky but by galactic gas emission lines. We thus incorporated those in the model without testing different masks. We obtained the ${{\sigma }^{\mathrm{LOS}}}_{P5}=123\pm 71\mathrm{km}{\mathrm{s}}^{-1}$. According to Knabel et al. (2025), spectra with relatively low S/N, such as this one, may be affected by systematic bias up to ∼4%, which is negligible compared to the uncertainty of the measurement.

Appendix C: Lens modelling, additional material

C.1. PSF reconstruction for HST imaging

The Fig.C.1 shows the PSF obtained after the iterative reconstruction procedure described in Section 6.2.4.

C.2. Image reconstruction

We now show the decomposition of the imaging reconstruction, i.e. the lens light and source light obtained with a power-law in Fig.C.2 and composite in Fig. C.3 models.

Fig. C.1. From left to right: Initial guess of the PSF from stacked stars in the field, final PSF obtained, and the difference between both.

Fig. C.2. HE 1104−1805 reconstructed light component with power-law. From left to right: The light component unconvolved, convolved with the PSF, and its corresponding feature in the imaging data. The top row shows the lens light profile, and the bottom row shows the lensed source light profile. This figure shows that the light-components model corresponds to expectations based on the observation. The centre of the PSF and of the lens are not ideally modelled, which justifies the noise map increase in the first region and the masking of the second region.

Fig. C.3. HE 1104−1805 reconstructed light component with power-law. From left to right: The light component unconvolved and then convolved with the PSF and its corresponding feature in the imaging data. The top row shows the lens light profile, and the bottom row shows the lensed source light profile. This figure shows that the light-components model corresponds to expectations based on the observation. The centre of the PSF and of the lens are not ideally modelled, which justifies the noise map increase in the first region and the masking of the second region.

Appendix D: Galaxy perturber groups members

Fig. D.1 displays the position and velocities of the galaxy groups included as perturbers in the lens model.

Fig. D.1. Galaxy groups identified in the spectroscopic sample of galaxies in the field of HE 1104−1805. For each group, the upper plot shows the positions of the candidate member galaxies associated with that group relative to the lens galaxy, with rejected group members represented as red squares and accepted group members represented as green circles. The lens galaxy (star) and group centroid (cross) are also displayed. The R200 radius of the group is represented by a solid line, while the dashed circle represents the angular separation cut of the group-finding algorithm in its final iteration. The lower plot shows the observer-frame velocity of individual member galaxies relative to the group centroid as a function of that galaxy’s angular distance from the centroid. Galaxies that passed the iterative algorithm described in Buckley-Geer et al. (2020), Sluse et al. (2017), and Sluse et al. (2019) are shown in green, while trial galaxy members that were cut through the algorithm are shown in red. Horizontal error bars represent the measurement error for each galaxy. The final observer-frame velocity dispersion and angular separation cuts from the group-finding algorithm are presented as dashed and dotted lines, respectively. We also show 1-D histograms and rug plots of the velocity and distance distributions of the member galaxies. The 1-D histograms are produced using a kernel density estimate (KDE) with a bandwidth determined using Scott’s Rule. In the 1D velocity histogram, the dashed blue line shows a Gaussian with a width equal to the observer-frame velocity dispersion of the group.

All Tables

All Figures

Fig. 1. Left panel: Stack of 280 WFI exposures of the field of view used for light curve extraction. Stars nominated PSF1 to PSF5 are used to model the PSF of each exposure, and stars designated as N1 to N5 are used to normalise the flux between each exposure. Right panel: HE 1104−1805 imaging in the filter F160W band using HST WFC3. The main lensing galaxy is denoted as G, whereas the main perturbers within 5″ of the lens are numbered from P1 to P9.

In the text

Fig. 2. Top panel: HE 1104−1805 R-band light curve obtained by joining three datasets: SMARTS (yellow and blue squares), ECAM (purple and green circles), and WFI (brown and pink triangles). For clarity, the light curve of B was shifted by –1.05 mag. The inset zooms on the WFI dataset span, showcasing its superior sampling to the ECAM dataset. Bottom panel: Example of a simultaneous fit of an intrinsic light curve with η = 45 days and nml = 15. The time shift obtained by this fit gives a point estimate of the system’s time delay. By repeating such measurements with randomised starting points 800 times, we obtained the time-delay measurement of this (η, nml) configuration.

In the text

Fig. 3. Left panel: Measurements of the time delay for different configurations of (η, nml) with the merged ECAM+SMART+WFI dataset. Right panel: Measurement of ΔtAB for the different datasets. Since the datasets overlap, they cannot be considered fully independent. We therefore use our estimate on the joint dataset ECAM+SMARTS+WFI (green data point) as our final estimate. In both panels, the ‘combined τ = 0.5’ value was obtained by marginalising over the measurements listed above it.

In the text

	Fig. 4. MUSE data cube summed over all wavelengths. The inset displays the point-source subtracted cube, revealing the lens galaxy light. The apertures used to extract individual spectra of the quasar images, lens centre, and perturbers shown in Fig. B.1 are represented by orange circles.
In the text

	Fig. 5. Velocity dispersion point-estimate in the three radial bins of HE 1104−1805’s lens galaxy. Left panel: Mean of the MUSE data cube after PSF subtraction and masking with overlaid bin numbers and contours. Right panel: Example of a pPXF fit of the integrated spectra of each bin, where data masked for the fit are marked in grey.
In the text

	Fig. 6. Wide-field imaging of HE 1104−1805 in the i-band. Red circles denote regions within 8″, 45″ and 120″ of the centre of the lens.
In the text

	Fig. 8. κext measurement showing that the line-of-sight of HE 1104−1805 is along a slightly underdense region of the Universe. Therefore, the fact that multiple perturbers are near the main lens is coincidental and the approximation κext ≈ κs holds.
In the text

	Fig. 7. Position and weight (size of the dot) of galaxies between us and the source. Black circles indicate the 45″ and 120″ apertures. Blue dots correspond to objects with i-magnitudes less than 24.
In the text

	Fig. 9. Reconstructed image of HE 1104−1805 using the best power-law and composite models. From left to right: The imaging data, reconstructed image, and normalised residuals are displayed for the power law model (PEMD + Shear + perturbers 2, 5, 6, 9 + groups 1, 2, 3, and 4) and composite models (NFW + Double chameleon + Shear + perturbers 2, 5, 6, 9 + groups 1, 2, 3, and 4).
In the text

	Fig. 10. HE 1104−1805 modelled convergence κ with the two families of models and different perturbers and groups inclusion strategies. The top row shows a wide-field including all groups and the bottom row zooms in on the region closer to the main lens; white lines indicate the area enclosing 25, 50, and 75% of the total mass.
In the text

Fig. 11. BIC weighted comparison of the power-law and composite models. For each family, we marginalised over the three posteriors given by the three perturber inclusion scenarios, which results in the multi-modality of some of the parameters such as γ. We obtained ${\theta }_{\mathrm{E}}^{\mathrm{PL}}=1.356\pm 0.115$, γpl = 1.88 ± 0.04, ${\gamma }_{\mathrm{ext}}^{\mathrm{PL}}=0.14\pm 0.01$, ${\theta }_{\mathrm{E}}^{\mathrm{Comp}}=1.37\pm 0.057$, ${\gamma }_{\mathrm{ext}}^{\mathrm{Comp}}=0.12\pm 0.02$, and γComp = 2.03 ± 0.08.

In the text

Fig. 12. Left panel: Kinematic fit using power-law and composite models. Red and blue central dots show the best-fit values, while error bars are the 16th and 84th percentiles of 100 predicted velocity dispersions. Right panel: Posterior of Dd and DΔt using kinematics based on the power-law and composite mass models. Using likelihood weighting we measured $D_{\mathrm{\Delta }t}={5083}_{487}^{+504}$ Mpc and $D_{\mathrm{d}}={1921}_{-236}^{+254}$ Mpc.

In the text

	Fig. 13. H0 measurement with HE 1104−1805 using the joint constraint of DΔt and Dd measured with the power-law, composite models, and both models combined. Our final measurement using the likelihood weighted combination of both models is $H_0=64.2_{-5.0}^{+5.8}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ for λint = 1.
In the text

	Fig. A.1. PSF parameters fitted across wavelength. The values obtained when using stacks of 12.5 Åare used as constraining prior for the wavelength-by-wavelength fit. The shaded area indicates the range of wavelength used for the σLOS measurement, and the dotted lines highlight the mean value of the parameters in this range.
In the text

	Fig. A.2. Top panel: three different masked MUSE cube cutouts. Bottom panel: measured velocity of the three bins in each case.
In the text

Fig. B.1. Integrated spectra from the MUSE data cube with apertures shown in Fig. 4. The binned flux is computed by imposing a minimal S/N of 5 per wavelength bin. From top to bottom, Lens centre: from left to right, the stellar absorption lines CaH, CaK, H[[INLINE557]], and G-band (orange dashed lines) are consistent with the previously determined zd = 0.729. P5: from left to right, the galactic emission lines Ne III, He I, H[[INLINE559]], and H[[INLINE560]] are highlighted with green dashed lines and allowed us to measure z5 = 0.3575. We also identify the stellar absorption lines H[[INLINE562]], HEI, CaH, and CaK allowing a velocity dispersion measurement. P6: from left to right the lines Mg I, Ca II O II, H[[INLINE563]], and O III doublet allowed us to measure z6 = 0.505. Quasar images: from left to right, we recognise the C IV, C III, C II, and Mg II emission lines.

In the text

	Fig. B.2. Measurement of the P5 stellar kinematics. Point-estimate of σLOSP5.
In the text

	Fig. C.1. From left to right: Initial guess of the PSF from stacked stars in the field, final PSF obtained, and the difference between both.
In the text

Fig. C.2. HE 1104−1805 reconstructed light component with power-law. From left to right: The light component unconvolved, convolved with the PSF, and its corresponding feature in the imaging data. The top row shows the lens light profile, and the bottom row shows the lensed source light profile. This figure shows that the light-components model corresponds to expectations based on the observation. The centre of the PSF and of the lens are not ideally modelled, which justifies the noise map increase in the first region and the masking of the second region.

In the text

Fig. C.3. HE 1104−1805 reconstructed light component with power-law. From left to right: The light component unconvolved and then convolved with the PSF and its corresponding feature in the imaging data. The top row shows the lens light profile, and the bottom row shows the lensed source light profile. This figure shows that the light-components model corresponds to expectations based on the observation. The centre of the PSF and of the lens are not ideally modelled, which justifies the noise map increase in the first region and the masking of the second region.

In the text

Fig. D.1. Galaxy groups identified in the spectroscopic sample of galaxies in the field of HE 1104−1805. For each group, the upper plot shows the positions of the candidate member galaxies associated with that group relative to the lens galaxy, with rejected group members represented as red squares and accepted group members represented as green circles. The lens galaxy (star) and group centroid (cross) are also displayed. The R200 radius of the group is represented by a solid line, while the dashed circle represents the angular separation cut of the group-finding algorithm in its final iteration. The lower plot shows the observer-frame velocity of individual member galaxies relative to the group centroid as a function of that galaxy’s angular distance from the centroid. Galaxies that passed the iterative algorithm described in Buckley-Geer et al. (2020), Sluse et al. (2017), and Sluse et al. (2019) are shown in green, while trial galaxy members that were cut through the algorithm are shown in red. Horizontal error bars represent the measurement error for each galaxy. The final observer-frame velocity dispersion and angular separation cuts from the group-finding algorithm are presented as dashed and dotted lines, respectively. We also show 1-D histograms and rug plots of the velocity and distance distributions of the member galaxies. The 1-D histograms are produced using a kernel density estimate (KDE) with a bandwidth determined using Scott’s Rule. In the 1D velocity histogram, the dashed blue line shows a Gaussian with a width equal to the observer-frame velocity dispersion of the group.

In the text