TDCOSMO XIV: Practical Techniques for Estimating External Convergence of Strong Gravitational Lens Systems and Applications to the SDSS J0924+0219 System

Time delay cosmography focuses on strong lensing of time-variable sources, usually a quasar. As the luminosity of the source varies, this variation appears in each of the several images of the source, but not at the same time. The time delay between any two images can be written as follows:

where $\overrightarrow{\theta}$ represents the angular position of an image on the sky, $\overrightarrow{\beta}$ represents the actual (unobservable) angular position of the source on the sky, and $\psi$ is the scaled lensing potential. The first two terms are a result of the different distances traveled by the light for the two images, while the later two terms are the difference in the Shapiro time delay. The goal of a complete cosmographic analysis is to infer the time-delay distance (see eq. 1). Using measurements of the time delay combined with a robust model of the lensing galaxy it is possible to measure the time delay distance and, therefore, the Hubble constant.

The multiple images observed in such a system are en example of strong gravitational lensing. Qualitatively, strong lensing is any lensing which produces multiple images of some background object. Quantitatively, strong lensing occurs whenever the local density of the lens is greater than the lensing critical density:

For a given perturber, the convergence at a given location is defined as the local density in units of the critical density:

For $\kappa<1$, strong lensing does not occur. Instead, mass distributions with $\kappa<1$ result in magnification or demagnification of the image of the background object. This is the situation for the perturbers along the line of sight in our system that are not included in the primary lens model, though we note that some perturbers that meet this criterion are included in the mass model based on their flexion shift (see section 2.2), and therefore excluded from our statistical analysis. However, the effect of the pertubers we do include in the estimate discussed here is not directly observable because the actual angular size of the background object is not known.

We quantify the cumulative effect of all perturbers (not including those incorporated into the primary lens model) and voids along the line of sight with the external convergence, $\kappa_{\rm{ext}}$. Conceptually, the value of $\kappa_{\rm{ext}}$ is the convergence of a mass sheet which, if placed coplanar to the primary lens, would produce the same magnification or demagnification as the perturbers do collectively. The external convergence is defined relative to that of a line of sight where the mass distribution is smoothly distributed with a density equal to the global mass density of the Universe. Therefore positive values of $\kappa_{\rm{ext}}$ represent lines of sight that are overdense with respect to the overall density of the Universe, while underdense lines of sight have negative values.

These perturbers can also produce shear, denoted $\gamma_{\rm{ext}}$, which results in stretching and distorting of the image of the source. However, this effect can be estimated in the primary lens model, as it affects the location of the the images and the shape of the Einstein ring (if present). We use the $\gamma_{\rm{ext}}$ constraint from the lens model in our analysis as an additional constraint (see section 2.6)

Assuming $\kappa_{\rm{ext}}$ can be measured, it serves as a correction factor to the computed value of the time delay distance:

where $D^{\prime}_{\Delta t}$ represents the uncorrected value of the time delay distance. This propagates directly to the inferred value of the Hubble constant by

We now review the techniques we use to compute the value of $\kappa_{\rm{ext}}$ for a generic lens system.

To start, we identify objects along the line of sight that contribute to $\kappa_{\rm{ext}}$. We separate objects into strong and weak perturbers, using an operational definition discussed below. Strong perturbers are generally close to the center of the field or agalaxy group (see for example Sluse et al., 2019; Fassnacht et al., 2006). Strong perturbers are included explicitly in the mass model of the lens, while weak perturbers are treated statistically. To separate these, we use the flexion shift formalism, first proposed in McCully et al. (2017). The flexion shift of an object is given by

where $\theta_{E}$ and $\theta_{E,p}$ are the Einstein radii of the main lens and perturber respectively, and $\theta$ is their angular separation. The quantity $f(\beta)$ is given by

where

and where $D_{\rm{dp}}$, $D_{\rm{s}}$, $D_{\rm{p}}$, and $D_{\rm{ds}}$ are the angular diameter distances from the deflector to the perturber, to the source, to the perturber, and from the deflector to the source. The flexion shift roughly measures the perturbations to the images of the source due to 3rd order terms from the perturber. Ultimately, what constitutes a ”weak” or ”strong” perturber is somewhat arbitrary, but there this is a clear trade-off between the improvement from including a particular perturber in the mass model and the amount of work required to do so. McCully et al. (2017) recommends using $\Delta_{3}x=10^{-4}\text{arcsec}$ as the cutoff between strong and weak perturbers to ensure a $<1\%$ bias on $H_{0}$.

The first step in determining the value of $\kappa_{\rm{ext}}$ is comparing the field of interest to a large number of lines of sight from some large reference field. This comparison gives us an empirical estimate of the relative matter density of the lens field as compared to all lines of sight in the universe. The reference field should be large enough to avoid sampling bias. For a small comparison field (on the order of a few $deg^{2}$) statistical overdensities or underdensities may occur (see for example Fassnacht et al., 2010). However modern survey datasets are available which cover hundreds to thousands of square degrees, allowing us to use a sufficiently large comparison field to avoid sampling bias (see section 4.2 for a discussion of our choices for this analysis). Ideally the data for the reference field and lens field are taken by the same instrument and processed by the same analysis pipeline. This turns out to be the case for the analysis of J0924 discussed later in this paper, but will not be true in general. At a very minimum we seek data with at least one band in common, with deep enough observations to produce meaningful results. When comparing the survey dataset to our line of sight, we set a magnitude cut and remove any objects from both datasets fainter than this cut. Previous work (see Collett et al., 2013) has suggested that $i<24$ represents a good limit, as setting a fainter limit does not appear to meaningfully impact the results of this style of analysis. This limit is also bright enough to be well above the detection limits of modern sky surveys, ensuring reliable photometry.

Once the appropriate data are in hand, it is important to consider which objects should be used in the comparison. It is important to set a magnitude cut that will remove objects too close to the detection limit of the instrument for photometry to be reliable, while still leaving enough data to make robust estimates of the relative density. However setting too bright of a limit will result in having too few objects to compare to.

Additionally, we cut out all objects with a redshift greater than the redshift of the source quasar, as these objects will not affect the path of the light as it travels from the quasar to our telescopes.

Because the value of $\kappa_{\rm{ext}}$ cannot be directly measured, we first define tracer quantities that can be computed directly from the available data. By comparing the value of these quantities in the lens field to the value of the identical quantity computed for a large number of reference fields, we obtain an empirical estimate of the relative density of the line of sight of interest. As a first approximation, we expect the greatest contribution to $\kappa_{\rm{ext}}$ from massive objects close to the center of the line of sight. The primary mass contribution in any line of sight will be dark matter halos, but these are not directly observable. To quantify the contribution from these weak perturbers, we compute weighted number counts of the visible structures (i.e. galaxies) in the lens field as compared to a large number of reference fields. This technique has been used extensively in previous work (Fassnacht et al., 2006; Suyu et al., 2010; Greene et al., 2013; Buckley-Geer et al., 2020) . To do this, we select a region of interest around the lens and compare it to a large number of identically-shaped fields selected at random from the reference survey. At each step, we compute the ratio of the weighted number counts for the galaxies in the lens field to the identical statistic computed in the given reference field. For the given step, the value of the weight is therefore:

Where $i$ indexes the reference fields, $j$ indexes the galaxies in a given field, and $w$ is the value of the weighted statistic for the given galaxy.

Following Rusu et al. (2017) we also consider a second style of weighting that improves our results. Instead of summing the value of weights for all objects in the field, we instead compute the weight for a given field as $w_{i,\rm{meds}}=n_{i}\overline{w_{j}}$ where $n_{i}$ is the number of galaxies in the given reference field and $\overline{w_{j}}$ is the median value of the weight for all galaxies in the reference field. Doing this helps avoid situations where single objects dominate the sum in a particular line of sight. This is especially important for weights involving stellar mass and the inverse seperation. In this scheme, the value of the weight for the given reference field is therefore

There has been discussion in the literature about which are the best weights to consider for the most robust determination of $\kappa_{\rm{ext}}$ (Greene et al., 2013; Rusu et al., 2017, 2019). In this work, we only consider a subset of the weights considered in previous works (see Table 1).

In order to determine the posterior distribution of $\kappa_{\rm{ext}}$ for the given system, it is necessary to compare the weighted number counts to similar counts obtained from a reference field for which $\kappa_{\rm{ext}}$ is known. We use a simulated dataset for this purpose. The simulation must contain several components in order to be suitable for this analysis. First, it must contain catalogs of galaxies with known luminosity and redshift. Second, values of the external convergence must be measured at a suitably large number of points to be representative of the universe at large. Rusu et al. (2017) examined possible biases from inferring $\kappa_{\rm{ext}}$ using the number counts method in the Millennium Simulation, and found these methods produced a good estimate of $\kappa_{\rm{ext}}$. We discuss our choice of simulation further in section 4.3. Because this technique involves ratios, much of the dependence on the simulation’s underlying cosmological parameters should cancel out. However, ensuring this would require a second simulation with the attributes described. An exciting development in this space are the initial results from the MillenniumTNG (Hernández-Aguayo et al., 2022). Full weak-lensing convergence maps are planned, but not yet available.

We now have weighted number counts for the lens field itself and for a large number of fields in a simulated dataset, each of which is associated with a value of $\kappa_{\rm{ext}}$. We seek to compute $p(\kappa_{\rm{ext}}|\textbf{d}$): the probability distribution of $\kappa_{\rm{ext}}$ given the data. We can replace this with the joint probability distribution of $\kappa_{\rm{ext}}$ and the data as follows:

After some work it can be shown (see Rusu et al., 2017)

Where $p(W|\textbf{d})$ is the probability distribution of given weighting scheme given the data, and $p_{\rm{sim}}(\kappa_{\rm{ext}}|W)$ is the probability distribution of $\kappa_{\rm{ext}}$ in the simulated dataset, given a particular value of the weight. Here, it is assumed that the simulated dataset is the correct prior for the observable universe. In other words, that it correctly relates values of weights to values of $\kappa_{\rm{ext}}$

Additionally, a constraint on $\kappa_{\rm{ext}}$ based on $\gamma_{\rm{ext}}$ can be included, which accounts for the expected correlation between these two values. In general, $\gamma$ is a two-dimensional vector on the plane of the sky, but we use only the overall magnitude in our analysis. $\gamma_{\rm{ext}}$ is a parameter that can be fitted in the primary mass model of the lens. We use $\gamma_{\rm{ext}}$ just as we would a weighted number count distribution. For any range of values of $\gamma_{\rm{ext}}$, we can construct a histogram of the value of $\kappa_{\rm{ext}}$ for all lines of sight with values of $\gamma_{\rm{ext}}$ that fall in this range. We then weight the contribution from these lines of sight by $p(\gamma|\mathbb{d})$, which is the posterior on $\gamma$ inferred from the mass model.This provides a prior on $\kappa_{\rm{ext}}$ that can meaningfully affect the final result. In section 5, we present results both with and without using $\gamma_{\rm{ext}}$ as a constraint.

In previous work (e.g. Rusu et al., 2017), this full probability distribution for each weight $p(W|\textbf{d})$ was replaced by a normal distribution centered on the median of the full weight distribution, with a width determined by examining measurement uncertainties. In practice, this width was much smaller than the width of the actual distribution. This is cheaper computationally, but ignores covariance between the various types of the weights. While the increase in computation time is significant, the majority of the important decisions in the analysis are made when we compute weighted number count ratios, and using this method does not meaningfully increase our time-to-result

Since all the individual weights are measured at the same sequence of randomly-drawn fields, we can construct a full m-dimensional probability distribution $p(W_{s}|\textbf{d})=p(W_{n},W_{1/r},...|\textbf{d})$ where m is the number of weights being considered. We then explore this probability distribution when implementing the formalism described above. Formally, the posterior on $\kappa_{\rm{ext}}$ becomes:

When computing this quantity, we split the m-dimensional probability distribution into $200^{m}$ m-dimensional bins. We have also tested this procedure $100^{m}$ bins and see consistent results for the J0924 field. With significantly more bins, the computational time balloons and the number of fields in each bin drops significantly, even near the center of the distribution. The value of $p(W_{s}|\textbf{d})$ is simply the number of lines of sight in the reference survey that fall in this bin divided by the total number of lines of sight being considered. Given the large numbers of lines of sight we consider, the distributions are smooth and it may be possible to model them explicitly and explore the distribution with an MCMC, but we save this for a future analysis. The exception to this is the pure unweighted number counts in a $45^{\prime\prime}$ aperture, due to the fact that the ”weight” is integer valued and the number of galaxies in the aperture is relatively modest.

For $p_{\rm{sim}}(\kappa_{\rm{ext}}|W_{s})$ we construct a histogram of the measured values of $\kappa_{\rm{ext}}$ for all lines of sight from the simulated datasets that fall within the given bin, normalized by the number of lines of sight in the bin. Without this normalization, $\kappa_{\rm{ext}}$ values near the mean of the simulation will always be weighted more heavily because there are comparatively more of them. Our choice of simulation is discussed in section 4.3

The inputs to our analysis code include the full weight distributions for the lens field, the weights computed from the simulated dataset, and the $\kappa_{\rm{ext}}$ maps for the simulated dataset for the source redshift. We iterate over the probability distribution discussed above, computing a histogram of the values of $\kappa_{\rm{ext}}$ in each m-dimensional bin. The overall histogram is the sum of these histograms, weighted by the value of the weight distribution in that bin.

There are other methods for estimating $\kappa_{\rm{ext}}$ which have been explored as discussed below. Ultimately many of these techniques involve a trade-off between speed of analysis and precision of the final result. Our goal is explicitly to design and implement techniques that allow us to analyze hundreds or even thousands of lens systems in a reasonable amount of time, with the goal of combining full results (including lens modeling, and time delay measurements) from many lenses to make some final statement about the value of $H_{0}$. Because of this, a technique with slightly less constraining power for a single lens is tolerable if it can be performed and iterated on rapidly. A few other techniques for estimating $\kappa_{\rm{ext}}$ are discussed below. While all show promise, and are interesting for their own sake, none show a significant enough improvement to justify the increased time and complexity of analysis, at least in the context of our stated large-scale goals.

One of the big challenges in doing analyses on these kinds of datasets is the need for the computing environments to be close to the data whenever possible. Querying over the internet is useful for assembling datasets, but is not a particularly good solution when analyzing a dataset at scale. Many researchers will not have access to sufficient storage to store these datsets, and it is impractical to expect individual survey teams to provide computing resource for general use. The size of these datasets will enable next-generation analyses, but only with the development of next-generation tools running at scale, which will require computing infrastructure that may not be readily available to many researcher.

We support using cloud computing services to fill this gap. Commercial providers have expanded and matured by leaps and bounds over the last decade, and routinely handle storage and analysis tasks on datasets orders of magnitude larger than the ones being discussed here. Additionally, cloud computing technologies are significantly more accessible then on-premise technology: they require far lower startup costs and can be quickly scaled (to accommodate more users, or bigger jobs) without the bottlenecks that slow down the expansion of on-premise infrastructure.

It is also possible to use cloud solutions as a supplement to already-existing on premise solutions. Holzman et al. (2017) demonstrated this by analyzing data from the Compact Muon Solenoid experiment at massive scale.

In the future, we plan to develop lenskappa and heinlein tools that could be easily deployed onto services like these, providing quick and easy access to large survey datasets in addition to techniques for processing and analyzing that data. The datasets would be stored in the cloud, allowing users to deploy their analyses without worrying about the connection to the underlying dataset. Such an approach has been demonstrated by Kusnierz et al. (2022), which enabled serverless access to ROOT for quick analysis tasks on high energy datasets without the user having to manually retrieve the data. With some work, heinlein will enable this type of analysis by serving as a bridge between computing infrastructure (which could be managed by individual researchers) and the underlying data lake (which could be managed by the survey team). Abstracting away data retrieval will allow researchers to focus on what matters most: designing the analyses they wish to perform and interpreting results.

SDSS J0924+0219 is a quadruply lensed quasar first reported in Inada et al. (2003). The quasar itself is at redshift $z=1.523$, while the lensing galaxy is located at $z=0.384$ (Eigenbrod, A. et al., 2006). Quadruply lensed quasars are particularly valuable for cosmographic analysis because it is in principle possible to measure 12 different time delays, though only three of these are independent. An image of the field can be seen in figure 2. This system was modeled in Chen et al. (2022)

One particularly important feature of the lens field is the bright star located in the lower left. The star covers a reasonably large fraction of the overall area in the 120” aperture, undoubtedly covering several background objects and making photometry for objects very near it unreliable. We apply techniques that have previously been used in this kind of analysis to correct for this. This technique is discussed in more detail in section 4.2.

The lens field falls within the Subaru Hyper Surpime-Cam Strategic Survey Program footprint, and full color information is available for all objects in the relevant catalogs. We use these catalogs, including photometric redshifts, for all objects inside the field. Additionally we use the bright star masks and photometric redshift PDFs provided by the HSC team.

Our analysis follows the same outline discussed in section 2. We discuss the details particular to this lens field below.

The Hyper Suprime-Cam Subaru Strategic Survey is a large survey program, aiming to cover roughly 1400 deg² of sky in five photometric bands (grizy) down to $i\sim$ 26, with deeper coverage expected in smaller regions of the sky (Aihara et al., 2017). We base our analysis on the roughly 400 deg² that had coverage in all five bands as of the second data release (Aihara et al., 2019). Data release 3 was made available while this paper was in preparation, but we do not consider it here.

The HSC Survey is a natural choice for a comparison dataset for this system because the J0924 field itself falls within the survey footprint. We therefore have both robust and comparable photometry. When computing weighted number count ratios, the basic approach is identical to the one outlined in the section 2.4, with some specific adjustments:

We use objects with r_extendness_value $=1.0$, which selects galaxies. Bosch et al. (2017) demonstrates that their algorithm for computing this quantity does a reasonably good job of selecting galaxies, though it may incorrectly classify some galaxies as stars, and vice-versa. However we note that the HSC-wide survey was intentionally selected to enable cosmological analyses by being in regions that are away from the galactic plane and low on dust extinction (Aihara et al., 2017). This, combined with robust bright-star masks and our large sample size ensures unmasked stars do not significantly impact the final weighted number count ratios.

The full photometric redshift PDFs for all objects in the data release have been made available by the HSC team (Nishizawa et al., 2020). They use two separate fitting algorithms, and results from both algorithms are included in the catalog. We compute weighted number count ratios using both sets of redshifts, and do not find a meaningful difference between the resultant distributions. As such, we use the ”DEmP” redshift and stellar masses for our analysis (Hsieh & Yee, 2014; Tanaka et al., 2017).

When computing weighted number counts, we remove all objects closer than 5 arcsec from the center of the field following Rusu et al. (2017). Objects this close to the center of the field are typically explicitly included in the mass model of the lens, and so we also remove them from the comparison fields to avoid biasing results.

The HSC survey team makes available masks that represent areas of the sky where photometry may be unreliable or lacking due to the presence of bright stars (Coupon et al., 2017). This is particularly important in our field due to the presence of the bright star that can be seen in Figure 2 When iterating over the reference survey, we retrieve the bright star masks for each region being considered. We apply both these masks and the masks for the lens field itself to both catalogs at each weighting step. Doing this ensures that results are not biased if a given field in the reference survey has significantly more or less of its area covered by bright star masks. This procedure was first used in Rusu et al. (2017)

The 400 deg² of sky we consider here is separated into seven disconnected regions. Initially, we computed weighted number count ratios in five of these regions. This produces nearly-identical distributions, with the difference between the lowest and highest median being $0.1\sigma$ the standard deviation of the distribution, suggesting our fields are large enough to to avoid sampling bias. Based on this, we restrict our subsequent analysis to 135 deg² of sky located in the region $332^{\circ}<\text{RA}<359^{\circ}$ and $-1.5^{\circ}<\text{Dec}<6^{\circ}$. For each combination of aperture and limiting magnitude, we compute weighted count ratios at 100,000 randomly selected fields.

The Millennium Simulation (Springel et al., 2005) is a dark matter only simulation split into 64 $4\times 4\deg^{2}$ fields. After the original run was completed, synthetic galaxy catalogs were painted into the resultant halos by several teams. Following (Rusu et al., 2017) we use the semi-analytic catalogues of De Lucia & Blaizot (2007). Additionally, Hilbert et al. (2009) split each $4\times 4\deg^{2}$ field into a grid of 4096x4096 points and used ray tracing to compute convergence and shear at each of these points in 63 redshift planes. These, combined with its large size, makes it an excellent choice for our analysis. In our analysis, we use redshift plane 36 with $z=1.504$.

First, we compute the weighted number counts at a large number (order $10^{6}$) of equally spaced grid points in the Millennium simulation. For the $45"$ aperture, we place the fields $90"$ apart (snapped to the nearest grid point), while for the $120"$ aperture we place fields $60"$ apart. Both cases result in over 1.5 million lines of sight across the simulation, each of which has an associated value of $\kappa_{\rm{ext}}$ and $\gamma_{\rm{ext}}$ This differs from the same calculation for the lens field itself in a key way: the values reported are the total value of the weights at every point considered, rather than a ratio of values. To normalize, we divide weighted number count in each field by the median value for all lines of sight in our sample. The median value of the resultant distribution is therefore unity.

A key difference between the synthetic catalogs and the real data from the HSC survey is that in the Millennium Simulation catalog redshifts are exact. In previous work (eg. Rusu et al., 2017) this difference was accounted for by computing photometric redshifts for all objects in the Millennium Simulation, using the same pipeline that was used to compute the photometric redshifts in the comparison dataset. We are unable to use this technique in this case because the HSC survey photo-z pipelines are not publicly available as of the preparation of this paper. Instead, we download the full catalog of training data used by the HSC. These are galaxies for which spectroscopic redshifts are available. We divide the objects in the test dataset into redshift and magnitude bins. For each galaxy in a bin, we compute the offset in the central value of the redshift $(z_{phot}-z_{spec})$, and take the median of these as our estimate of the redshift bias in that bin. Additionally we take the median value of $\sigma_{z}$ (as reported by the HSC photometric redshift pipeline) for all photometric redshifts in that bin. We take the the median value of $\sigma_{z}$ as our estimate of photometric redshift uncertainty. For each object in the Millennium Simulation catalog, we construct an artificial photemetric redshift PDF. The center of the distribution is the ”actual” redshift of the object, offset by the the amount computed previously for the appropriate bin. The width of the distribution is the value of $\sigma_{z}$ computed for the same bin.

At each weighing step in the Millennium Simulation, we then sample from these ”photometric redshift” PDFs when computing weighted number counts. For each line of sight, we sample from each of these PDFs 50 times. This produces 50 separate catalogs (all with slightly different values for object ”redshifts”), and we compute weighted number counts for each one. We find that this process produces no meaningful change to our final inference for $\kappa_{\rm{ext}}$, even when a weight that depends on the redshift is taken into account. However this process massively increases the amount of data output by our code, and therefore significantly increases the amount of time required to do the final inference on $\kappa_{\rm{ext}}$. This does suggest that photometric redshift uncertainties do not have a significant impact on the inferred value of $\kappa_{\rm{ext}}$, but we plan to explore this more completely in the future.

Our primary goal for this project has been to build a software tool that can quickly and reliably analyze weak perturbers along lines of sight to strong gravitational lenses. We have accomplished this goal, but we plan to extend the capabilities of lenskappa to include tools for analyzing strong perturbers and coherent structures (such as galaxy groups). These objects must be handled individually, and require very different analysis tools. However we see a significant advantages to being able to do full environment analyses within a single software package.

Leveraging the capabilities of lenskappa, we hope to better understand lines of sight to gravitational lenses on a population level. Fassnacht et al. (2010) and Wong et al. (2018) have shown that lenses seem to fall in preferentially overdense lines of sight, but that this overdensity seems to be confined to the immediate surroundings of the lens itself. Lenskappa gives us the tools to perform these population-level analyses quickly, with the freedom to adjust and re-run the analysis as needed. As a first step, we hope to complete an analysis with 4-5 new lens systems, as well as 2-3 systems that have previously been analyzed to check our results.

As an additional check, we would like to analyze a large number of non-lens lines of sight. This would ensure there are no biases introduced by comparing distributions based on real galaxy catalogs to those obtained from the synthetic catalogs in the Millennium Simulation.

Recently, Park et al. (2022) demonstrated the use of Bayesian Graph Neural Network to estimate the value of $\kappa_{\rm{ext}}$ in a simulated dataset. Their method out-performs a simplified version of the analysis performed here that uses only a single weight. Further work may demonstrate the ability of the technique to match or even outperform the weighted number counts technique, but a more complete comparision will need to be performed to assess this.