The impact of line-of-sight structures on measuring H 0 with strong lensing time-delays
MMNRAS , 1–9 (20**) Preprint 24 June 2020 Compiled using MNRAS L A TEX style file v3.0
The impact of line-of-sight structures on measuring H with strong lensing time-delays Nan Li , (cid:63) , Christoph Becker , Simon Dye School of Physics and Astronomy, University of Nottingham, University Park, Nottingham, NG7 2RD, UK. National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Road, Beijing 100012, China. Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK.
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Measurements of The Hubble-Lemaitre constant from early- and local-universe obser-vations show a significant discrepancy. In an attempt to understand the origin of thismismatch, independent techniques to measure H are required. One such technique,strong lensing time delays, is set to become a leading contender amongst the myriadmethods due to forthcoming large strong lens samples. It is therefore critical to un-derstand the systematic effects inherent in this method. In this paper, we quantify theinfluence of additional structures along the line-of-sight by adopting realistic lightconesderived from the CosmoDC2 semi-analytical extra-galactic catalogue. Using multiplelens plane ray-tracing to create a set of simulated strong lensing systems, we haveinvestigated the impact of line-of-sight structures on time-delay measurements andin turn, on the inferred value of H . We have also tested the reliability of existingprocedures for correcting for line-of-sight effects. We find that if the integrated contri-bution of the of line-of-sight structures is close to a uniform mass sheet, the bias in H can be adequately corrected by including a constant external convergence κ ext in thelens model. However, for realistic line-of-sight structures comprising many galaxies atdifferent redshifts, this simple correction over-estimates the bias by a factor of approx-imately three. We therefore conclude that lens modelling must incorporate multiplelens planes to account for line-of-sight structures for accurate and precise inference of H . Key words: gravitational lensing: strong – cosmology: cosmological parameters
The Hubble-Lemaitre constant, H , is a cornerstone of thestandard cosmological model, setting the distance scale, ageand critical density of the Universe. Accurate estimation ofthe value of H is therefore critical for constraining cosmo-logical models in the era of precision cosmology. However,presently, there is a significant mismatch between H de-termined from early- and late-universe probes (Riess 2019;Verde et al. 2019), for instance, measurements of the Cos-mic Microwave Background (CMB; see Bennett et al. 2013;Planck Collaboration et al. 2018) and Baryon Acoustic Os-cillations (BAO; see Addison et al. 2018; DES Collaborationet al. 2020) and those made in the more local Universe us-ing supernovae (SNe;see Dhawan et al. 2018; Macaulay et al.2019), the tip of the red giant branch (TRGB; see Freedmanet al. 2019; Yuan et al. 2019) and Cepheid variables (Riess (cid:63) E-mail: [email protected] et al. 2019; Pietrzy´nski et al. 2019). Independent from any ofthe aforementioned methods, strong lensing time delays pro-vide valuable measurements of H (e.g., Wong et al. 2019;Shajib et al. 2019) which may assist in the understanding ofthese discrepancies once systematic uncertainties in the tech-nique are fully calibrated. With such systematics in mind, inthis paper we focus on the effects of line-of-sight structure,one of the most dominant sources of error in the lens timedelay method.Strong lensing time delays are observed when a vari-ation in flux of a strongly-lensed background source suchas a quasar, supernova or a gravitational wave event is de-tected at different times between its multiple images. Thedeflection of the light path from the source due to the grav-itational potential of a lens, as well as the structures alongthe line-of-sight, leads to both a geometrical and a gravita-tional delay of the arrival time of the light from the source.The geometrical delays are sensitive to H (see Schneideret al. 1992). Therefore, measuring the time delays and re- c (cid:13) a r X i v : . [ a s t r o - ph . C O ] J un Nan Li et al. constructing the mass distribution of the lens accurately al-lows H to be estimated. The existing relative paucity ofstrong-lens systems suitable for this method and the neces-sary long monitoring campaigns has somewhat limited theuse of this technique but good progress has already beenmade with only a handful of systems (e.g., Suyu et al. 2010,2013; Birrer et al. 2016; Wong et al. 2017; Bonvin et al. 2017;Wong et al. 2019; Chen et al. 2019). However, this is set todramatically change (Oguri & Marshall 2010; Collett 2015)with the advent of the Rubin Observatory Legacy Survey ofSpace and Time (LSST), which will give rise to about 400well-measured time delay systems to constrain H to withinonly a few percent (Liao et al. 2015; Dobler et al. 2015).Even with precise time delay measurements, the relia-bility of estimates of H depends on how faithfully the lensmass model follows the true lensing mass. Degeneracies andinadequacies in the parameterisation of the lens mass modelcan directly propagate into the inferred value of H (e.g., seeSchneider & Sluse 2013; Sereno & Paraficz 2014; Xu et al.2016; Mu˜noz & Kamionkowski 2017; Tie & Kochanek 2018;Tagore et al. 2018; Wertz et al. 2018) as can selection effectswithin the lens sample (see Collett & Cunnington 2016). Inaddition, perturbative effects from sub-structure within themain lens and from structure along the line-of-sight can sig-nificantly modify time delays which can bias measurementsof H if not properly taken into account. One approach toaccount for these effects is to directly characterise perturb-ing 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 useexternal shear, γ ext , and external convergence, κ ext , in thelens model. By connecting cosmological simulations and realobservations, an estimate of the distribution function of theamplitude of these external lensing effects can be obtained(e.g., Suyu et al. 2010, 2013; Greene 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 ofreal perturbing structures. Motivated by this, more sophis-ticated approaches have been developed using multiple lensplanes 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 alongthe line-of-sight on measurements of H by using multiplelens plane ray-tracing simulations. To obtain simulated timedelays we construct the lightcone of each lens from a state-of-the-art semi-analytic model (CosmoDC2 ; Korytov et al.2019) based upon the large Outer Rim cosmological N-bodysimulation (Heitmann et al. 2019). By modelling these timedelays with the same methods used for real data, we directlyassess the biases introduced by line-of-sight effects and theefficacy with which these can be accounted for using externalcorrections such as γ ext and κ ext .The paper is structured as follows. We outline themethodology used for determining strong lensing time de-lays in the cases of the single-lens plane and multiple-lensplanes in Section 2. Details of the simulations and the pro-cess of estimating H from the simulated data are given in https://portal.nersc.gov/project/lsst/cosmoDC2 Sections 3 and 4 respectively. We present our findings inSection 5, then conclude with a summary and discussion inSection 6. The cosmological model adopted in this paper isthe same as used in CosmoDC2: ΛCDM with Ω Λ = 0 . M = 0 . h = 0 .
71, which is closed to the best-fitWMAP-7 parameter set (Komatsu et al. 2011).
In this section, we present a basic description of the theory oftime-delays in strong lensing systems with multiply-lensedpoint sources we have used in this work, for the cases ofsingle and multiple lens planes. Throughout the paper, wehave applied the thin lens approximation. For more details,we refer the reader to Schneider et al. (1992) and Narayan& Bartelmann (1996).
For the case of a lensing system with a single deflector, ad-hering to the thin lens approximation, one can project thethree-dimensional mass distribution to a two-dimensionalmass sheet normal to the line-of-sight from the observer tothe source. The dimensionless surface mass density of a thinlens plane can be written as a function of the lens planeangular position vector, θ , as κ ( θ ) = Σ( θ D d ) / Σ crit , (1)with the critical surface mass densityΣ crit = c πG D s D d D ds , (2)where D s and D d are the angular diameter distances fromthe source and lens to the observer respectively, D ds is theangular diameter distance from the lens to the source, andΣ( θ D d ) is the surface mass density of the lens. The lensingpotential is given by ψ ( θ ) = 1 π (cid:90) d θ (cid:48) κ ( θ (cid:48) )ln | θ − θ (cid:48) | , (3)and the deflection angle vector is given by α ( θ ) = 1 π (cid:90) d θ (cid:48) κ ( θ (cid:48) ) θ − θ (cid:48) | θ − θ (cid:48) | . (4)Once the deflection field at the lens plane is known, wecan construct the lensing equation for a given set of sourceplanes. For example, in the case of a single lens plane and asingle source plane, the lensing equation is simply β = θ − α ( θ ) , (5)where β is the angular source plane position vector thatmaps to θ in the image plane (or, equivalently, “lens plane”for the case of single lens-plane). Based on Eq. 5, ray-tracingsimulations can be performed from the observer, crossingthe lens plane to the source plane to produce lensed images.For extended source-like galaxies, to create distorted lensedimages, interpolation can be used in the source plane tomap spatially varying surface brightness back to the imageplane. However, for the point sources used in this work, onehas to adopt triangle mapping and a barycentric coordinate MNRAS , 1–9 (20**) with L.O.S. halos system to solve the lensing equations numerically. Details ofthe approach are discussed in Sec. 3.3.In the case of a single lens plane, the delay of the arrivaltime of a light ray from the source to the observer is τ ( θ , β ) = (1 + z d ) c D d D s D ds (cid:20) ( θ − β ) − ψ ( θ ) (cid:21) , (6)where z d is the redshift of the lens. The last term in Eq. 6is also known as the Fermat potential,Φ( θ , β ) ≡ (cid:20) ( θ − β ) − ψ ( θ ) (cid:21) . (7)This delay is undetectable, the true observable being thedifference between the arrival time of two separate lensedimages (say, image A and image B), t AB ≡ τ A − τ B . FromEq. 6, the time difference can be written t AB = D ∆ τ c ∆Φ AB , (8)where, D ∆ τ ≡ (1 + z d ) D d D s D ds (9)and∆Φ ≡ Φ( θ A , β ) − Φ( θ B , β ) . (10)Note that D a ( z ) = cH (1 + z ) (cid:90) z dz (cid:48) E ( z (cid:48) ) (11)where E ( z ) = (cid:112) Ω r (1 + z ) + Ω m (1 + z ) + Ω k (1 + z ) + Ω Λ . (12)These equations show that t AB ∝ D ∆ τ ∝ H (13)and thus H can be measured from t AB if the mass distri-bution of the lens is reconstructed accurately. In the case of multiple lens planes, the lens equation mustbe modified to account for multiple deflections; β = θ − N (cid:88) i =1 α i ( θ i ) , (14)where the quantities retain their definition from the singlelens plane case but now take on a subscript referring to aspecific lens plane. We consider N mass distributions, eachcharacterised by a surface mass density Σ i , at redshift z i ,ordered such that z i < z j for i < j and such that the sourcehas a redshift z s > z N . The physical distance, ξ j , of theintersections on the lens planes from the optic axis (i.e., theimpact parameters) are then ξ j = D j D ξ − j − (cid:88) i =1 D ij ˆ α i ( ξ i ) , (15)where D i is the angular diameter distance from the observerto each lens plane, D ij (such that i < j ) is the angular diameter distance from the i th lens plane to the j th lensplane and ˆ α i is the deflection angle at the i th lens plane(see Fig. 1). For simplicity, we convert the physical distanceto angular positions on the sky θ i = ξ i /D i and the deflectionangles to effective movements on the sky α i = D is D s ˆ α i , (16)where D is is the angular diameter distance from the i th lensplane to the source plane. By defining a factor Γ ij Γ ij = D ij D s D j D is , (17)eq. 15 becomes θ j = θ − j − (cid:88) i =1 Γ ij α i ( θ i ) . (18)In particular, for j = N + 1 = s , Γ is = 1, thus, β ≡ θ N +1 = θ − N (cid:88) i =1 α i ( θ i ) . (19)The delay of the arrival time of a deflected light pathcompared to a straight light path is the integral of the timedifference along the line-of-sight though all lens planes. Forinstance, the time delay created by lens plane i and j is τ ij ( θ i , θ j ) = 1 + z i c D i D j D ij (cid:20)
12 ( θ i − θ j ) − Γ ij ψ ( θ i ) (cid:21) , (20)where the first term is the geometric delay and the second isthe gravitational delay. Replacing j with i + 1 and summingover all time delays gives the total time delay through thewhole line-of-sight, τ ( θ , ..., θ N , β ) = N (cid:88) i =1 τ i,i +1 ( θ i , θ i +1 ) . (21)Therefore, similar to the case of a single lens plane, the timedelay between two separate lensed images A and B can begiven by t AB ≡ τ A − τ B = N (cid:88) i =1 τ i,i +1 ( θ A ,i , θ A ,i +1 ) − N (cid:88) i =1 τ i,i +1 ( θ B ,i , θ B ,i +1 ) , (22)which means that deflection fields, lensing potentials andthe angular positions of the intersections on the lens planesare all required for the calculation of time delays in multiplelens plane systems. In section 3, we discuss how we constructa lightcone and model the lenses to obtain the informationrequired to implement time-delay simulations with multiplelens planes. To quantify the influence of galaxies along the line-of-sighton measuring H with strong lensing time-delays, we gen-erated simulated images following the formalism in Sec. 2for both single and multiple lens planes with a strong lens-ing simulation pipeline named PICS (Li et al. 2016). In thissection, we describe the simulations used and how the lensequations were solved using a triangle-mapping algorithm.
MNRAS000
MNRAS000 , 1–9 (20**)
Nan Li et al. (cid:18) (cid:4) (cid:18) (cid:4) (cid:18) i-1 (i-1)i i(i+1)1 Figure 1.
A schematic view of the multi-plane formalism, as described in Section 2.2. A light ray (solid black line) experiences adeflection only when it passes through a lens plane (vertical solid grey lines). The deflection angle ˆ α i is the actual deflection of a raypassing through the i th lens plane, ˆ α i can be calculated from the surface density Σ i on the i th lens plane. Using the deflection angleˆ α i and the position of the intersection of the light ray at the ( i − ξ i − and that on the i th lens plane ξ i , the physicalposition of the intersection in the ( i + 1)th plane ξ i +1 can be obtained based on ξ i − and ξ i . For creating lightcones with realistic spatial and redshift dis-tributions of the galaxies, we extract lightcones from theCosmoDC2 (Korytov et al. 2019). Designed for an LSSTdata challenge project, it is established upon a large cosmo-logical simulation called The Outer Rim Simulation run bythe Argonne Cosmology Group using the Hybrid/HardwareAccelerated Cosmology Code (HACC, Habib et al. 2016).CosmoDC2 covers 500 square degrees in the redshift range0 . ≤ z ≤ . (cid:48)(cid:48) × (cid:48)(cid:48) , and thecorresponding simulated images are 512 ×
512 pixels in size.To focus on the impact of line-of-sight galaxies, we selectthe lightcones with the primary lens located in the range of z d = 0 . ± .
01 and we assume a fixed source redshift of z s = 2 .
0. We calculate the Einstein radius of the primarylens of each lightcone and then discard lightcones that yieldEinstein radii outside the range of [1 . (cid:48)(cid:48) , . (cid:48)(cid:48) ]. The lowerlimit avoids resolution issues encountered by ground-basedtelescopes/surveys (such as CFHT, DES, and LSST) andthe upper limit discards systems which give year-like timedelays. In total, we selected five hundred lightcones adher-ing to these criteria. Furthermore, within each lightcone, weremoved any deflectors with Einstein radii larger than 0 . (cid:48)(cid:48) to concentrate our study on the effects of secondary pertur-bations to the lensing potential. The substructures of the primary lens were also removed to focus on the influence ofline-of-sight structures only. For each lightcone, we run two sets of simulations for gen-erating the lens time delays. The first set includes only asingle lens plane containing the primary lens galaxy. In thisset, the omitted line-of-sight halos are approximated with aconstant external convergence, κ ext , in the lens model whencomputing deflection angles. For each lightcone, we estimateits value of κ ext by tracing multiple rays throughout it asdescribed in more detail below. In the second set of simula-tions, we use one lens plane per halo for all halos, includingthe primary lens, in the lightcone.In the simulations, we assume a singular isothermal el-lipsoid (SIE) density profile for all halos in the lightcone,since this provides a realistic model for the total mass pro-file of real elliptical galaxies (Koopmans et al. 2006; Boltonet al. 2012; Shu et al. 2016). The deflection angles of the SIEmodel are given by Kormann et al. (1994); Keeton (2001), α x ≡ ψ x = bq (cid:112) (1 − q ) tan − (cid:34) (cid:112) − q θ x φ (cid:35) , (23) α y ≡ ψ y = bq (cid:112) (1 − q ) tanh − (cid:34) (cid:112) − q θ y φ (cid:35) , (24)where φ = q x + y , q is the minor to major axis ratio and b is an effective factor to represent Einstein radius, b = 4 π √ q (cid:16) σc (cid:17) D ls D s . (25)In the case of circular lenses, b can be calculated from thevelocity dispersion. The lensing potential can be computedaccording to the relationship between the lensing potential MNRAS , 1–9 (20**) with L.O.S. halos RaytracingInterpolation
Figure 2.
The Interpolation scheme used for determining image positions of point sources. The regular grid of rays in the image plane(left filled circles) is used to partition the image plane into triangles (grey lines in the left panel). The image positions (the open whitecircle in the left panel) of a source inside a triangle (the grey triangle in the right panel) formed by the backtraced rays on the sourceplane (grey filled circles in the right panel) is then determined by using linear interpolation in the barycentric coordinates. and the deflection field of SIE model (Keeton 2001), ψ ( θ x , θ y ) = θ x ψ x + θ y ψ y . (26)The complete parameter set required by equations (23 −
26) is { x , x , σ v , q, Θ , z d } , where ( x , x ) is the angularposition of the SIE profile centre with respect to the centreof the field of view, σ v is the velocity dispersion of the lens, q is the ellipse axis ratio, Θ is the position angle of the ellip-soid and z d is the redshift of the deflector. The parameters x , x , q, Θ , z d are taken directly from the cosmoDC2 cata-logue. σ v is derived from the L − σ scaling relation from thebright sample of Parker et al. (2007) given by σ v = 142 (cid:18) LL (cid:63) (cid:19) (1 / km s − , (27)where, log ( L/L (cid:63) ) = − . mag r − mag r(cid:63) ), and mag r isthe apparent r -band magnitude of the galaxy given by thecosmoDC2 catalogue. We adopt the assumption in Moreet al. (2016) that mag r(cid:63) evolves with redshift as mag r(cid:63) =+1 . z − . − .
44 (Faber et al. 2007).Sources are described by the parameter set { y , y , m s , z s } , where ( y , y ) is the angular positionof the source with respect to the image centre, m s isthe apparent r -band magnitude of the source and z s isthe redshift, fixed to z s = 2. The angular positions arerandomly sampled in the source plane in the vicinity of thecaustic structures and only those creating quadruply-lensedimages are retained.With a fully parametrically-defined lightcone, the sim-ulated lensed images can be produced by ray-tracing andimage-finding. For our single lens-plane simulations, we de-termine κ ext in the following manner. First, we trace raysthrough a given lightcone from the image plane, comput-ing the deflections caused by all halos (including the pri-mary 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 pri-mary lens halo. This external halo convergence ignores the divergence caused by voids and so we must apply a correc-tion to obtain κ ext . The correction uses the results of Collettet al. (2013) who showed that κ ext can be obtained by sub-tracting the median convergence along random sight linesfrom the external halo convergence. The resulting κ ext hasan uncertainty associated with it due to the scatter in therelationship between the two quantities, but negligible bias.Firing rays along random lines-of-sight in our lightcones andcomputing the convergence, again using Eq. 1, yields a valueof κ corr = 0 . κ corr across all lens planes according tothe lensing weights ( D ds D d /D s ) for each plane and subtractthem separately.Figure 3 shows the probability distribution functions ofthe mean and median values of κ ext across all lightcones ob-tained in the manner described. We note that our peak of κ ext (cid:39) . .
075 and 0 .
05 in Suyu et al. (2013) andMcCully et al. (2017) respectively. We attribute this mainlyto our selection of BCGs from cosmoDC2 and their loca-tion within more over-dense galaxy groups. Secondary effectsalso likely include a difference in mass models and simulatedlightcones.With κ ext determined, we include it in the primary lensmodel for the single-plane simulations and calculate maps ofthe deflection angle and the lensing potential. The lensingequation in Eq. 5 is used to map the image plane back tothe source plane. Since the sources in this paper are pointsources, we have to adopt a triangle-mapping algorithm tosolve the lensing equation. This is described further in Sec-tion 3.3.For the case of multiple lens-planes, we ray-tracethrough the whole lightcone in the same manner as out-lined above when computing the external halo convergence,placing each halo on its own lens plane. As Eq. 20 shows, tocalculate the total time delay, the deflection map and lensingpotential for every lens plane must be computed. The inter- MNRAS000
05 in Suyu et al. (2013) andMcCully et al. (2017) respectively. We attribute this mainlyto our selection of BCGs from cosmoDC2 and their loca-tion within more over-dense galaxy groups. Secondary effectsalso likely include a difference in mass models and simulatedlightcones.With κ ext determined, we include it in the primary lensmodel for the single-plane simulations and calculate maps ofthe deflection angle and the lensing potential. The lensingequation in Eq. 5 is used to map the image plane back tothe source plane. Since the sources in this paper are pointsources, we have to adopt a triangle-mapping algorithm tosolve the lensing equation. This is described further in Sec-tion 3.3.For the case of multiple lens-planes, we ray-tracethrough the whole lightcone in the same manner as out-lined above when computing the external halo convergence,placing each halo on its own lens plane. As Eq. 20 shows, tocalculate the total time delay, the deflection map and lensingpotential for every lens plane must be computed. The inter- MNRAS000 , 1–9 (20**)
Nan Li et al. ext P D F Mean ext
Median ext
Figure 3.
The probability distribution function (PDF) of themean and median values of the fully ray-traced convergence mapof each lens lightcone used in this work. X-axis is the κ ext , y-axisis the PDF of κ ext . Blue and red histograms stand for the distri-bution of the mean and median values of the external convergencemaps separately. Blue and red curves present the correspondingkernel kernel density estimates (KDE). sections of the light rays traced from the image plane (givenby Eq. 18) are required for the calculation of the time delaybetween two lens planes. These are summed over all neigh-bouring pairs of lens planes to obtain the total time delayaccording to Eq. 21. Again, for our adopted point source,we have to apply triangle mapping and barycentric inter-polation to obtain the position of lensed images for a givensource position on the source plane (see Section 3.3). Thesame image-finding process is applied to locate the intersec-tions of the light rays between neighbouring lens planes (seeEq. 20).Note that we have not introduced any errors to oursimulated times delays as would usually be present with ob-served delays since we are concerned purely with the effectsof line-of-sight structure in this study. However, narrow pri-ors are used for the time delays during modelling to allow asmall degree of freedom. More details are given in Section 4. Since we are concerned with multiply-imaged point-likesources, e.g. AGNs or SNe, in this work, solving the lensingequation 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 inSchneider et al. (1992). First, a set of Delaunay trianglesis constructed from a regular grid of image plane positionswhich define the intersections of light rays from the source(see Fig. 2). These image plane vertices are then mapped tothe source plane. Any image plane triangles which map toa triangle in the source plane containing the source positionare identified. For each of these identified image-plane tri-angles, we compute the barycentric coordinate of the sourceposition inside the corresponding source-plane-mapped tri-angle using the relation x x x y y y λ λ λ = x P y P (28) where, ( x P , y P ) are the Cartesian coordinates of the pointsource inside its triangle of vertices ( x , y ), ( x , y ),and ( x , y ); the corresponding barycentric coordinatesare ( λ , λ , λ ). We then assume that the barycentric co-ordinates are conserved between the image and source planesand use them, with the vertices of the image-plane triangleto determine the position of each image of the source.For the case of multiple lens planes, the intersectionsbetween the light rays from the source and the lens planesare required for the calculation of total time-delays. Hence,we need to ascertain all the intersections. If there are N lensplanes plus one source plane in the lensing system, thereare N parent triangles for the triangle on the source plane.Also, we assume the barycentric coordinates of the sourceare conserved in the source triangle and all parent triangles.Then the intersections can be obtained. The intersections onthe first lens plane (0th plane in Fig. 1) are the positions ofthe lensed images. We use the multi-purpose open source lensing package lenstronomy (Birrer et al. 2015a; Birrer & Amara 2018)to measure H from our simulated data. In our lens mod-elling, we use an SIE profile since this was used for creatingthe simulated lens catalogue. To investigate additional biasand uncertainties in the comparison between the input andreconstructed parameters, we also try two other SIE -basedmass 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 arethe four image positions, the three flux ratios and three timedelays. For optimisation of the lens model parameters and H , we use lenstronomy ’s particle swarm optimiser (PSO)(Eberhart & Kennedy 1995) since this technique performswell in lower dimensional parameter spaces such as ours (seeBirrer et al. 2015b).The simulated lensing observables are taken as inputsfor the modelling process. We take very optimistic errorbounds on these lensing observables since we want to explorethe lens modelling rather than observational techniques thatprovide the input. We suppose that all input uncertaintiesfollow normal priors centred on the values provided by thesimulation. The time delays of all multiple images obtain astandard deviation of 2 days, whereas the multiple imagepositions acquire a 1 − σ astrometric uncertainty of 0 . − σ Gaussian erroron the flux ratio between multiple images is set to 0 . κ ext lens models both have the sameeight free parameters, namely the Einstein radius, θ E , ellip-ticity, e , , lens centre of mass co-ordinates, θ , , source po-sition, β , and H . In the SIE + κ ext model, the parameter κ ext is held fixed at the value determined from the light-cone hosting the lens being modelled (see section 3.2). TheSIE + γ ext model has the additional two parameters γ ext and θ γ, ext for the magnitude and direction of external shearrespectively. The priors on all parameters are shown in Ta-ble 1. https://github.com/sibirrer/lenstronomy MNRAS , 1–9 (20**) with L.O.S. halos Parameter PriorModel constraints
Multiple image pos. RA, DEC (arcsec) N ( θ sim , . F − , − , − N ( F sim , . t − , − , − (days) N (∆ t sim , Model component
Lens, SIE θ E (arcsec) U (0 . , e , U ( − . , . θ , (arcsec) U ( − , γ ext U (0 . , . θ γ, ext (rad) U ( − π, π )Source, Point β , (arcsec) U ( − , H (km/s/Mpc) U (20 , Table 1.
Priors applied to parameters in the lens modelling.These are either uniform ( U ) with upper and lower limits quotedor normal ( N ) with mean and standard deviation quoted. In the process of fitting of each lens system, the particle-scatter in PSO is set to 1, the number of particles is 200,and the upper limit on the number of iterations is 500. Thesechoices yield an acceptable computation time whilst still al-lowing a thorough exploration of the model parameter space.
As described in Section 3, we simulate five hundredquadruply-imaged time-delay systems with two differentconfigurations; the first approximates LOS structure witha constant κ ext and the second uses the full lightcones con-taining halos. For each of these two simulation configura-tions, we use two different lens models in lenstronomy thus giving four measurements of h for each of the five hun-dred lenses. These four different combinations of simulationand 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 + LOS - galaxies 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 tothe simulated data generated with full lightcones (i.e. whatwould be the combination SL | SK ) because this is identicalto the application of a corrective factor of 1 − κ ext we in-vestigate later (see below). Instead, we test the efficacy ofincluding external shear in the modelling.Not all measurements of h obtained are valid due to theaccuracy of the simulations. For instance, when a source isalmost coincident with the caustic in the source plane, themagnifications of the simulated lensed images are inaccu-rate because of the limits of grid size and linear barycen-tric interpolation leading to unsuccessful model fits. Hence,we discard poor fits by setting a likelihood threshold oflog( L ) > − SK | S , SK | SK , SL | S , and SL | SG respectively.Note that, by applying this threshold in likelihood, we alsoremove poor fits arising from large perturbations from sub-structures not caught by the 0 . (cid:48)(cid:48) cut in Einstein radius.First we consider our analysis of the simulations cre-ated with LOS structure approximated by a constant ex-ternal convergence. Fig 4 shows the PDFs of h obtainedfor the two different lens models applied, i.e. the SIE-only h P D F SK|SSK|SKSK|S|K
Figure 4.
The distribution of h determined by modelling thesimulated strong lensing data generated with constant externalconvergence. The black vertical line shows the input h used inthe simulations. The blue histogram shows the distribution of h obtained by modelling without external convergence ( SK | S )whereas the orange histogram corresponds to modelling with ex-ternal convergence ( SK | SK ). The green histogram shows the dis-tribution of h obtained after correcting the SK | S case by thefactor (1 − κ ext ) (see main text for details). model and the SIE+ κ ext model. As expected, the PDF corre-sponding to the modelling that uses κ ext peaks at the input h = 0 .
71. However, when external convergence is not used inthe modelling, as the figure shows, values of measured h arebiased high, with the PDF peak lying ∼
10 per cent higherthan the input value.Following the procedure commonly used in the litera-ture to correct for external convergence effects (see, for ex-ample Suyu et al. 2017), we apply a correction of 1 − κ ext to the biased measurements of h from the SK | S configura-tion. The orange histogram shown in Fig 4 shows the resultsof this correction. Clearly, the correction in this simplifiedcase works well, recovering a mean value of h that is almostidentical to the input value.Second, we consider our modelling of the simulationscreated with the full lightcones containing halos (i.e., thecases of SL | S and SL | SG ). In Fig. 5, the blue and redhistograms show the distribution of measured h with andwithout γ ext in the modelling. Without any correction forexternal convergence, h is only biased by approximately +3per cent on average, although the scatter is larger at around10 per cent. Furthermore, the effect of including externalshear in the modelling is negligible on average. The greenand purple histograms show the SL | S and SL | SG cases cor-rected with (1 − κ ext ). As the figure shows, the correctionbiases h significantly, the peak of both distributions shiftingto smaller values and leading to an underestimate of h of ∼ − κ ext correctioncan not be reliably used to account for real clumpy externalconvergence. MNRAS000
10 per cent higherthan the input value.Following the procedure commonly used in the litera-ture to correct for external convergence effects (see, for ex-ample Suyu et al. 2017), we apply a correction of 1 − κ ext to the biased measurements of h from the SK | S configura-tion. The orange histogram shown in Fig 4 shows the resultsof this correction. Clearly, the correction in this simplifiedcase works well, recovering a mean value of h that is almostidentical to the input value.Second, we consider our modelling of the simulationscreated with the full lightcones containing halos (i.e., thecases of SL | S and SL | SG ). In Fig. 5, the blue and redhistograms show the distribution of measured h with andwithout γ ext in the modelling. Without any correction forexternal convergence, h is only biased by approximately +3per cent on average, although the scatter is larger at around10 per cent. Furthermore, the effect of including externalshear in the modelling is negligible on average. The greenand purple histograms show the SL | S and SL | SG cases cor-rected with (1 − κ ext ). As the figure shows, the correctionbiases h significantly, the peak of both distributions shiftingto smaller values and leading to an underestimate of h of ∼ − κ ext correctioncan not be reliably used to account for real clumpy externalconvergence. MNRAS000 , 1–9 (20**)
Nan Li et al. h P D F SL | SSL | SGSL | S | KSL | SG | K Figure 5.
The distribution of h determined by modelling thesimulated strong lensing data generated from the full lightcones.Red and blue histograms show the distribution of h obtained bymodelling with SIE only ( SL | S ) and SIE+shear ( SL | SG ) respec-tively. The green and purple histograms show the distribution of h obtained after correcting the cases of ( SL | S ) and ( SL | SG ) bythe factor (1 − κ ext ) (see main text for details). To quantify the influence of secondary deflectors on the mea-surement of H with strong lensing time delays, we havesimulated one thousand galaxy-scale strong lensing systemswith quadruply-lensed variable point objects; five hundred ofthese were created with a primary lens and line-of-sight ha-los and five hundred with a primary lens plus a constant ex-ternal convergence. In the simulations with line-of-sight ha-los, the lightcones are extracted from a semi-analytic modelbased on the Outer Rim large-scale cosmological simulation.The lightcones are centred on the location of central galax-ies of groups of galaxies. In the simulations constructed withexternal convergence, we used a single lens-plane located atthe redshift of the primary lens galaxy whereas in the sim-ulations containing halos, each halo has its own lens plane.Using an SIE mass profile for the primary lens galaxy andthe halos, and an interpolative mapping method to refinethe location of the lensed point source images, we generatedtime delay data. This time-delay data is then modelled us-ing lenstronomy to estimate H using different SIE-basedlens mass models and this is compared to the known inputvalue.Our main conclusion is that incorporating constant ex-ternal convergence in the modelling only works reliably if thelensed time delays are subjected to a constant external con-vergence. If time-delays are subjected to perturbations dueto halos lying along the line-of-sight as expected in the realUniverse and no correction for external convergence is madein the modelling, the inferred value of H is over-estimatedby approximately 3 per cent on average. However, if a con-stant external convergence is incorporated in the lens modelwith a normalisation set by the mean or median convergenceof the line-of-sight halos, then an over-correction of H oc-curs such that it is biased low by ∼ H from strong lensing time de-lays are typically quoted as being lower than this (Bonvinet 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 futuretime domain large scale surveys, e.g., Mephisto and LSST,the effect becomes even more problematic.Qualitatively, our conclusions are consistent with thoseof McCully et al. (2017) in the sense that line-of-sight struc-tures significantly affect the accuracy of the measurement of H . Quantitatively, we find a larger external convergence of κ ext (cid:39) . r -band apparent mag-nitude of 28, compared to the i -band limit of 21 . h . Not unexpectedly, we also find thatcorrecting this SIE+ γ ext model with the average constantexternal convergence also leads to a ∼ H using lensed time delays can not be reliedupon. Time delay studies opt for lens systems which are ap-parently free of strong perturbers in an attempt to excludeline-of-sight effects. Nevertheless, our simulations have mim-icked this selection by removing any halo from all of ourlightcones that produces a deflection with an Einstein ra-dius larger than 0.3”. Our work reveals that more sophis-ticated lens modelling methods, likely including multiplelens planes, are required to reliably measure H from thehundreds of well-measured time-delay systems anticipatedin forthcoming large strong lens samples. MNRAS , 1–9 (20**) with L.O.S. halos ACKNOWLEDGEMENTS
The authors thank Sherry Suyu and Thomas Collett forinspiring discussion and suggestions. We are grateful toCharles R. Keeton and Masamune Oguri for taking theirtime to answer our questions. We are especially thank-ful to Simon Birrer for having made it possible to use lenstronomy for this project. NL and CB acknowledgessupport by the UK Science and Technology Facilities Coun-cil (STFC). SD is supported by the UK’s STFC ErnestRutherford Fellowship scheme. This research made use ofCosmoDC2 and GCR-Catalogs-Reader REFERENCES
Addison G. E., Watts D. J., Bennett C. L., Halpern M., HinshawG., Weiland J. L., 2018, ApJ, 853, 119Bennett C. L., et al., 2013, ApJS, 208, 20Birrer S., Amara A., 2018, ] 10.1016/j.dark.2018.11.002Birrer S., Amara A., Refregier A., 2015a, Astrophys. J., 813, 102Birrer S., Amara A., Refregier A., 2015b, ApJ, 813, 102Birrer S., Amara A., Refregier A., 2016, J. Cosmology Astropart.Phys., 2016, 020Birrer S., Welschen C., Amara A., Refregier A., 2017, J. Cosmol-ogy Astropart. Phys., 4, 049Birrer S., et al., 2019, MNRAS, 484, 4726Bolton A. S., et al., 2012, ApJ, 757, 82Bonvin V., et al., 2017, MNRAS, 465, 4914Chen G. C. F., et al., 2019, MNRAS, 490, 1743Collett T. E., 2015, ApJ, 811, 20Collett T. E., Cunnington S. D., 2016, MNRAS, 462, 3255Collett T. E., et al., 2013, MNRAS, 432, 679DES Collaboration et al., 2020, arXiv e-prints, p.arXiv:2002.11124Dhawan S., Jha S. W., Leibundgut B., 2018, A&A, 609, A72Dobler G., Fassnacht C. D., Treu T., Marshall P., Liao K., HojjatiA., Linder E., Rumbaugh N., 2015, ApJ, 799, 168Eberhart R., Kennedy J., 1995, in Proceedings of the IEEE in-ternational conference on neural networks. pp 1942–1948Faber S. M., et al., 2007, ApJ, 665, 265Freedman W. L., et al., 2019, ApJ, 882, 34Greene Z. S., et al., 2013, ApJ, 768, 39Habib S., et al., 2016, New Astron., 42, 49Heitmann K., et al., 2019, ApJS, 245, 16Keeton C. R., 2001, ArXiv Astrophysics e-prints,Komatsu E., et al., 2011, ApJS, 192, 18Koopmans L. V. E., Treu T., Bolton A. S., Burles S., MoustakasL. A., 2006, ApJ, 649, 599Kormann R., Schneider P., Bartelmann M., 1994, A&A, 284, 285Korytov D., et al., 2019, ApJS, 245, 26Li N., Gladders M. D., Rangel E. M., Florian M. K., Bleem L. E.,Heitmann K., Habib S., Fasel P., 2016, ApJ, 828, 54 https://portal.nersc.gov/project/lsst/cosmoDC2 https://github.com/LSSTDESC/gcr-catalogs Liao K., et al., 2015, Astrophys. J., 800, 11Macaulay E., et al., 2019, MNRAS, 486, 2184McCully C., Keeton C. R., Wong K. C., Zabludoff A. I., 2014,MNRAS, 443, 3631McCully C., Keeton C. R., Wong K. C., Zabludoff A. I., 2017,ApJ, 836, 141Momcheva I. G., Williams K. A., Cool R. J., Keeton C. R.,Zabludoff A. I., 2015, ApJS, 219, 29More A., et al., 2016, MNRAS, 455, 1191Mu˜noz J. B., Kamionkowski M., 2017, Phys. Rev. D, 96, 103537Narayan R., Bartelmann M., 1996, ArXiv Astrophysics e-prints,Oguri M., Marshall P. J., 2010, MNRAS, 405, 2579Parker L. C., Hoekstra H., Hudson M. J., van Waerbeke L., Mel-lier Y., 2007, ApJ, 669, 21Pietrzy´nski G., et al., 2019, Nature, 567, 200Planck Collaboration et al., 2018, arXiv e-prints,Riess A. G., 2019, Nature Reviews Physics, 2, 10Riess A. G., Casertano S., Yuan W., Macri L. M., Scolnic D.,2019, ApJ, 876, 85Rusu C. E., et al., 2017, MNRAS, 467, 4220Rusu C. E., et al., 2019, arXiv e-prints, p. arXiv:1905.09338Schneider P., Sluse D., 2013, A&A, 559, A37Schneider P., Ehlers J., Falco E. E., 1992, Gravitational Lenses,doi:10.1007/978-3-662-03758-4.Sereno M., Paraficz D., 2014, MNRAS, 437, 600Shajib A. J., et al., 2019, arXiv e-prints, p. arXiv:1910.06306Shu Y., et al., 2016, ApJ, 833, 264Sluse D., et al., 2017, MNRAS, 470, 4838Suyu S. H., Marshall P. J., Auger M. W., Hilbert S., BlandfordR. D., Koopmans L. V. E., Fassnacht C. D., Treu T., 2010,ApJ, 711, 201Suyu S. H., et al., 2013, ApJ, 766, 70Suyu S. H., et al., 2017, MNRAS, 468, 2590Tagore A. S., Barnes D. J., Jackson N., Kay S. T., Schaller M.,Schaye J., Theuns T., 2018, MNRAS, 474, 3403Tie S. S., Kochanek C. S., 2018, MNRAS, 473, 80Tihhonova O., et al., 2018, MNRAS, 477, 5657Verde L., Treu T., Riess A. G., 2019, Nature Astronomy, 3, 891Wertz O., Orthen B., Schneider P., 2018, A&A, 617, A140Wong K. C., Keeton C. R., Williams K. A., Momcheva I. G.,Zabludoff A. I., 2011, ApJ, 726, 84Wong K. C., et al., 2017, MNRAS, 465, 4895Wong K. C., et al., 2018, ApJ, 867, 107Wong K. C., et al., 2019, arXiv e-prints, p. arXiv:1907.04869Xu D., Sluse D., Schneider P., Springel V., Vogelsberger M., Nel-son D., Hernquist L., 2016, MNRAS, 456, 739Yuan W., Riess A. G., Macri L. M., Casertano S., Scolnic D. M.,2019, ApJ, 886, 61MNRAS000