The impact of line-of-sight structures on measuring $H_0$ with strong lensing time-delays

Measurements of The Hubble-Lemaitre constant from early- and local-universe observations show a significant discrepancy. In an attempt to understand the origin of this mismatch, independent techniques to measure $H_0$ are required. One such technique, strong lensing time delays, is set to become a leading contender amongst the myriad methods due to forthcoming large strong lens samples. It is therefore critical to understand the systematic effects inherent in this method. In this paper, we quantify the influence of additional structures along the line-of-sight by adopting realistic lightcones derived from the \textit{CosmoDC2} semi-analytical extra-galactic catalogue. Using multiple lens plane ray-tracing to create a set of simulated strong lensing systems, we have investigated the impact of line-of-sight structures on time-delay measurements and in turn, on the inferred value of $H_0$. We have also tested the reliability of existing procedures for correcting for line-of-sight effects. We find that if the integrated contribution of the of line-of-sight structures is close to a uniform mass sheet, the bias in $H_0$ can be adequately corrected by including a constant external convergence $\kappa_{\rm ext}$ in the lens model. However, for realistic line-of-sight structures comprising many galaxies at different redshifts, this simple correction over-estimates the bias by a factor of approximately three. We therefore conclude that lens modelling must incorporate multiple lens planes to account for line-of-sight structures for accurate and precise inference of $H_0$.


INTRODUCTION
The Hubble-Lemaitre constant, H0, is a cornerstone of the standard cosmological model, setting the distance scale, age and critical density of the Universe. Accurate estimation of the value of H0 is therefore critical for constraining cosmological models in the era of precision cosmology. However, presently, there is a significant mismatch between H0 determined from early-and late-universe probes Verde et al. 2019), for instance, measurements of the Cosmic Microwave Background (CMB; see Bennett et al. 2013;Planck Collaboration et al. 2018) and Baryon Acoustic Oscillations (BAO; see Addison et al. 2018;DES Collaboration et al. 2020) and those made in the more local Universe using supernovae (SNe;see Dhawan et al. 2018;Macaulay et al. 2019), the tip of the red giant branch (TRGB; see Freedman et al. 2019;Yuan et al. 2019) and Cepheid variables (Riess E-mail: nan.li@nottingham.ac.uk et al. 2019; Pietrzyński et al. 2019). Independent from any of the aforementioned methods, strong lensing time delays provide valuable measurements of H0 (e.g., Wong et al. 2019;Shajib et al. 2019) which may assist in the understanding of these discrepancies once systematic uncertainties in the technique are fully calibrated. With such systematics in mind, in this paper we focus on the effects of line-of-sight structure, one of the most dominant sources of error in the lens time delay method.
Strong lensing time delays are observed when a variation in flux of a strongly-lensed background source such as a quasar, supernova or a gravitational wave event is detected at different times between its multiple images. The deflection of the light path from the source due to the gravitational potential of a lens, as well as the structures along the line-of-sight, leads to both a geometrical and a gravitational delay of the arrival time of the light from the source. The geometrical delays are sensitive to H0 (see Schneider et al. 1992). Therefore, measuring the time delays and re-constructing the mass distribution of the lens accurately allows H0 to be estimated. The existing relative paucity of strong-lens systems suitable for this method and the necessary long monitoring campaigns has somewhat limited the use of this technique but good progress has already been made with only a handful of systems (e.g., Suyu et al. 2010Suyu et al. , 2013Birrer et al. 2016;Wong et al. 2017;Bonvin et al. 2017;Wong et al. 2019;Chen et al. 2019). However, this is set to dramatically change (Oguri & Marshall 2010;Collett 2015) with the advent of the Rubin Observatory Legacy Survey of Space and Time 1 (LSST), which will give rise to about 400 well-measured time delay systems to constrain H0 to within only a few percent Dobler et al. 2015).
Even with precise time delay measurements, the reliability of estimates of H0 depends on how faithfully the lens mass model follows the true lensing mass. Degeneracies and inadequacies in the parameterisation of the lens mass model can directly propagate into the inferred value of H0 (e.g., see Schneider & Sluse 2013;Sereno & Paraficz 2014;Xu et al. 2016;Muñoz & Kamionkowski 2017;Tie & Kochanek 2018;Tagore et al. 2018;Wertz et al. 2018) as can selection effects within the lens sample (see Collett & Cunnington 2016). In addition, perturbative effects from sub-structure within the main lens and from structure along the line-of-sight can significantly modify time delays which can bias measurements of H0 if not properly taken into account. One approach to account for these effects is to directly characterise perturbing structures identified in observations (e.g., Wong et al. 2011;Momcheva et al. 2015;Rusu et al. 2017;Sluse et al. 2017;Wong et al. 2018). Another common technique is to use external shear, γext, and external convergence, κext, in the lens model. By connecting cosmological simulations and real observations, an estimate of the distribution function of the amplitude of these external lensing effects can be obtained (e.g., Suyu et al. 2010Suyu et al. , 2013Greene et al. 2013;Collett et al. 2013;Rusu et al. 2017;Birrer et al. 2017;Tihhonova et al. 2018). However, the corrections provided by γext and κext are isotropic and cannot properly capture the complexity of real perturbing structures. Motivated by this, more sophisticated approaches have been developed using multiple lens planes or approximations thereof (e.g., McCully et al. 2014;Birrer et al. 2017;McCully et al. 2017).
In this work, we investigate the influence of halos along the line-of-sight on measurements of H0 by using multiple lens plane ray-tracing simulations. To obtain simulated time delays we construct the lightcone of each lens from a stateof-the-art semi-analytic model (CosmoDC2 2 ; Korytov et al. 2019) based upon the large Outer Rim cosmological N-body simulation (Heitmann et al. 2019). By modelling these time delays with the same methods used for real data, we directly assess the biases introduced by line-of-sight effects and the efficacy with which these can be accounted for using external corrections such as γext and κext.
The paper is structured as follows. We outline the methodology used for determining strong lensing time delays in the cases of the single-lens plane and multiple-lens planes in Section 2. Details of the simulations and the process of estimating H0 from the simulated data are given in 1 https://www.lsst.org/ 2 https://portal.nersc.gov/project/lsst/cosmoDC2 Sections 3 and 4 respectively. We present our findings in Section 5, then conclude with a summary and discussion in Section 6. The cosmological model adopted in this paper is the same as used in CosmoDC2: ΛCDM with ΩΛ = 0.735, ΩM = 0.265, and h = 0.71, which is closed to the best-fit WMAP-7 parameter set (Komatsu et al. 2011).

STRONG LENSING TIME DELAYS
In this section, we present a basic description of the theory of time-delays in strong lensing systems with multiply-lensed point sources we have used in this work, for the cases of single and multiple lens planes. Throughout the paper, we have applied the thin lens approximation. For more details, we refer the reader to Schneider et al. (1992) and Narayan & Bartelmann (1996).

Time Delays in Single Lens Planes
For the case of a lensing system with a single deflector, adhering to the thin lens approximation, one can project the three-dimensional mass distribution to a two-dimensional mass sheet normal to the line-of-sight from the observer to the source. The dimensionless surface mass density of a thin lens plane can be written as a function of the lens plane angular position vector, θ, as with the critical surface mass density where Ds and D d are the angular diameter distances from the source and lens to the observer respectively, D ds is the angular diameter distance from the lens to the source, and Σ(θD d ) is the surface mass density of the lens. The lensing potential is given by and the deflection angle vector is given by Once the deflection field at the lens plane is known, we can construct the lensing equation for a given set of source planes. For example, in the case of a single lens plane and a single source plane, the lensing equation is simply where β is the angular source plane position vector that maps to θ in the image plane (or, equivalently, "lens plane" for the case of single lens-plane). Based on Eq. 5, ray-tracing simulations can be performed from the observer, crossing the lens plane to the source plane to produce lensed images. For extended source-like galaxies, to create distorted lensed images, interpolation can be used in the source plane to map spatially varying surface brightness back to the image plane. However, for the point sources used in this work, one has to adopt triangle mapping and a barycentric coordinate system to solve the lensing equations numerically. Details of the approach are discussed in Sec. 3.3.
In the case of a single lens plane, the delay of the arrival time of a light ray from the source to the observer is where z d is the redshift of the lens. The last term in Eq. 6 is also known as the Fermat potential, This delay is undetectable, the true observable being the difference between the arrival time of two separate lensed images (say, image A and image B), tAB ≡ τA − τB. From Eq. 6, the time difference can be written where, and Note that where E(z) = Ωr(1 + z) 4 + Ωm(1 + z) 3 + Ω k (1 + z) 2 + ΩΛ .
These equations show that and thus H0 can be measured from tAB if the mass distribution of the lens is reconstructed accurately.

Time Delays in Multiple Lens Planes
In the case of multiple lens planes, the lens equation must be modified to account for multiple deflections; where the quantities retain their definition from the single lens plane case but now take on a subscript referring to a specific lens plane. We consider N mass distributions, each characterised by a surface mass density Σi, at redshift zi, ordered such that zi < zj for i < j and such that the source has a redshift zs > zN . The physical distance, ξj, of the intersections on the lens planes from the optic axis (i.e., the impact parameters) are then where Di is the angular diameter distance from the observer to each lens plane, Dij (such that i < j) is the angular diameter distance from the ith lens plane to the jth lens plane andαi is the deflection angle at the ith lens plane (see Fig. 1). For simplicity, we convert the physical distance to angular positions on the sky θi = ξi/Di and the deflection angles to effective movements on the sky where Dis is the angular diameter distance from the ith lens plane to the source plane. By defining a factor Γij eq. 15 becomes In particular, for j = N + 1 = s, Γis = 1, thus, The delay of the arrival time of a deflected light path compared to a straight light path is the integral of the time difference along the line-of-sight though all lens planes. For instance, the time delay created by lens plane i and j is where the first term is the geometric delay and the second is the gravitational delay. Replacing j with i + 1 and summing over all time delays gives the total time delay through the whole line-of-sight, Therefore, similar to the case of a single lens plane, the time delay between two separate lensed images A and B can be given by which means that deflection fields, lensing potentials and the angular positions of the intersections on the lens planes are all required for the calculation of time delays in multiple lens plane systems. In section 3, we discuss how we construct a lightcone and model the lenses to obtain the information required to implement time-delay simulations with multiple lens planes.

SIMULATIONS
To quantify the influence of galaxies along the line-of-sight on measuring H0 with strong lensing time-delays, we generated simulated images following the formalism in Sec. 2 for both single and multiple lens planes with a strong lensing simulation pipeline named PICS (Li et al. 2016). In this section, we describe the simulations used and how the lens equations were solved using a triangle-mapping algorithm. Figure 1. A schematic view of the multi-plane formalism, as described in Section 2.2. A light ray (solid black line) experiences a deflection only when it passes through a lens plane (vertical solid grey lines). The deflection angleα i is the actual deflection of a ray passing through the ith lens plane,α i can be calculated from the surface density Σ i on the ith lens plane. Using the deflection anglê α i and the position of the intersection of the light ray at the (i − 1)th lens plane ξ i−1 and that on the ith lens plane ξ i , the physical position of the intersection in the (i + 1)th plane ξ i+1 can be obtained based on ξ i−1 and ξ i .

Semi-Analytic Lightcones
For creating lightcones with realistic spatial and redshift distributions of the galaxies, we extract lightcones from the CosmoDC2 (Korytov et al. 2019). Designed for an LSST data challenge project, it is established upon a large cosmological simulation called The Outer Rim Simulation run by the Argonne Cosmology Group using the Hybrid/Hardware Accelerated Cosmology Code (HACC, Habib et al. 2016).
CosmoDC2 covers 500 square degrees in the redshift range 0.0 ≤ z ≤ 3.0 and is complete to a magnitude depth of 28 in the r-band. Each galaxy is characterised by a multitude of properties including stellar mass, morphology, spectral energy distributions, broadband filter magnitudes, host halo information and weak lensing shear.
The lightcones for each of our strong lensing simulations are cut out from the full lightcone of CosmoDC2. Each extracted lightcone is centred on a bright central galaxy (BCG) identified in the cosmoDC2 since these massive central elliptical galaxies are likely strong lensing candidates. Each BCG forms the primary lens mass in its corresponding lightcone. The field of view of the lightcones is 20 × 20 , and the corresponding simulated images are 512 × 512 pixels in size. To focus on the impact of line-of-sight galaxies, we select the lightcones with the primary lens located in the range of z d = 0.5 ± 0.01 and we assume a fixed source redshift of zs = 2.0. We calculate the Einstein radius of the primary lens of each lightcone and then discard lightcones that yield Einstein radii outside the range of [1.3 , 2.4 ]. The lower limit avoids resolution issues encountered by ground-based telescopes/surveys (such as CFHT, DES, and LSST) and the upper limit discards systems which give year-like time delays. In total, we selected five hundred lightcones adhering to these criteria. Furthermore, within each lightcone, we removed any deflectors with Einstein radii larger than 0.3 to concentrate our study on the effects of secondary perturbations to the lensing potential. The substructures of the primary lens were also removed to focus on the influence of line-of-sight structures only.

Ray-tracing Simulations
For each lightcone, we run two sets of simulations for generating the lens time delays. The first set includes only a single lens plane containing the primary lens galaxy. In this set, the omitted line-of-sight halos are approximated with a constant external convergence, κext, in the lens model when computing deflection angles. For each lightcone, we estimate its value of κext by tracing multiple rays throughout it as described in more detail below. In the second set of simulations, we use one lens plane per halo for all halos, including the primary lens, in the lightcone.
In the simulations, we assume a singular isothermal ellipsoid (SIE) density profile for all halos in the lightcone, since this provides a realistic model for the total mass profile of real elliptical galaxies (Koopmans et al. 2006;Bolton et al. 2012;Shu et al. 2016). The deflection angles of the SIE model are given by Kormann et al. (1994); Keeton (2001), where φ 2 = q 2 x 2 + y 2 , q is the minor to major axis ratio and b is an effective factor to represent Einstein radius, In the case of circular lenses, b can be calculated from the velocity dispersion. The lensing potential can be computed according to the relationship between the lensing potential and the deflection field of SIE model (Keeton 2001), The complete parameter set required by equations (23 − 26) is {x1, x2, σv, q, Θ, z d }, where (x1, x2) is the angular position of the SIE profile centre with respect to the centre of the field of view, σv is the velocity dispersion of the lens, q is the ellipse axis ratio, Θ is the position angle of the ellipsoid and z d is the redshift of the deflector. The parameters x1, x2, q, Θ, z d are taken directly from the cosmoDC2 catalogue. σv is derived from the L − σ scaling relation from the bright sample of Parker et al. (2007) given by where, log 10 (L/L ) = −0.4(magr − magr ), and magr is the apparent r-band magnitude of the galaxy given by the cosmoDC2 catalogue. We adopt the assumption in More et al. (2016) that magr evolves with redshift as magr = +1.5(z − 0.1) − 20.44 (Faber et al. 2007). Sources are described by the parameter set {y1, y2, ms, zs}, where (y1, y2) is the angular position of the source with respect to the image centre, ms is the apparent r-band magnitude of the source and zs is the redshift, fixed to zs = 2. The angular positions are randomly sampled in the source plane in the vicinity of the caustic structures and only those creating quadruply-lensed images are retained.
With a fully parametrically-defined lightcone, the simulated lensed images can be produced by ray-tracing and image-finding. For our single lens-plane simulations, we determine κext in the following manner. First, we trace rays through a given lightcone from the image plane, computing the deflections caused by all halos (including the primary lens), each in their own lens plane. Along each ray, we compute an 'external halo convergence' by summing κ as given by Eq. 1 for all secondary halos excluding the primary lens halo. This external halo convergence ignores the divergence caused by voids and so we must apply a correction to obtain κext. The correction uses the results of Collett et al. (2013) who showed that κext can be obtained by subtracting the median convergence along random sight lines from the external halo convergence. The resulting κext has an uncertainty associated with it due to the scatter in the relationship between the two quantities, but negligible bias. Firing rays along random lines-of-sight in our lightcones and computing the convergence, again using Eq. 1, yields a value of κcorr = 0.048. When correcting the external halo convergence, we distribute κcorr across all lens planes according to the lensing weights (D ds D d /Ds) for each plane and subtract them separately. Figure 3 shows the probability distribution functions of the mean and median values of κext across all lightcones obtained in the manner described. We note that our peak of κext 0.1 is higher than that of previous studies, for example, peaks of 0.075 and 0.05 in Suyu et al. (2013) and McCully et al. (2017) respectively. We attribute this mainly to our selection of BCGs from cosmoDC2 and their location within more over-dense galaxy groups. Secondary effects also likely include a difference in mass models and simulated lightcones.
With κext determined, we include it in the primary lens model for the single-plane simulations and calculate maps of the deflection angle and the lensing potential. The lensing equation in Eq. 5 is used to map the image plane back to the source plane. Since the sources in this paper are point sources, we have to adopt a triangle-mapping algorithm to solve the lensing equation. This is described further in Section 3.3.
For the case of multiple lens-planes, we ray-trace through the whole lightcone in the same manner as outlined above when computing the external halo convergence, placing each halo on its own lens plane. As Eq. 20 shows, to calculate the total time delay, the deflection map and lensing potential for every lens plane must be computed. The inter- sections of the light rays traced from the image plane (given by Eq. 18) are required for the calculation of the time delay between two lens planes. These are summed over all neighbouring pairs of lens planes to obtain the total time delay according to Eq. 21. Again, for our adopted point source, we have to apply triangle mapping and barycentric interpolation to obtain the position of lensed images for a given source position on the source plane (see Section 3.3). The same image-finding process is applied to locate the intersections of the light rays between neighbouring lens planes (see Eq. 20).
Note that we have not introduced any errors to our simulated times delays as would usually be present with observed delays since we are concerned purely with the effects of line-of-sight structure in this study. However, narrow priors are used for the time delays during modelling to allow a small degree of freedom. More details are given in Section 4.

Image Finding
Since we are concerned with multiply-imaged point-like sources, e.g. AGNs or SNe, in this work, solving the lensing equation for point sources is a critical issue in the simulation. To determine the apparent positions of our point-sources, we make use of a triangle mapping technique described in Schneider et al. (1992). First, a set of Delaunay triangles is constructed from a regular grid of image plane positions which define the intersections of light rays from the source (see Fig. 2). These image plane vertices are then mapped to the source plane. Any image plane triangles which map to a triangle in the source plane containing the source position are identified. For each of these identified image-plane triangles, we compute the barycentric coordinate of the source position inside the corresponding source-plane-mapped triangle using the relation where, (xP , yP ) are the Cartesian coordinates of the point source inside its triangle of vertices (x1, y1), (x2, y2), and (x3, y3); the corresponding barycentric coordinates are (λ1, λ2, λ3). We then assume that the barycentric coordinates are conserved between the image and source planes and use them, with the vertices of the image-plane triangle to determine the position of each image of the source. For the case of multiple lens planes, the intersections between the light rays from the source and the lens planes are required for the calculation of total time-delays. Hence, we need to ascertain all the intersections. If there are N lens planes plus one source plane in the lensing system, there are N parent triangles for the triangle on the source plane. Also, we assume the barycentric coordinates of the source are conserved in the source triangle and all parent triangles. Then the intersections can be obtained. The intersections on the first lens plane (0th plane in Fig. 1) are the positions of the lensed images.

STRONG LENS MODELLING
We use the multi-purpose open source lensing package lenstronomy 3 (Birrer et al. 2015a;Birrer & Amara 2018) to measure H0 from our simulated data. In our lens modelling, we use an SIE profile since this was used for creating the simulated lens catalogue. To investigate additional bias and uncertainties in the comparison between the input and reconstructed parameters, we also try two other SIE-based mass models in the lens-modelling process. These are the SIE + γext (SG) and SIE + κext (SK) models.
The simulated data that we fit with lenstronomy are the four image positions, the three flux ratios and three time delays. For optimisation of the lens model parameters and H0, we use lenstronomy's particle swarm optimiser (PSO) (Eberhart & Kennedy 1995) since this technique performs well in lower dimensional parameter spaces such as ours (see Birrer et al. 2015b).
The simulated lensing observables are taken as inputs for the modelling process. We take very optimistic error bounds on these lensing observables since we want to explore the lens modelling rather than observational techniques that provide the input. We suppose that all input uncertainties follow normal priors centred on the values provided by the simulation. The time delays of all multiple images obtain a standard deviation of 2 days, whereas the multiple image positions acquire a 1 − σ astrometric uncertainty of 0.004" for both RA and Dec directions. The 1 − σ Gaussian error on the flux ratio between multiple images is set to 0.1.
The SIE and SIE + κext lens models both have the same eight free parameters, namely the Einstein radius, θE, ellipticity, e1,2, lens centre of mass co-ordinates, θ1,2, source position, β1,2 and H0. In the SIE + κextmodel, the parameter κext is held fixed at the value determined from the lightcone hosting the lens being modelled (see section 3.2). The SIE + γext model has the additional two parameters γext and θγ,ext for the magnitude and direction of external shear respectively. The priors on all parameters are shown in Table 1.

Parameter Prior Model constraints
Multiple image pos.
In the process of fitting of each lens system, the particlescatter in PSO is set to 1, the number of particles is 200, and the upper limit on the number of iterations is 500. These choices yield an acceptable computation time whilst still allowing a thorough exploration of the model parameter space.

RESULTS
As described in Section 3, we simulate five hundred quadruply-imaged time-delay systems with two different configurations; the first approximates LOS structure with a constant κext and the second uses the full lightcones containing halos. For each of these two simulation configurations, we use two different lens models in lenstronomy thus giving four measurements of h for each of the five hundred lenses. These four different combinations of simulation and lens model are as follows: 1) simulated SIE + κext and SIE-Only lens model (SK|S); 2) simulated SIE + κext and SIE + κext lens model (SK|SK); 3) simulated SIE + LOSgalaxies and SIE-Only lens model (SL|S); 4) simulated SIE + LOS-galaxies and SIE + γext lens model (SL|SG). Note that we do not apply the SIE + κext lens model to the simulated data generated with full lightcones (i.e. what would be the combination SL|SK) because this is identical to the application of a corrective factor of 1 − κext we investigate later (see below). Instead, we test the efficacy of including external shear in the modelling.
Not all measurements of h obtained are valid due to the accuracy of the simulations. For instance, when a source is almost coincident with the caustic in the source plane, the magnifications of the simulated lensed images are inaccurate because of the limits of grid size and linear barycentric interpolation leading to unsuccessful model fits. Hence, we discard poor fits by setting a likelihood threshold of log(L) > −1000 to give 471, 471, 307, and 279 constraints in the cases of SK|S, SK|SK, SL|S, and SL|SG respectively. Note that, by applying this threshold in likelihood, we also remove poor fits arising from large perturbations from substructures not caught by the 0.3 cut in Einstein radius.
First we consider our analysis of the simulations created with LOS structure approximated by a constant external convergence. Fig 4 shows  model and the SIE+κext model. As expected, the PDF corresponding to the modelling that uses κext peaks at the input h = 0.71. However, when external convergence is not used in the modelling, as the figure shows, values of measured h are biased high, with the PDF peak lying ∼ 10 per cent higher than the input value. Following the procedure commonly used in the literature to correct for external convergence effects (see, for example Suyu et al. 2017), we apply a correction of 1 − κext to the biased measurements of h from the SK|S configuration. The orange histogram shown in Fig 4 shows the results of this correction. Clearly, the correction in this simplified case works well, recovering a mean value of h that is almost identical to the input value.
Second, we consider our modelling of the simulations created with the full lightcones containing halos (i.e., the cases of SL|S and SL|SG). In Fig. 5, the blue and red histograms show the distribution of measured h with and without γext in the modelling. Without any correction for external convergence, h is only biased by approximately +3 per cent on average, although the scatter is larger at around 10 per cent. Furthermore, the effect of including external shear in the modelling is negligible on average. The green and purple histograms show the SL|S and SL|SG cases corrected with (1 − κext). As the figure shows, the correction biases h significantly, the peak of both distributions shifting to smaller values and leading to an underestimate of h of ∼ 8 per cent on average. This is to be expected given the anticipated size of the correction from Fig. 4 and the small 3 per cent over-prediction of the non-corrected cases. We therefore conclude that statistically, the 1 − κext correction can not be reliably used to account for real clumpy external convergence.

DISCUSSION AND CONCLUSIONS
To quantify the influence of secondary deflectors on the measurement of H0 with strong lensing time delays, we have simulated one thousand galaxy-scale strong lensing systems with quadruply-lensed variable point objects; five hundred of these were created with a primary lens and line-of-sight halos and five hundred with a primary lens plus a constant external convergence. In the simulations with line-of-sight halos, the lightcones are extracted from a semi-analytic model based on the Outer Rim large-scale cosmological simulation. The lightcones are centred on the location of central galaxies of groups of galaxies. In the simulations constructed with external convergence, we used a single lens-plane located at the redshift of the primary lens galaxy whereas in the simulations containing halos, each halo has its own lens plane.
Using an SIE mass profile for the primary lens galaxy and the halos, and an interpolative mapping method to refine the location of the lensed point source images, we generated time delay data. This time-delay data is then modelled using lenstronomy to estimate H0 using different SIE-based lens mass models and this is compared to the known input value. Our main conclusion is that incorporating constant external convergence in the modelling only works reliably if the lensed time delays are subjected to a constant external convergence. If time-delays are subjected to perturbations due to halos lying along the line-of-sight as expected in the real Universe and no correction for external convergence is made in the modelling, the inferred value of H0 is over-estimated by approximately 3 per cent on average. However, if a constant external convergence is incorporated in the lens model with a normalisation set by the mean or median convergence of the line-of-sight halos, then an over-correction of H0 occurs such that it is biased low by ∼ 8 per cent on average. This effect can not be ignored, since the uncertainties of current measurements of H0 from strong lensing time delays are typically quoted as being lower than this (Bonvin et al. 2017;Chen et al. 2019;Wong et al. 2019;Rusu et al. 2019;Birrer et al. 2019). With the forthcoming large sample of strong lensing time delay systems observed by the future time domain large scale surveys, e.g., Mephisto 4 and LSST, the effect becomes even more problematic.
Qualitatively, our conclusions are consistent with those of McCully et al. (2017) in the sense that line-of-sight structures significantly affect the accuracy of the measurement of H0. Quantitatively, we find a larger external convergence of κext 0.1 compared to the value of 0.05 from McCully et al. We attribute this to the fact that we have selected BCGs as the primary lensing galaxies in our lightcones and because we have included more line-of-sight structures; we include galaxies from cosmoDC2 down to an r-band apparent magnitude of 28, compared to the i-band limit of 21.5 adopted by McCully et al.
We have also investigated the effects of incorporating external shear in the lens model. In the simulations using line-of-sight halos, adding an external shear term to the SIE lens model makes a negligible impact on the distribution of recovered values of h. Not unexpectedly, we also find that correcting this SIE+γext model with the average constant external convergence also leads to a ∼ 8 per cent underestimation, which implies that the influence of external shear is negligible in the case of our study. This conclusion differs from that of McCully et al., most likely because we cleaned our lens sample by removing secondary halos that give rise to an Einstein radius of greater than 0.3 arcsec.
The Outer Rim simulations used to populate our lensing lightcones with halos include only dark matter. As such, we have used SIE profiles in place of identified halos to better represent the total mass (baryons + dark matter) profiles of real lens galaxies. One effect this may have is that the lensing strength of any lower mass halos, which in the real Universe may not have accrued baryons, could be artificially enhanced by the more efficient isothermal profile. In addition, our simulated datasets do not include any large scale structure such as filaments although this is expected to be a small effect. Finally, we have ignored the effects of environmental structure in the simulations in the sense that our assumed smooth SIE profiles for the primary lens do not include substructure. We will leave consideration of these additional effects for future work.
To summarise, simple corrections for line-of-sight structure such as external shear or external convergence in estimations of H0 using lensed time delays can not be relied upon. Time delay studies opt for lens systems which are apparently free of strong perturbers in an attempt to exclude line-of-sight effects. Nevertheless, our simulations have mimicked this selection by removing any halo from all of our lightcones that produces a deflection with an Einstein radius larger than 0.3". Our work reveals that more sophisticated lens modelling methods, likely including multiple lens planes, are required to reliably measure H0 from the hundreds of well-measured time-delay systems anticipated in forthcoming large strong lens samples.