© ESO 2020

1. Introduction

Time delay cosmography with multiply imaged quasars delivers measurements of the Hubble constant H₀ to percent-level precision (Chen et al. 2019; Birrer et al. 2019; Wong et al. 2020; Shajib et al. 2019), and substructure lensing studies use strong gravitational lenses to measure the halo mass function below 10⁹ M_⊙, thereby constraining the nature of dark matter (Dalal & Kochanek 2002; Nierenberg et al. 2014; Hezaveh et al. 2016a; Birrer et al. 2017; Vegetti et al. 2018; Hsueh et al. 2020; Gilman et al. 2020a,b). To date, these two lines of research are pursued independently. Dark substructure is omitted from the lens models used in time delay cosmography, while substructure lensing analyses do not explicitly marginalize over the uncertainties of cosmological parameters. In this paper, we move to bridge the gap between these fields by quantifying how dark substructure affects the precision and accuracy of measurements of H₀ obtained from strong lenses.

The impetus to carry out this analysis now stems from a recent challenge to the standard cosmological model including cold dark matter and a cosmological constant Λ (ΛCDM). Strong-lensing measurements of H₀ in the local Universe presented by the H0LiCOW collaboration (Suyu et al. 2017), and from Type Ia supernova carried out by the SH0ES collaboration (Riess et al. 2019), differ from the value of H₀ derived from the cosmic microwave background (CMB; Planck Collaboration I 2020). The tension between the CMB data and six strongly lensed quasars alone is 3.1σ, and exceeds 5σ when lensing is combined with supernova data (Wong et al. 2020; Riess 2019). A resolution to this tension may require new physics. Before reaching this conclusion it is necessary to quantify potential sources of bias and measurement uncertainties relevant for time delay cosmography.

Oguri (2007) pointed out that dark matter subhalos surrounding the main deflector in galaxy-scale strong lenses alter the arrival times between lensed images. The arrival times provide a direct, one-step method to measure H₀, so it is conceivable that perturbations from substructure to these data could bias, or contribute to the error budget, of cosmographic analyses. Keeton & Moustakas (2009; see also Congdon et al. 2010; Cyr-Racine et al. 2016) carried out the first comprehensive analysis of time delay perturbations from substructure, and showed that populations of dark subhalos impart arrival time perturbations on the scale of a few hours. This level of perturbation lies safely below current measurement precision, which measures image arrival times down to a few-days precision (e.g., Bonvin et al. 2018; Courbin et al. 2018; Millon et al. 2020a). On these grounds, substructure is assumed to be negligible and it is omitted from the strong lens models used for cosmography.

However, as the precision of time delay cosmography improves it is timely to revisit this assumption and quantify the potential impact of substructure both in terms of accuracy and precision. First, Keeton & Moustakas (2009) did not account for field halos along the line of sight to strong lenses, which can significantly outnumber subhalos. Because the number of line-of-sight halos depends on the lens and source redshifts, the amplitude of time delay perturbations should evolve with redshift, as should any resulting uncertainties and/or bias in the H₀ value inferred from the time delays. Second, it is possible that perturbations from dark matter halos on lensing observables propagate to H₀ inferences through indirect channels, not directly through the time delays. The collective lensing effects of a population of halos could be absorbed by the global mass model of the lens system, or the “macromodel”, in which case substructure would impart stochastic perturbations to the macromodel.

In this work, we quantify these effects by applying the same modeling assumptions used in cosmographic analyses to mock lenses that include full populations of subhalos and line-of-sight halos. Using the image positions and the surface brightness of lensed arcs to constrain the macromodel, we infer H₀ from the mock lenses with lens models that do not include substructure, as is currently done in the literature. Comparing the inferred H₀ with the “true” value of H₀ used to create the mocks, we quantify to what degree dark substructure would bias or increase the uncertainties reported in cosmographic inferences were it included in the lens model. The mock lenses are based on six real quadruply imaged quasar systems modeled by the H0LiCOW, STRIDES, and SHARP teams (from now on, we refer to this collaboration as TDCOSMO) that have delivered measurements of H₀, allowing us to estimate the effects of substructure for each individual system. By construction, our simulations isolate the role of substructure from other sources of uncertainty, such as the mass sheet transformation (Schneider & Sluse 2013) and systematic uncertainties associated with the lens macromodel (e.g., Kochanek 2020).

A necessary component of our analysis is the creation of hundreds of mock lenses with full line-of-sight and subhalo populations. The dark matter halos in the mock systems simultaneously perturb the arrival times and the flux ratios. In addition, halos occasionally produce gravitational imaging residuals in extended images that lead to the discovery of individual dark matter halos in strong lens systems, a technique known as gravitational imaging (Koopmans 2005; Vegetti et al. 2010, 2018; Hezaveh et al. 2016a). Based on our simulations, we comment on results relevant for gravitational imaging, and in particular recent suggestions to measure the convergence power spectrum of dark matter halos (e.g., Hezaveh et al. 2016a; Cyr-Racine et al. 2016). As the focus of this work is time delay cosmography, we relegate discussion of these topics to Appendices B and D.

This paper is organized as follows: Section 2 describes our approach for estimating H₀ uncertainties from simulations that include substructure in mock lenses. Section 3 describes the specific modeling choices implemented in the simulations. Section 4 describes the results of the analysis, which we summarize concisely in Sect. 5. We make concluding remarks in Sect. 6.

We perform gravitational lensing computations with lenstronomy¹ (Birrer & Amara 2018). Computations that involve the halo mass function and the matter power spectrum are executed with colossus (Diemer 2018). We assume a flat ΛCDM cosmology using the parameters from WMAP9 (Hinshaw et al. 2013) (Ω_m = 0.28, σ₈ = 0.82), and set h = 0.733 (see Sect. 2). The conclusions we present are independent of the choice of cosmological parameters.

2. Cosmographic inference on mock lenses that include substructure

This section describes our method to estimate the effect of dark substructure on cosmographic inferences. A brief outline is as follows: First, we create mock datasets with substructure included in the lens models. Next, we analyze the data as would be done in a modern cosmographic inference, ignoring substructure in the model. Finally, we investigate how the omission of dark substructure affects cosmographic inferences through comparison with a reference dataset that contains no substructure.

We begin by introducing the time delay distance D_Δt, and describe the observables from strong lens systems used to constrain it, in Sect. 2.1. In Sect. 2.2, we discuss the task of inferring H₀ in the presence of full populations of dark matter halos in a Bayesian context. Section 2.3 describes how we create mock datasets with substructure included the lens model. Section 2.4 details our procedure for using these mock datasets to quantify any bias or additional uncertainties associated with the presence of substructure in strong lenses.

2.1. The time delay distance and strong lensing observables

Strong lensing constrains the Hubble constant through a measurement of an absolute distance scale D_Δt, referred to as the time delay distance $D_{\mathrm{\Delta }t}=(1+z_{\mathrm{d}})\frac{D_{\mathrm{d}}{D}_{\mathrm{s}}}{D_{\mathrm{ds}}}$. The factors D_d, D_s, D_d, s are angular diameter distances to the deflector, the source, and from the deflector to the source², and z_d is the deflector redshift. The combination of cosmological distances entering D_Δt depends on cosmology, with a particularly strong dependence on H₀, with $D_{\mathrm{\Delta }t}\propto \frac{1}{H_0}$.

The data vector d = (θ_img,t,s) that constrains quantities such as D_Δt includes the four lensed quasar image positions θ_img = (θ_A,θ_B,θ_C,θ_D), the time delays between the images t = (t₁₂,t₁₃,t₁₄), and the observed surface brightness s(θ) at angular position θ in the image plane. All matter in the lens system, dark or baryonic, affects the time delays, image positions, and the arc surface brightness. Two equations encode the dependence of these observables on the deflection angles α, and on angular diameter distances and cosmology. First, the ray tracing equation describes the propagation of light rays through lines of sight populated by dark matter halos (Blandford & Narayan 1986)

$$\begin{array}{c}\hfill {\mathit{\theta }}_K\left(\mathit{\theta }\right)=\mathit{\theta }-\frac{1}{D_{\mathrm{s}}}{\displaystyle \sum _{k=1}^{K-1}{D}_{k,\mathrm{s}}{\mathit{\alpha }}_k\left(D_k{\mathit{\theta }}_k\right).}\end{array}$$(1)

This is a recursive relation for angular coordinates on the kth lens plane, given a coordinate on the sky θ, the deflection angle α, and the angular diameter distance from the kth lens plane to the source plane, D_k, s. Equation (1) gives the surface brightness in the image plane s(θ) by ray-tracing through the lens model to the source plane, and the multiple images of the background quasar appear at the coordinates θ that each map to the same position on the source plane.

The time delay between images located at coordinates θ_A and θ_B is

$$\begin{array}{c}\hfill t_{\mathit{AB}}=\frac{D_{\mathrm{\Delta }t}^{\mathrm{eff}}}{c}({\varphi }^{\mathrm{eff}}\left({\mathit{\theta }}_A\right)-{\varphi }^{\mathrm{eff}}\left({\mathit{\theta }}_B\right)).\end{array}$$(2)

where the effective Fermat potential is given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\varphi }^{\mathrm{eff}}\left({\mathit{\theta }}_X\right)\equiv & {\displaystyle \sum _{k=1}^{K-1}\frac{1+z_k}{1+z_{\mathrm{d}}}\frac{D_k{D}_{k+1}{D}_{\mathrm{ds}}}{D_{\mathrm{d}}{D}_{\mathrm{s}}{D}_{k,k+1}}}\hfill \\ \hfill & \times [\frac{{({\mathit{\theta }}_k^X-{\mathit{\theta }}_{k+1}^X)}^2}{2}-{\xi }_{k,k+1}\mathrm{\Psi }\left({\mathit{\theta }}_k^X\right)],\hfill \end{array}\end{array}$$(3)

where Ψ(θ) is the projected gravitational potential at position θ, ${\xi }_{i,j}\equiv \frac{D_{i,j}{D}_{\mathrm{s}}}{D_i{D}_{i,\mathrm{s}}}$, and where we use shorthand notation ${\mathit{\theta }}_K\left({\mathit{\theta }}_X\right)\equiv {\mathit{\theta }}_K^X$ for Eq. (1) evaluated at a specific angular position θ_X. The subscript d refers to the main deflector lens plane.

Assuming a cosmology, there is a unique distance scale $D_{\mathrm{\Delta }t}^{\text{eff}}$, given a value of H₀. Following the standard approach (e.g., Wong et al. 2020; Shajib et al. 2020), we can define an effective time delay distance $D_{\mathrm{\Delta }t}^{\text{eff}}$ that normalizes the measured time delays and the model-predicted Fermat potential, and compute H₀ from this effective distance scale.

A typical cosmographic analysis includes nuisance parameters n, such as the lens mass profile, the lens light profile, the source light profile, and other factors that can affect the inference such as the microlensing time delay (Tie & Kochanek 2018), and the external convergence associated with large-scale structure near the lens and along the line of sight (Schneider & Sluse 2013; Birrer et al. 2016). Obtaining a D_Δt posterior from a strong lens system requires marginalizing the likelihood function over these nuisance parameters

$$\begin{array}{c}\hfill p\left(\mathit{d}|D_{\mathrm{\Delta }t}\right)={\displaystyle \int }\mathcal{L}(\mathit{d}|D_{\mathrm{\Delta }t},\mathit{n})p\left(\mathit{n}\right)\mathrm{d}\mathit{n}.\end{array}$$(4)

The evaluation of Eq. (4) involves simultaneously reconstructing the lensed background source and the lensed image.

2.2. Directly computing H₀ posteriors with substructure

To incorporate substructure, we describe dark field halos and subhalos with population-level statistical properties specified by a set of hyper-parameters q_cdm. These hyper-parameters describe, for example, the form of the field halo and subhalo mass functions, the average concentration-mass relation of halos, the probability distribution for their spatial distribution, etc., assuming CDM. Typically, cosmographic analyses ignore q_cdm altogether when measuring D_Δt, assuming p(D_Δt|d,q_cdm) ≈ p(D_Δt|d). In other words, the assumption in the literature is that inferences on H₀ do not depend on the nature of small-scale dark matter structure in the lens system, or at least that other factors in the analyses, such as the nuisance parameters n, dominate the D_Δt error budget such that dark substructure can be safely ignored.

The direct route toward testing this assumption involves a daunting calculation of a likelihood function. We can account for dark substructure in Eq. (4) by generating realizations of dark matter halos specified by parameters m_sub from the model specified by hyper-parameters q_cdm. The likelihood that accounts for q_cdm includes an additional integral over all possible configurations m_sub

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{L}(\mathit{d}|D_{\mathrm{\Delta }t},{\mathit{q}}_{\mathrm{cdm}})=& {\displaystyle \int }p(\mathit{d}|D_{\mathrm{\Delta }t},{\mathit{m}}_{\mathrm{sub}},\mathit{n})\hfill \\ \hfill & \times p\left({\mathit{m}}_{\mathrm{sub}}|{\mathit{q}}_{\mathrm{cdm}}\right)p\left(\mathit{n}\right)\mathrm{d}{\mathit{m}}_{\mathrm{sub}}\mathrm{d}\mathit{n}.\hfill \end{array}\end{array}$$(5)

Incorporating substructure by directly evaluating the integral over m_sub is not feasible because Eq. (5) is computationally intractable. The problem stems mainly from integrating p(d|D_Δt,m_sub,n) over n, because m_sub includes potentially tens of thousands of halos whose deflection angles need to be continuously re-evaluated while performing the integral over n.

2.3. Forward modeling substructure perturbations on lensing data

To circumvent the direct evaluation of the intractable integral in Eq. (5), we instead create mock datasets computed with substructure in the lens model using a a known true value of H₀. We then apply same assumptions used in the literature (ignore substructure) to analyze these systems with the goal of recovering the true H₀ value. The degree to which we fail to recover the true H₀, and the degree to which the uncertainties increase relative to a control sample with no substructure in the lens system, allowing us to to quantify any bias or increased uncertainties associated with substructure.

The mock lenses we generate are analogs of the lens systems RXJ1131−1231, PG1115+080, HE0435−1223, B1608+656, WFI2033−4723, and DESJ0408+5354, each of which has been used to measure H₀ (see the discussion in Sect. 3.5 for details). Through the procedure outlined in the remainder of this section, we determine on a case-by-case basis to what degree omitting substructure from the lens model results in under-estimated uncertainties and/or bias in the inferred value of H₀.

To begin, we describe two sets of mock data created for each lens system: the control data and the substructure-perturbed data:

– The control data: The control data includes no substructure, and the same parameterization for the mass model is used to create the data as is used to fit the resulting mock lens. The uncertainty in H₀ inferred from this dataset is determined entirely by the measurement precision of the time delays, and the degree to which the lensed images and extended arc constrain the mass profile. We tune the signal to noise in the mock data until the H₀ inference we obtain from this dataset has the same error budget from time delays and lens modeling as those quoted in Table 3 of Wong et al. (2020), or the uncertainties reported by Shajib et al. (2020) in the case of DESJ0408.

– The substructure-perturbed data: We first generate a realization of dark matter field halos and subhalos. Next, starting from the lens model used to create the control dataset, we solve for a new set of macromodel parameters that map the observed image positions to the same coordinate on the source plane through Eq. (1). The new macromodel is initialized from the macromodel used to create the control data, so the new macromodel is nearby in parameter space to the one used to create the control dataset. We place the same extended background source used to create the control data at the new source position obtained from adjusting the macromodel, and proceed to compute the time delays and the surface brightness of the extended images. The mock lens obtained through this process is indistinguishable by eye from the control data, but the time delays and lensed arc include the subtle lensing signatures of dark matter substructure.

Both the control and substructure-perturbed models are created using the same value of H₀ = 73.3 km s⁻¹ Mpc, the value quoted by Wong et al. (2020), and we apply the same Gaussian point spread function (PSF) with a full width at half maximum (FWHM) of 50 milli-arcseconds to the simulated lenses, comparable to the imaging resolution from the Hubble Space Telescope. We stress that our results are independent of the specific choice of H₀.

2.4. Estimating H₀ uncertainties from un-modeled dark substructure

We fit the mock datasets using smoothly parameterized mass profiles for the main deflector mass profile while ignoring substructure in the lens system. The model applied to the control data is the same as the model used to create the data, so we should fit the mock data down to the noise; this is a perfect model scenario. We then apply the same modeling assumptions to the substructure-perturbed datasets, fitting them with the same smoothly-parameterized mass profile. To isolate the role of substructure, we assume perfectly measured image positions and do not add astrometric uncertainties to the lensed point sources. The flux ratios are not used to constrain the lens model, although we can compute them jointly with the time delays (see Appendix D). We perform the lens modeling on 200 mock datasets, each of which corresponds to one realization of dark matter halos. This gives 200 corresponding posterior distributions of H₀, which we combine by adding them together.

In our procedure, mimicking the actual measurements, we include additional model complexity in the mock datasets, rather than in the model applied to the data. The requirement that must be met in order to safely make this approximation is that each of the substructure-perturbed datasets wield the same constraining power over H₀. We take two measures to ensure this is the case: first, because the mock “measured” time delays change slightly from realization to realization in the substructure-perturbed datasets, a fixed uncertainty in units of days translates to different percent precision in H₀. To keep the same constraining power of the time delays the same for each individual realization, we operate with relative uncertainties on the time delays, using the same uncertainties as those measured from the real lens systems studied in the literature.

Second, by re-weighting posterior samples we ensure that we operate at fixed uncertainty in the logarithmic slope of the main deflector mass profile γ, and the Einstein radius R_Ein. Operating at fixed uncertainty in γ and R_Ein serves two purposes: First, power-law mass models fit to the lens system HE0435−1223 (to use a specific example) have logarithmic slopes γ = 1.9 ± 0.05 (Chen et al. 2019). Substructure realizations that lead to global mass profiles with profile slopes substantially different than this are therefore not analogs of the lens in question. Quantitatively, the realizations that lead to this mismatch between the model and the data do not contribute to the integral in Eq. (5), because the likelihood associated with these realizations would be extremely low. This argument applies to both γ and R_Ein. The second reason for operating at fixed γ and R_Ein is that both of these parameters are, in principle, constrained by the stellar kinematics of the main deflector. Enforcing equality between the γ and R_Ein posteriors ensures that the modeled-predicted velocity dispersion matches the velocity dispersion that would have been “measured” from the mock datasets. By incorporating information that could, in principle, be constrained by kinematic data, we isolate the role of substructure in the lens system.

After this series of steps, we are ready to make comparisons between the H₀ inference from the control data and the inferences from the substructure-perturbed data. If substructure is irrelevant for time delay cosmography, then we expect the H₀ posteriors obtained from the control data and the substructure-perturbed dataset to be identical. The degree to which the H₀ posteriors differ allows us to quantify how the reported values of H₀ from the TDCOSMO analyses would have been affected, had substructure been included in the lens model. Compared to the control data, the H₀ inferences from the substructure-perturbed data will be inflated by some amount σ_sub, independent of the other sources of uncertainty σ_reported. The total uncertainty is therefore

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{total}}^2={\sigma }_{\mathrm{reported}}^2+{\sigma }_{\mathrm{sub}}^2.\end{array}$$(6)

A bias in H₀ will manifest as different medians between the control and substructure-perturbed posteriors.

3. Simulation setup

The procedure outlined in the previous section requires a prescription for rendering populations of dark matter halos, and models for the lens mass distribution, lens light profile, and source light profile. We describe the parameterizations and modeling choices for the subhalo and line-of-sight halo populations in Sect. 3.1. Section 3.2 describes how subhalos and line-of-sight halos are generated symmetrically around the light rays throughout the lensing volume, and the addition of negative convergence profiles to correct for the additional mass added in substructure. Sections 3.3 and 3.4 describe the models we implement for the main deflector mass profile and the background source light profile, respectively. We reiterate that the lens models we use omit complexity in the form of external convergence or more complicated mass profiles in order to isolate the perturbative effects of substructure.

3.1. The CDM subhalo and field halo populations

We use the same parameterizations for the CDM mass function as those implemented by (Gilman et al. 2020a,b); we briefly review the essential features in the following paragraphs. Both subhalos and field halos along the line of sight, are modeled as truncated Navarro-Frenk-White (NFW) profiles (Navarro et al. 1996; Baltz et al. 2009) defined with a density contrast of 200 with respect to the critical density at z = 0. We use the mass-concentration relation presented by Diemer & Joyce (2019). The truncation radius for subhalos depends on the mass of the subhalo, and the position inside the host halo, while field halos are truncated at r₅₀. The choice of r₅₀ for field halos is motivated by the splashback radius, a physical halo boundary that is typically larger than, but the same order of magnitude, as the virial radius (Diemer 2020). The truncation of field halos at radii > r₂₀₀ keeps the total mass in line-of-sight halos finite with a negligible impact on their lensing properties.

We draw field halo masses from a Sheth-Tormen mass function (Sheth et al. 2001) with the addition of the two-halo term ξ_2halo(M_halo,z)

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{N}_{\mathrm{los}}}{\mathrm{d}m\mathrm{d}V}=(1+{\xi }_{2\mathrm{halo}}(M_{\mathrm{host}},z))\frac{{\mathrm{d}}^{2}N}{\mathrm{d}m\mathrm{d}V}{|}_{\mathrm{Sheth}\text{-}\mathrm{Tormen}}.\end{array}$$(7)

The two-halo term boosts the amount of structure near the host halo, which has a mass M_host. We render subhalos from a power-law mass function defined in projection

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{N}_{\mathrm{sub}}}{\mathrm{d}m\mathrm{d}A}=\frac{{\mathrm{\Sigma }}_{\mathrm{sub}}}{m_0}{\left(\frac{m}{m_0}\right)}^{\alpha }\mathcal{F}(M_{\mathrm{host}},z)\end{array}$$(8)

with logarithmic slope α = −1.9 (Springel et al. 2008) and a pivot mass m₀ = 10⁸ M_⊙. The function ℱ(M_host,z) encodes the evolution of the projected number density of subhalos with host halo mass and redshift (see Gilman et al. 2020a), and is calibrated using the semi-analytic modeling code galacticus (Benson 2012). It is given by ${log}_{10}\left(\mathcal{F}\right)=k_1{log}_{10}\left(\frac{M_{\mathrm{host}}}{{10}^{13}{M}_{\odot }}\right)+k_2{log}_{10}(1+z)$, with k₁ = 0.88 and k₂ = 1.7. We assume a host halo mass of 10^13.3 M_⊙ for each mock lens, based on the population mean halo mass for strong lensing galaxies inferred by Lagattuta et al. (2010). We use a normalization Σ_sub = 0.025 kpc⁻² based on the measurements from (Gilman et al. 2020a). These modeling choices are intended to be representative of the subhalo populations around massive elliptical galaxies.

We render both halos and subhalos in the mass range 10⁶ − 10⁹ M_⊙. We have carried out tests with a lower minimum mass threshold and verified that our results do not depend on this choice. We assume halos more massive than 10⁹ would host a visible galaxy, and would therefore be modeled explicitly in cosmographic analyses. The choice of 10⁹ solar masses is conservative in the sense that a larger mass would increase the uncertainties we attribute to un-detected substructure. The rendering volume is a double-cone geometry with an opening angle of 15 times the Einstein radius of each lens. We have tested larger opening angles, and verify that our results are unchanged. The total number of halos and subhalos included in lens models ranges from 50 000 − 200 000 depending on the lens and source redshifts and the Einstein radius.

3.2. Lens cone symmetry and convergence corrections

At each lens plane along the line of sight³, we generate halos uniformly distributed in two dimensions around a coordinate center that tracks the path of a lensed light ray backward from the observer toward the lens mass centroid. Behind the main deflector, the rendering area tracks the light as it is deflected by background galaxies, as is the case in HE0435, WFI2033, and DESJ0408. Without shifting the rendering area to track lensed light rays, substructure is not distributed symmetrically around the path of the light, and the asymmetric mass distribution biases the inferred value of H₀. The shifting of substructure along a lensed light path is particularly important for lens systems with large satellites and/or line-of-sight galaxies, as the source coordinate in these systems can be offset from the lens center by as much as 1 arcsec, comparable to a typical lens Einstein radius.

We add negative convergence at each lens plane to subtract the mean mass rendered in (sub)halos, ensuring that the average dark matter density in the rendering volume is equal to the mean dark matter density of the Universe. For line-of-sight halos, the correction is in the form of a uniform (negative) convergence sheet. The correction is more complicated for subhalos, as they are spatially distributed following a cored NFW profile. The cored spatial distribution, with a central core radius set to half the scale radius of the host, is intended to mimic the tidal disruption of subhalos that pass close to the central galaxy. The net convergence profile that needs to be subtracted for the subhalos is formally given by a convolution of the subhalo density profile with their spatial distribution function, but we are able to approximate this profile with a simple cored NFW profile normalized such that the residual convergence from substructure near the Einstein radius is zero on average.

Omitting the mass sheet corrections trivially biases H₀ by re-scaling each time delay through the addition of a net positive convergence in the form of dark matter halos. This effect is similar to the external convergence that is accounted for in cosmographic inferences, which affects the inferred time delay distance as

$$\begin{array}{c}\hfill D_{\mathrm{\Delta }t}=\frac{{D_{\mathrm{\Delta }t}}_{\mathrm{model}}}{1-{\kappa }_{\mathrm{ext}}}·\end{array}$$(9)

While the ensemble of substructure-perturbed datasets have a residual convergence of zero on average, on a realization-by-realization basis the net residual convergence from substructure fluctuates around the mean dark matter density in the Universe. This introduces a mass-sheet transformation from random populations of dark matter halos that is distinct from the external convergence κ_ext, which comes from large-scale structure in the lensing volume. The fluctuations of the substructure convergence cause fluctuations in the model-predicted time delays, and therefore contribute to the H₀ error budget.

3.3. The main deflector mass and light profiles

Since our goal is to isolate the effects of substructure, we use the same mass and light profiles to create the control and substructure-perturbed datasets as we use in the lens models fit to these data. The main deflector is modeled as a power-law ellipsoid with logarithmic mass profile slope γ, plus external shear γ_ext. To keep the analogs as similar to the real lenses as possible, we choose a logarithmic profile slope γ consistent with the mean value inferred in the lens modeling by TDCOSMO. While the degree to which the macromodel absorbs the effects of substructure depends in principle on how it is parameterized, we do not expect the specific choice of a power law ellipsoid to significantly affect our results. Substructure imparts lensing effects on much smaller scales that those that can be subsumed into most reasonable parameterizations of the macromodel⁴.

All of the lens systems we consider have either a nearby satellite galaxy in the main lens plane, a massive galaxy along the line of sight near the main deflector (in projection), or both. We model these objects as singular isothermal spheres, using Einstein radii and positions inferred by the lens modeling from the TDCOSMO teams. We add an elliptical Sérsic light profile for the main deflector light, and circular Sérsic light profiles for the satellite or nearby galaxies. The parameters describing the macromodels used to to create the mock lenses are listed in Table A.1.

3.4. Background source light profiles

The mock datasets are all created using an elliptical Sérsic model for the background source. The exact parameters describing these profiles are unique to each lens, and are chosen such that the lensed arc in each mock resembles the extended features in the real lens it is based on.

Modeling lenses for cosmography involves a simultaneous reconstruction of the background source light and the lensed image. Using an elliptical Sérsic profile to model the background source (the same model as is used to create the mock lens), we detect significant residuals in the reconstructed lensed images that result from substructure in the lens system. Following the approach taken in the literature, when we detect these residuals we add additional complexity to the background source. We add complexity to the background source model by adding shapelets (Refregier 2003; Birrer et al. 2015), with the amount small-scale structure in the source controlled by the parameter n_max, which we set to n_max = 8. In the case of DESJ0408+5354, we follow Shajib et al. (2020) and add additional set of shapelets with n_max = 3 to one of the additional lensed sources. Adding complexity to the background source and fitting the lens model to the mock data removes the residuals in the image plane. This suggests the substructure-induced imaging residuals depend explicitly on the parameterization of the background light source. We elaborate on this topic in Appendix B.

3.5. Discussion of individual lens systems

We construct realistic analogs of six strong lens systems analyzed by the TDCOSMO collaboration. We briefly comment on each system in following paragraphs, highlighting both cosmographic inferences and substructure lensing analyses performed with each system. The control lens models used for each system are summarized in Table A.1.

RXJ1131−1231. This quad was discovered by Sluse et al. (2003), and consists of source at z = 0.658 lensed by a galaxy at z = 0.295 (Sluse et al. 2003, 2007). The time delays are measured at ∼1.6% precision by Tewes et al. (2013), and exhibit a relatively long maximum delay of ∼90 days. This system has been modeled extensively with the aim of constraining cosmology (Suyu et al. 2013, 2014; Birrer et al. 2016; Chen et al. 2016, 2018). The lens mass model includes a satellite galaxy inside the Einstein radius. Perturbations from substructure in the lensed arc of this system were analyzed by Birrer et al. (2017), and were used to place a lower limit on the mass of a thermal relic dark matter particle of 2.2 keV.

PG1115+080. This quad was discovered by Weymann et al. (1980), and consists of a source at redshift z = 1.722 lensed by a galaxy at z_d = 0.310 (Tonry 1998). The time delays were measured at ∼6.4% precision by Bonvin et al. (2018), and the longest delay is ∼18 days. This system was recently modeled with both adaptive optics and HST imaging by (Chen et al. 2019). We explicitly model the nearby cluster of galaxies as a single SIS profile projected ∼10 arcsec away from the lens center. This system was first used to measure the abundance of dark subhalos by Dalal & Kochanek (2002), and more recently was used to constrain the free-streaming length of dark matter Hsueh et al. (2020), Gilman et al. (2020a), and the mass-concentration relation of CDM halos (Gilman et al. 2020b).

HE0435+1223. This quad was discovered by Wisotzki et al. (2002), with lens/source redshifts of z = 0.45 and z_s = 1.69, respectively (Morgan et al. 2005). The time delays were measured at ∼5.3% precision (Bonvin et al. 2017), with a maximum delay of ∼14 days. This system was analyzed for cosmography by Bonvin et al. (2017), Wong et al. (2017), Chen et al. (2019). We explicitly model the galaxy located behind the main deflector at z = 0.78 as an SIS profile, placing it at an angular position corrected by ∼0.65 arcsec due to foreground lensing effects from the main deflector⁵. Nierenberg et al. (2017) measured narrow-line flux ratios for this system (see also Nierenberg et al. 2020), and used them to look for substructure near the images. The narrow-line flux ratio measurements for this system were also used measure the halo mass function and mass-concentration relation of CDM halos (Gilman et al. 2020a,b).

B1608+656. This lens system was discovered by Myers et al. (1995), with time delay measurements and early cosmographic inferences carried out by Fassnacht et al. (2002) and Koopmans et al. (2003). The lens/source redshifts for this system are z_d = 0.45, and z_s = 1.69 (Fassnacht et al. 1996; Myers et al. 1995). The longest delay in this system is ∼76 days, measured at ∼2% precision. This system is unique in that main deflector consists of two merging galaxies, although (Suyu et al. 2009) show it is still well-described by a single power law profile. We follow Suyu et al. (2009) and model the larger object as a power law ellipsoid, and included smaller object as an SIS profile.

WFI2033−4723. This quad was discovered by Morgan et al. (2004). The lens and source lie at z = 0.66 and z = 1.66, respectively (Sluse et al. 2019). The time delays were measured by Bonvin et al. (2019) with a precision of ∼2%, and a maximum delay of ∼60 days. The main lensing galaxy has a small satellite located ∼2 arcsec away from the lens centroid. In addition, there is a galaxy located behind the main deflector at z = 0.745 situated ∼3.5 arcsec away from the lens center⁶. This system was modeled for cosmography by Rusu et al. (2020), and for substructure lensing inferences by Gilman et al. (2020a,b).

DESJ0408−5354. This lens systems stands out from the rest due to the presence of two strongly lensed sources, in addition to the quadruply imaged quasar at z = 2.355, discovered by Lin et al. (2017). The time delays for this system were measured by Courbin et al. (2018) at a precision of 1% with the longest delay ∼150 days. Shajib et al. (2020) modeled this system for cosmography, and exploited the additional information provided by the multiple lensed sources to measure H₀ at 3.8% precision (not including substructure). The main deflector is accompanied by a satellite galaxy that splits one of the lensed images, and several other galaxies along the line of sight. Shajib et al. (2020) modeled six of these line-of-sight galaxies, but in our mock lens models we only include the nearest two (G3 and G4, using the naming convention from Shajib et al. 2020). This simplification does not affect our conclusions, since we are matching the precision of the original inference by construction.

We do not explicitly model doubly imaged quasar SDSS 1206+4332, which was analyzed for cosmography by Birrer et al. (2019), because the analysis pipeline we develop is designed for quadruple-image lenses. However, based on the results obtained for the six systems we do model, we can estimate the fractional uncertainty in H₀ from substructure in SDSS 1206+4332 (see Sect. 4.2).

4. Results

In this section we apply the procedure outlined in Sect. 2 to the six strong lens systems described in the previous section. In Sect. 4.1, we show examples of the mock and reconstructed lensed images for each system, as well as maps of the residual time delay surface and convergence. Section 4.2 presents the main results of this analysis, in which we quantify how un-modeled dark substructure contributes to uncertainties and/or bias in H₀.

4.1. Reconstructed lenses: Imaging residuals, time delay anomalies, and convergence

Figures 1 and 2 show examples of the mock observed and reconstructed lenses for each of the six systems we analyze. The panels on the far left show the mock observed lens system, which includes substructure. The reconstructed lens shown in the second from left panels includes no substructure. The central panels show the imaging residuals between the mock observed lens and the reconstruction, and displays the reduced χ² obtained from the fit to the imaging data. We discard realizations with reduced χ² <  −1.1, as a structure massive enough to produced large imaging residuals would likely be detected and explicitly modeled. In Appendix B we comment on the implications of these imaging residuals for gravitational imaging of dark substructure, and in particular on the implications for recent proposals to measure the convergence power spectrum of halos through gravitational imaging residuals.

Fig. 1. Mock observed lens that includes substructure in the lens model (far left), the reconstructed lensed image (second from left), the normalized imaging residuals (center), the residual time delay surface (second from right), and the residual convergence (far right). The residual time delay surface is the true time delay surface from the mock lens minus the model-predicted time delay surface from the lens model fit to the data. Similarly, the residual convergence is the full multi plane convergence (Eq. (10)) minus the convergence from the model fit to the mock data. The lensed quasar image positions are marked as points (crosses) in the center (right) panels. Each row corresponds to one substructure-perturbed dataset for the lens systems based on RXJ1131−1231, PG1115+080, and HE0435−1223.

Fig. 2. The figure is the same as Fig. 1, but shows examples of the mock observed and reconstructed lens (left), the imaging residuals (center), the residual time delay surface (second from right), and the residual convergence (far right) maps for the lens systems B1608+656, WFI2033−4723, and DESJ0408−5354.

The second from right panels show the residual time delay surface, which is the difference between the true time delay surface of the control dataset, and the model-predicted time delay surface from the reconstructed lens system. Dark matter halos imprint small perturbations on the time delay surface, which, if they were somehow detected, could be interpreted as “time delay anomalies”, similar to the “flux ratio anomalies”, attributed to dark substructure.

The far right panels of Figs. 1 and 2 show the residual convergence⁷ between the full lens system that includes substructure, and the smooth model fit to the mock data. We use the definition of an effective convergence in terms of the deflection angles α (Gilman et al. 2019)

$$\begin{array}{c}\hfill {\kappa }_{\mathrm{eff}}=\frac{1}{2}(\frac{\partial \mathit{\alpha }}{\partial \mathit{x}}+\frac{\partial \mathit{\alpha }}{\partial \mathit{y}}),\end{array}$$(10)

which generalizes lensing convergence to systems with line-of-sight structure⁸. The residual convergence shown in the right-most panels of Figs. 1 and 2 is simply κ_eff − κ_{reconstructed}, where κ_{reconstructed} is the convergence from the smooth model fit to the mock data and κ_eff is computed from the mock lens system that includes substructure. Structures in the convergence maps that are tangentially sheared around the Einstein radius come from halos located between the main deflector and the source that produce highly nonlinear lensing effects.

4.2. H₀ inferences with dark substructure

Figures 3–8 show the joint posterior distribution for the inferred value of H₀ and the residual time delays Δt, defined as the true time delay from the mock data minus the model-predicted time delay. The posterior from the control data is shown in black, and the posterior obtained from the substructure-perturbed datasets are shown in blue. The joint posterior distributions of H₀ and the lens macromodel parameters are shown in Appendix C.

Fig. 3. The joint posterior distribution of the residual time delays Δt, defined as the model-predicted time delay minus the true time delay, and the Hubble constant, computed for the lens system RXJ1131−1231. The contours show the 68% and 95% confidence intervals for the control posterior that includes no substructure (black), and the combined posteriors from the inferences on lens models that include substructure (blue). Dashed lines in the histograms show 95% confidence intervals. The median and 68% confidence intervals of the H0 inference are printed above the figure.

Fig. 4. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system PG1115+080.

Fig. 5. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system HE0435−1223.

Fig. 6. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system B1608+656.

Fig. 7. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system WFI2033−4723.

Fig. 8. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system DESJ0408−5354.

The median inferred value of H₀ can differ from the true value of H₀ = 73.3 km s⁻¹ Mpc⁻¹ if the H₀ posterior is asymmetric, as is the case for PG1115+080. This does not suggest the presence of a bias in the inference, because a bias would manifest as different medians between the control and substructure-perturbed posteriors. We plot the difference between the medians of the two posteriors as a function of the square root of the lensing volume divided by the longest time delay (the significance of these units will be discussed in the following paragraph) in Fig. 9. The median values of the two posteriors are statistically consistent to within the error bars, which are given the uncertainty of the control dataset divided by $\sqrt{N}$, where N = 200 is the total number of realizations. We conclude that the omission of substructure from strong lens models used for cosmography does not introduce a bias greater than 0.3%.

Fig. 9. The y-axis shows the bias in H0, defined as the median of the H0 posterior obtained from the substructure-perturbed datasets minus the median of the posterior from the control datasets. The mock lenses are sorted on the x-axis by the square root of the lensing volume divided by the longest time delay, a dimensional quantity that determines the additional uncertainty caused by substructure in the lens system. The error bars show the statistical uncertainty on the median, which is given by $\frac{{\sigma }_{\mathrm{control}}}{\sqrt{N}}$, where σcontrol is the uncertainty in H0 from the control dataset and N = 200 is the number of realizations. We find no evidence for bias in the inferred value of H0 stemming from substructure that is not modeled in cosmographic inferences.

While we find no evidence for a bias in H₀, substructure alters the lens model and changes the Fermat potential at the image positions, which leads to an increased scatter in the time delays that factors directly into the H₀ error budget. The additional uncertainty caused by substructure evolves with redshift and geometry in a predictable way. From Eq. (2), the fractional uncertainty in H₀ scales as

$$\begin{array}{c}\hfill \frac{\delta H_0}{H_0}\propto \frac{\delta t}{t_{\max }}\end{array}$$(11)

where t_max is the longest time delay, which is usually the most precisely measured and therefore the most constraining observable, and δt is a perturbation to the time delays from dark matter halos. Using analytic arguments, Keeton & Moustakas (2009) and Cyr-Racine et al. (2016) showed that δt is proportional to $\sqrt{f_{\mathrm{sub}}}$, where f_sub is proportional to the total number of subhalos. We can generalize these results to realizations with line-of-sight halos by replacing f_sub by N_halo, where N_halo is just the total number of halos. The total number of halos is proportional to the lensing volume V, so the scaling in the fractional uncertainty in H₀ is

$$\begin{array}{c}\hfill \frac{\delta H_0}{H_0}\propto \frac{\sqrt{V}}{t_{\max }}·\end{array}$$(12)

The volume V is defined in terms of the comoving distances to the lens d_lens, to the source d_src, and the Einstein radius (converted to radians) R_Ein

$$\begin{array}{c}\hfill V=\pi {\left(d_{\mathrm{lens}}{R}_{\mathrm{Ein}}\right)}^2{d}_{\mathrm{src}}.\end{array}$$(13)

Figure 10 shows the trend predicted by Eq. (12) alongside the fractional uncertainties obtained for each lens. The dashed line fit in Fig. 10 is given by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{\delta H_0}{H_0}=0.045{\left(\frac{V}{1{\mathrm{Mpc}}^3}\right)}^{0.5}\left(\frac{10\mathrm{days}}{t_{\max }}\right)·\end{array}\end{array}$$(14)

Fig. 10. Root mean square fractional uncertainty in H0 (y-axis) is plotted as a function of the number of halos divided by the longest time delay, which is proportional to the lensing volume divided by the longest time delay (x-axis). The color labeling of the points is the same as in Fig. 9. The points show the fractional uncertainty for each of the six systems analyzed in this work, and error bars show the 68% bootstrap confidence intervals. The dashed line shows the fitting function given by Eq. (14). We use the fitting function in Eq. (14) to estimate the contribution of substructure to the uncertainty in H0 for the lens system SDSS 1206+4332 (purple star).

This simple relation captures the leading-order effects of substructure, but there are additional factors that could affect the uncertainties and produce scatter around the $\frac{\sqrt{V}}{t_{\max }}$ curve. For example, Oguri (2007) noted that lenses with symmetric image configurations are more sensitive to gravitational potential perturbations than lenses with asymmetric image configurations. Also, the Sheth-Tormen halo mass function in the range 10⁶ − 10⁹ M_⊙ has a mild redshift dependence that we do not account for in the fitting function.

We can use Eq. (14) to estimate the fractional uncertainty in H₀ given the lens redshift, the source redshift, the Einstein radius, and the longest time delay. For the doubly imaged system SDSS 1206+4332, Eq. (14) gives an uncertainty of 1.0%. This system was not analyzed explicitly in this work since the analysis pipeline is designed for quadruple-image lenses, but the physical arguments that motivate the scaling $\frac{\sqrt{V}}{t_{\max }}$ should still apply to SDSS 1206+4332.

We summarize the total error budget for each lens system, including the additional source of uncertainty from dark substructure, in Table 1. For a joint inference from a sample of lenses, we compute the contribution to the error budget from substructure σ_net, given by ${\sigma }_{\text{net}}^2={\left(\sum {\sigma }_i^{-2}\right)}^{-1}$, where σ_i are the uncertainties from substructure of each lens system included in the joint inference. Based on the results for the lenses we explicitly model summarized in Table 1, plus the 1.0% uncertainty estimated using Eq. (14) for SDSS 1206+4332, we find σ_net = 0.5%. This is an independent source of uncertainty that should be summed in quadrature with other contributions.

5. Summary

We have carried out an analysis of how dark matter subhalos and field halos along the line of sight affect inferences of the Hubble constant H₀ obtained through time delay cosmography. Using mock lenses based on six quadruply imaged quasars modeled by the TDCOSMO collaboration, we quantify on a case-by-case basis how dark substructure would have affected the precision and accuracy of cosmographic inferences from each system, had substructure been included in the model. Our main results are summarized as follows:

• We show that excluding dark substructure from lens models used to infer cosmological parameters does not bias inferences on H₀ at the level of the 0.3%, the numerical and statistical precision of our work. This rules out the possibility that explicitly modeling substructure in cosmographic inferences would alleviate the H₀ tension between early and late-Universe probes. It also suggests that substructure cannot explain the anticorrelation between the lens redshift and the inferred time delay distance noted by Wong et al. (2020).

• Substructure perturbs the image arrival time delays and the global mass model, which translates into uncertainties on the value of H₀ derived from the time delays and the Fermat potential predicted by the lens model. The additional uncertainties from substructure listed in Table 1 range from 0.7% in the case of RXJ1131−1231 to 2.4% in the case of HE0435−1223. This source of uncertainty is an independent contribution to the error budget that should be added in quadrature with other contributions. This addition source of uncertainty does not significantly alter the current error budget per lens, but it will not be completely negligible as the precision per system improves, e.g. with the addition of spatially resolved kinematics (Shajib et al. 2018). Combining our results for the uncertainties in single lenses, the overall contribution of substructure to the error budget of the TDCOSMO sample, which should be added in quadrature with the other sources of uncertainty, is 0.5%.

• We show that the perturbation to H₀ scales as the square root of the lensing volume divided by the longest time delay. We provide a simple fitting function that uses this scaling to predict the level of perturbation for a strong lens system given the lens and source redshifts, the Einstein radius, and the longest delay.

6. Discussion and conclusions

Our results are consistent with previous analysis that show H₀ inferences are internally consistent (Millon et al. 2020b; Wong et al. 2020), and are consistent with cosmological distances inferred from Type Ia supernova (Pandey et al. 2020). Since the level of perturbation from line-of-sight halos evolves with redshift, if un-modeled halos biased H₀ inferences significantly it would lead to inconsistencies between lenses at different redshifts, and between strong-lensing based inferences on H₀ with independent probes of the distance ladder.

While we find no evidence for bias, we conclude that the omission of substructure from strong lens models results in a slight underestimate of the uncertainties in the measured values of H₀. The additional uncertainty term, which should be summed in quadrature with the other contributions to the error budget, ranges from less than one percent in the case of RXJ1131−1231, to 2.4% in the case of HE0435−1223, increasing the total uncertainties in these systems from 4.3% to 4.4% and 6.5% to 6.9%, respectively. When inferences from the seven lenses analyzed by TDCOSMO are combined, we estimate that substructure contributes an additional 0.5% uncertainty. This is a subdominant contribution to the overall error budget that we reiterate should be summed in quadrature with other contributions.

To make use of the results presented in this paper for future cosmographic inferences, we suggest two possibilities. First, one could obtain a rough estimate of the level of perturbation using the fitting function plotted in Fig. 10, and possibly incorporate this additional uncertainty when quoting results on H₀. Second, the analysis performed in this work could be repeated for new lens systems. The first approach sacrifices precision for speed, while the second approach is more time consuming, but gives a more precise estimate for the additional uncertainty from substructure.

Finally, we note that dark matter halos should also alter the arrival times between images of lensed supernova, such as supernova “Refsdal”, which is multiply imaged by an early-type galaxy in the cluster MACS J1149.5+2223 (Kelly et al. 2015). We can extrapolate the results of this work to the Refsdal system to approximate the impact of substructure on future cosmographic inferences. The Refsdal system consists of four images arranged in an Einstein cross configuration separated by ∼1″ with a maximum delay of 24 days (Rodney et al. 2016), and a fifth image ∼30″ away with a time delay of ∼1 yr. Coincidentally, both pairs of Einstein radii and time delays give an uncertainty of 2.1% in H₀ from substructure, per Eq. (14), because the long time delay between S_X and the Einstein cross counteracts the larger image separation.

Acknowledgments

We thank Anowar Shajib for answering many questions, and for sharing the MCMC chains computed for DESJ0408-5354. We also thank Adriano Agnello, Frederic Courbin, and Dominique Sluse for comments on the manuscript, and the anonymous referee for constructive feedback. DG, TT, and SB acknowledge support by the US National Science Foundation through grant AST-1714953 and AST-1906976. This work used computational and storage services associated with the Hoffman2 Shared Cluster provided by the UCLA Institute for Digital Research and Education’s Research Technology Group.

References

Appendix A: Baseline macro lensmodels

Table A.1 lists the parameters describing the lens macromodels for the mock lens analogs of the six strong lens systems considered in this work.

Appendix B: Hiding dark matter halos in extended images with shapelets

As a byproduct of this analysis, we obtain hundreds of mock lenses with realistic line-of-sight and substructure populations that we model with state of the art image reconstruction techniques. These simulations are relevant to the field of gravitational imaging of dark substructure in lensed arcs (Vegetti et al. 2010, 2018; Hezaveh et al. 2016a). Dark substructure occasionally produces imaging residuals in lensed arcs that betray the presence of individual halos, or provide a means to infer quantities such as the convergence power spectrum of halos in the lens system (Hezaveh et al. 2016b; Diaz Rivero et al. 2018; Cyr-Racine et al. 2019; Çağan Şengül et al. 2020).

As the convergence itself is not observable, it must be derived from the surface brightness residuals from the lens modeling (Cyr-Racine et al. 2019). Keeping all else fixed, adding additional complexity to the background source light partially removes gravitational imaging residuals in extended images. The gravitational imaging residuals therefore depend explicitly on how the unknown background source is modeled.

To illustrate, the two rows in Fig. B.1 show examples of imaging residuals computed for the mock analogs of B1608+656 and PG1115+080. The far left panels in Fig. B.1 show the imaging residuals without shapelets added to the background source, and the second from left panels show the same mock lens system, with the same population of dark matter halos, modeled with shapelets included in the background source. The population of line-of-sight halos and subhalos that perturbs the mock images are shown in the third panels, using the definition of the multi plane convergence described in Sect. 4.1. The far right panel shows the difference between the source surface brightness between the models that include and exclude shapelets. In these two cases, adding shapelets to the background source removes imaging residuals from the extended images.

Fig. B.1. Examples of imaging residuals, convergence, and residual source surface brightness for one mock lens system based on B1608+656 (top row) and PG1115+080 (bottom row). The far left panel shows the gravitational imaging residuals when only an elliptical Sérsic profile, the same profile used to create the mock data, is used to model the background source. The second from left panel shows the imaging residuals when shapelets are added to the source light profile. The dark matter halos that perturb the lensed arc in the mock data are shown in the third from left panel, and the residual source surface brightness is shown in the far right panel.

In Fig. B.2, we show the stacked power spectra of imaging residuals from 100 different realizations. The shaded regions represent show 68% bootstrap confidence intervals around the median in each distance bin plotted on the x-axis. Without shapelets included in the background source, substructure leads to enhanced power on small scales, as expected. The addition of shapelets, which in this example act as a proxy for any additional source complexity that is typically required to model lensed images, removes the signal.

Fig. B.2. The stacked power spectrum of gravitational imaging residuals computed from 100 residual maps like those in Fig. B.1. The manner in which substructure affects the power spectrum of imaging residuals on small scales depends on how the background source is modeled.

While the effects of very massive halos that produce direct hits on extended images evidently cannot be completely absorbed by the background source (e.g., Vegetti et al. 2010; Hezaveh et al. 2016a), the collective effects of the many smaller halos that affect the convergence power spectrum can be partially absorbed by a reconfiguration of the lens and source model. The degree to which this occurs depends on the specifics of the source model, and on the quality of the imaging data. In principle, the fact that the source is multiply imaged breaks degeneracies between the source structure and imaging residuals (Cyr-Racine et al. 2019), but with finite resolution and signal to noise the effect persists. Higher-resolution images, targeting highly magnified portions of extended images, and the identification of what scales are dominated by substructure perturbations rather than source morphology, may help deal with this systematic challenge.

Appendix C: Joint posterior distributions of H₀ with lens macromodel parameters

Figures C.1–C.6 show the joint posterior distributions of macromodel parameters with the Hubble constant for the six lens systems analyzed in this work.

Fig. C.1. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the control dataset (black) and the substructure-perturbed dataset (blue). The parameters are (reading the x-axis labels from the left) the Einstein radius, the logarithmic slope of the main deflector, the ellipticity position angle, the ellipticity, the external shear position angle, the external shear strength, and the Hubble constant. Vertical lines show 95% confidence intervals and the contours show 68% and 95% confidence intervals. The distributions show are computed for the lens system RXJ1131−1231.

Fig. C.2. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system PG1115+080. See Fig. C.1 for a description of the parameters.

Fig. C.3. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system HE0435−1223. See Fig. C.1 for a description of the parameters.

Fig. C.4. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system B1608+656. See Fig. C.1 for a description of the parameters.

Fig. C.5. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system WFI2033−4723. See Fig. C.1 for a description of the parameters.

Fig. C.6. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system DESJ0408−5354. See Fig. C.1 for a description of the parameters.

Appendix D: Joint posterior distributions of time delays and flux ratios

The initial motivation for considering the effects of substructure on time delays was to probe the nature of dark matter through these data. The thinking was that time delay “anomalies”, or measured time delays that differ from the model-predicted time delays, encode the abundance and spatial distribution of dark substructure (e.g., Keeton & Moustakas 2009; Cyr-Racine et al. 2016).

The Shapiro delay associated with dark (sub)halos depends on the projected gravitational potential of halos, and is a long-range interaction. This could, in principle, complement the information encoded by image flux ratios, which are highly localized probes of dark matter structure and provide an avenue to probe small-scale structure in strong lenses (Dalal & Kochanek 2002; Nierenberg et al. 2014; Hsueh et al. 2020; Gilman et al. 2020a,b). Figures D.1–D.6 show the joint distributions of the residual time delays and residual flux ratios for the six lens systems analyzed. The flux ratios are computed with a background source parameterized as a Gaussian with a FWHM of 25 pc.

Fig. D.1. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system RXJ1131−1231.

Fig. D.2. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system PG1115+080.

Fig. D.3. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system HE0435−1223.

Fig. D.4. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system B1608−656.

Fig. D.5. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system WFI2033−4723.

Fig. D.6. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system DESJ0408−5354.

All Tables

All Figures

Fig. 1. Mock observed lens that includes substructure in the lens model (far left), the reconstructed lensed image (second from left), the normalized imaging residuals (center), the residual time delay surface (second from right), and the residual convergence (far right). The residual time delay surface is the true time delay surface from the mock lens minus the model-predicted time delay surface from the lens model fit to the data. Similarly, the residual convergence is the full multi plane convergence (Eq. (10)) minus the convergence from the model fit to the mock data. The lensed quasar image positions are marked as points (crosses) in the center (right) panels. Each row corresponds to one substructure-perturbed dataset for the lens systems based on RXJ1131−1231, PG1115+080, and HE0435−1223.

In the text

	Fig. 2. The figure is the same as Fig. 1, but shows examples of the mock observed and reconstructed lens (left), the imaging residuals (center), the residual time delay surface (second from right), and the residual convergence (far right) maps for the lens systems B1608+656, WFI2033−4723, and DESJ0408−5354.
In the text

Fig. 3. The joint posterior distribution of the residual time delays Δt, defined as the model-predicted time delay minus the true time delay, and the Hubble constant, computed for the lens system RXJ1131−1231. The contours show the 68% and 95% confidence intervals for the control posterior that includes no substructure (black), and the combined posteriors from the inferences on lens models that include substructure (blue). Dashed lines in the histograms show 95% confidence intervals. The median and 68% confidence intervals of the H0 inference are printed above the figure.

In the text

	Fig. 4. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system PG1115+080.
In the text

	Fig. 5. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system HE0435−1223.
In the text

	Fig. 6. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system B1608+656.
In the text

	Fig. 7. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system WFI2033−4723.
In the text

	Fig. 8. Same as Fig. 3, but showing the posteriors of the residual time delays and H0 computed for the lens system DESJ0408−5354.
In the text

Fig. 9. The y-axis shows the bias in H0, defined as the median of the H0 posterior obtained from the substructure-perturbed datasets minus the median of the posterior from the control datasets. The mock lenses are sorted on the x-axis by the square root of the lensing volume divided by the longest time delay, a dimensional quantity that determines the additional uncertainty caused by substructure in the lens system. The error bars show the statistical uncertainty on the median, which is given by $\frac{{\sigma }_{\mathrm{control}}}{\sqrt{N}}$, where σcontrol is the uncertainty in H0 from the control dataset and N = 200 is the number of realizations. We find no evidence for bias in the inferred value of H0 stemming from substructure that is not modeled in cosmographic inferences.

In the text

Fig. 10. Root mean square fractional uncertainty in H0 (y-axis) is plotted as a function of the number of halos divided by the longest time delay, which is proportional to the lensing volume divided by the longest time delay (x-axis). The color labeling of the points is the same as in Fig. 9. The points show the fractional uncertainty for each of the six systems analyzed in this work, and error bars show the 68% bootstrap confidence intervals. The dashed line shows the fitting function given by Eq. (14). We use the fitting function in Eq. (14) to estimate the contribution of substructure to the uncertainty in H0 for the lens system SDSS 1206+4332 (purple star).

In the text

Fig. B.1. Examples of imaging residuals, convergence, and residual source surface brightness for one mock lens system based on B1608+656 (top row) and PG1115+080 (bottom row). The far left panel shows the gravitational imaging residuals when only an elliptical Sérsic profile, the same profile used to create the mock data, is used to model the background source. The second from left panel shows the imaging residuals when shapelets are added to the source light profile. The dark matter halos that perturb the lensed arc in the mock data are shown in the third from left panel, and the residual source surface brightness is shown in the far right panel.

In the text

	Fig. B.2. The stacked power spectrum of gravitational imaging residuals computed from 100 residual maps like those in Fig. B.1. The manner in which substructure affects the power spectrum of imaging residuals on small scales depends on how the background source is modeled.
In the text

Fig. C.1. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the control dataset (black) and the substructure-perturbed dataset (blue). The parameters are (reading the x-axis labels from the left) the Einstein radius, the logarithmic slope of the main deflector, the ellipticity position angle, the ellipticity, the external shear position angle, the external shear strength, and the Hubble constant. Vertical lines show 95% confidence intervals and the contours show 68% and 95% confidence intervals. The distributions show are computed for the lens system RXJ1131−1231.

In the text

	Fig. C.2. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system PG1115+080. See Fig. C.1 for a description of the parameters.
In the text

	Fig. C.3. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system HE0435−1223. See Fig. C.1 for a description of the parameters.
In the text

	Fig. C.4. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system B1608+656. See Fig. C.1 for a description of the parameters.
In the text

	Fig. C.5. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system WFI2033−4723. See Fig. C.1 for a description of the parameters.
In the text

	Fig. C.6. Joint distribution of selected macromodel parameters and the inferred Hubble constant for the lens system DESJ0408−5354. See Fig. C.1 for a description of the parameters.
In the text

	Fig. D.1. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system RXJ1131−1231.
In the text

	Fig. D.2. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system PG1115+080.
In the text

	Fig. D.3. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system HE0435−1223.
In the text

	Fig. D.4. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system B1608−656.
In the text

	Fig. D.5. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system WFI2033−4723.
In the text

	Fig. D.6. Joint distribution of flux ratio and time delay residuals (model minus data) computed for the lens system DESJ0408−5354.
In the text