Issue 
A&A
Volume 667, November 2022



Article Number  A86  
Number of page(s)  17  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202244324  
Published online  09 November 2022 
TDCOSMO
VIII. A key test of systematics in the hierarchical method of timedelay cosmography
^{1}
STAR Institute, Quartier Agora, Allée du Six Août, 19c, 4000 Liège, Belgium
email: mgomer@uliege.be
^{2}
Kavli Institute for Particle Astrophysics and Cosmology and Department of Physics, Stanford University, Stanford, CA 94305, USA
^{3}
SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
^{4}
Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA
^{5}
Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
Received:
22
June
2022
Accepted:
1
September
2022
The largest source of systematic errors in the timedelay cosmography method likely arises from the lens model mass distribution, where an inaccurate choice of model could in principle bias the value of H_{0}. A Bayesian hierarchical framework has been proposed which combines lens systems with kinematic data, constraining the mass profile shape at a population level. The framework has been previously validated using a small sample of lensing galaxies drawn from hydrosimulations. The goal of this work is to expand the validation to a more general set of lenses consistent with observed systems, as well as confirm the capacity of the method to combine two lens populations: one which has time delay information and one which lacks time delays and has systematically different image radii. For this purpose, we generated samples of analytic lens mass distributions made of baryons+dark matter and fit the subsequent mock images with standard powerlaw models. Corresponding kinematics data were also emulated. The hierarchical framework applied to an ensemble of timedelay lenses allowed us to correct the H_{0} bias associated with model choice to find H_{0} within 1.5σ of the fiducial value. We then combined this set with a sample of corresponding lens systems which have no time delays and have a source at lower z, resulting in a systematically smaller image radius relative to their effective radius. The hierarchical framework has successfully accounted for this effect, recovering a value of H_{0} which is both more precise (σ ∼ 2%) and more accurate (0.7% median offset) than the timedelay set alone. This result confirms that nontimedelay lenses can nonetheless contribute valuable constraining power to the determination of H_{0} via their kinematic constraints, assuming they come from the same global population as the timedelay set.
Key words: gravitational lensing: strong / cosmological parameters / methods: numerical / galaxies: elliptical and lenticular, cD
© M. R. Gomer et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
A precise determination of the Hubble constant H_{0} is a widely sought after goal, now with multiple mature measurement methods. Some methods measure the early universe and recover H_{0} (e.g., Aiola et al. 2020; Planck Collaboration VI 2020; Abbott et al. 2018), while some use standard candles to measure H_{0} locally (e.g., Riess et al. 2022; Freedman et al. 2020). There is a discrepancy between these measurements which may require the addition of new physics into the standard cosmological model (e.g., Poulin et al. 2019), as extensive considerations of possible unaccountedfor systematic effects in one or more methods of H_{0} determination have not provided an explanation yet (e.g., Freedman 2021; Huang et al. 2020; Kenworthy et al. 2019; Addison et al. 2018). Additional independent measurement methods may help diagnose the tension, and should be explored before making claims of new physics. Many such novel methods have been developed, such as using water megamasers (Pesce et al. 2020) or gravitational sirens (Abbott et al. 2017) to directly measure distance. One such distanceladderindependent method is the gravitational lensing timedelay method.
The gravitational lensing method measures distances through the use of time delays between different images of the same variable source. The combination of the observed time delays with the lens mass model gives a timedelay distance:
where z_{d} is the redshift of the lens (deflector), and D_{d}, D_{s}, and D_{ds} are the angular diameter distances from the observer to the lens, observer to the source, and deflector to the source, respectively. Distance (i.e., H_{0}) measurements can be combined between numerous systems to increase the statistical precision. The method requires highquality imaging data, longterm timedelay monitoring, and precise spectroscopy for kinematic measurements. Large collaborations leading the efforts with this method include COSMOGRAIL (Courbin et al. 2018; Bonvin et al. 2019; Millon et al. 2020a), SHARP (Spingola et al. 2018; Chen et al. 2019), STRIDES (Treu et al. 2018; Anguita et al. 2018; Shajib et al. 2019, 2020), and H0LiCOW (Suyu et al. 2017; Wong et al. 2020; Rusu et al. 2020), which together form the TDCOSMO umbrella collaboration (Millon et al. 2020b).
The primary challenges for the lensing method are degeneracies, whereby the same lensing observables (e.g., image positions and relative fluxes) are recovered from a set of mass distributions. Lensing information alone has no ability to constrain which distribution within the set is the truth. Of particular interest is the mass sheet degeneracy (MSD, Falco et al. 1985), a specific case of the more general source position transformation (SPT, Schneider & Sluse 2013, 2014; Unruh et al. 2017; Wertz et al. 2018). The MSD is extensively studied and welldocumented in the literature (e.g., Saha & Williams 2006; Xu et al. 2016; Birrer et al. 2016; Sonnenfeld 2018; Birrer et al. 2020; Kochanek 2020, 2021; Gomer & Williams 2020; Blum et al. 2020; Yıldırım et al. 2021; Birrer 2021). The MSD is described through the mass sheet transformation (MST), in which a mass distribution κ(θ) can be scaled by a factor of λ with the addition of a corresponding constant, analogous with the addition of a uniform sheet of mass, into any one realization from a set of transformed distributions:
where we have used the λ subscript to indicate the quantity has been transformed by an MST. The physical nature of a mass sheet is often conceptually separated into two different phenomena: an “external” mass sheet coming from additional mass along the line of sight, and an “internal” mass sheet which encodes the possibility of an alternative lens profile shape. This work concerns the internal mass sheet. Together with a corresponding transformation of the unobserved source position β_{λ} = λβ, the MST leaves image positions and relative fluxes unaffected, while leading to a multiplicative constant on the timedelay distance measurement D_{Δtλ} = λ^{−1}D_{Δt}. Since the recovered value of H_{0} is inversely proportional to the timedelay distance, a determination of H_{0} from a masssheettransformed model results in a biased value of H_{0λ} = λH_{0}. As such, the capacity of lens modeling to recover an accurate measure of H_{0} is inherently a problem of constraining λ.
Since lensing information alone is unable to constrain λ, one must turn to external information. The impact of the MST can be mitigated via stellar velocity dispersion measurements, which provide an independent measurement of mass (Treu & Koopmans 2002a,b; Koopmans et al. 2003; Barnabe & Koopmans 2007). This addition requires that a given mass distribution reproduces both the lensing data and kinematic data, which places constraints on the profile shape, and therefore the range of allowable λ.
To further constrain λ, it is common to incorporate existing knowledge about galaxy structure through the implementation of a parametric profile. Earlytype galaxy mass profiles are known to be wellapproximated as a power law with respect to radius with approximately isothermal slope, a phenomenon sometimes referred to as the bulgehalo conspiracy (Gavazzi et al. 2007; Suyu et al. 2009; Koopmans et al. 2009; Humphrey & Buote 2010; Auger et al. 2010; Barnabè et al. 2011; Dutton & Treu 2014). By specifically using a power law to fit a lens, a modeler fixes λ by restricting the fit to the κ_{λ} which matches the powerlaw profile form over the relevant radii. As such, the choice of a specific family of model, powerlaw or otherwise, restricts the range of possible MSTs which fit the data. If a modeler uses a power law to fit a lens which in reality differs from a power law, the MST which causes the profile to best match a power law is implicitly applied, which can result in a biased value of H_{0}. One approach to account for this effect is to use multiple different families of models for the same lens, which is why H0LiCOW uses both a powerlaw and a twocomponent composite model to fit timedelay lenses for cosmography (Suyu et al. 2014; Rusu et al. 2020). TDCOSMO’s exploration of systematics found agreement in H_{0} when using a power law or composite model in the fit, providing supporting evidence that this MST effect is not severe in real lenses (Millon et al. 2020b). When using a powerlaw model for 6 lenses supplemented with aperture kinematic measurements, H0LiCOW found km s^{−1} Mpc^{−1} (Wong et al. 2020).
Birrer et al. (2020, hereafter TDC4) introduced a technique to supplement the parametric method by explicitly folding the effect of the MSD into a hierarchical analysis of the lens systems. This method is able to combine sets of lens modeling results together, provided that the sets of lens galaxies come from the same overall galaxy population. The method applies an MST to the lens modeling results on a population level, serving to capture systemic departure from the initial parametric model. To account for the difference of lensing impact parameter arising from having different lens and source redshifts, the amplitude of the MST is allowed to vary as a function of the ratio between the effective radius and Einstein radius. By directly parameterizing and including this effect, this approach allows the velocity dispersion data of the collective set to constrain the space of possible MSTs. Specifically, systems from the Sloan Lens ACS Survey (SLACS, Bolton et al. 2006, 2008; Shajib et al. 2021), which have lens models and velocity dispersions, but lack time delays, are able to be included in the analysis to place constraints on the populationscale MST, given the hypothesis holds that this population comes from the same global population as that of the timedelay lenses. We detail this method in the following section.
The hierarchical framework has been validated for a population of timedelay lenses by TDC4 based on the simulated lenses from Rung 3 of the time delay lens modeling challenge (TDLMC; Ding et al. 2018, 2021a). The value of H_{0} was inferred by modeling hydrodynamically simulated lensed systems where the lensing galaxy was drawn from a set of 16 massive elliptical galaxies from IllustrisTNG and from zoomin simulations. The resulting value of H_{0} ( km s^{−1} Mpc^{−1} when using a uniform prior for anisotropy radius) was consistent with the fiducial input value (65.413 km s^{−1} Mpc^{−1}) with most of the uncertainty coming from limited constraints on the anisotropy of the velocity dispersion. However, the TDLMC Rung 3 lenses are known to have numerical effects in the central region, for which the authors caution against their strict interpretation as realistic lenses (Ding et al. 2021a). We also note that the TDLMC comparison does not test the inclusion of a corresponding nontimedelay lens population. With this in mind, the efficacy of the hierarchical method for real systems is supported, but not confirmed. A test of the method using analytically created mock systems which do not have such artificial cores would supplement this result, which is one goal of this work.
Combining the seven TDCOSMO systems, TDC4 recovered a value of km s^{−1} Mpc^{−1}, which includes the MSD in the error budget. This is consistent with the Wong et al. (2020) value, but with widened uncertainty. Seeking better constraints, 33 nontimedelay lenses from SLACS (Shajib et al. 2021) were implemented into the hierarchical framework, under the assumption that the SLACS and TDCOSMO lenses are part of the same population. Interestingly, the resulting value of H_{0} moves downward to km s^{−1} Mpc^{−1}. While these two values are statistically consistent, the hierarchical method nevertheless relies on different assumptions than other methods and so it must be further tested, motivating the experiment in this paper.
The main concern is that if there were some systematic difference between the SLACS and TDCOSMO lens populations, the primary assumption of the hierarchical framework would be flawed, and the result would be unreliable. Unfortunately, we do not have access to the true mass distributions of these galaxies, and as such we cannot confirm or reject the possibility that both samples are drawn from the same population. We can, however, perform targeted tests of systematics; a natural place to begin comes from the fact that the SLACS and TDCOSMO lenses have different redshifts, and therefore probe different radii relative to the effective radius. Ergo, the second goal of this paper is to test the capacity of the hierarchical framework to correctly capture and account for this change in local lensing effect which stems from a difference in redshift. This goal can be addressed by by creating mock lenses consistent with the SLACS and TDCOSMO lens populations and combining them together hierarchically.
In summary, in this work we further test the hierarchical framework by modeling a large sample of mock lensed systems meant to represent those of the TDCOSMO galaxies, constructed as analytical twocomponent light+dark matter mass distributions. Specifically, we test the capacity of the framework in two respects. Firstly, we wish to confirm that the framework still performs well with intrinsic mass distributions which are less cored in their center than galaxies from hydrosimulations. Secondly, we wish to confirm that the inclusion of the comparison sample of nontimedelay lens systems is correctly accounted for without introducing a bias arising from the change in redshift (i.e., lensing impact parameter). To address this second respect, we construct our lens population such that it is consistent with the TDCOSMO population at one source redshift, and consistent with the SLACS population when the source redshift is changed.
We first describe the hierarchical framework and experiment setup in more detail in Sect. 2, then outline our mock creation in Sect. 3, present our results in Sect. 4 and discuss these results in Sect. 5.
2. Hierarchical framework and experiment setup
The hierarchical framework is detailed in TDC4, but we summarize the main components here. The key idea of the hierarchical framework is that information on the amplitude of the MST can be retrieved by considering a population of lens systems with kinematic data. While singleaperture kinematics can place limited constraints on λ for individual systems, the framework is able to leverage the additional statistical constraining power of a population of systems.
To gain a physical intuition of the method, consider that when an individual lens is fit as a power law (PL), the actual constraint is that the lens profile is a PL+MST with an unknown λ. However, these galaxies as a population are thought to follow approximately the same global mass profile (with a characteristic radius that scales with the luminous effective radius, θ_{eff}^{1}), meaning the λ necessary to match a given lens to a PL is purely a function of which part of the profile the Einstein radius (θ_{Ein}) lands upon. Kinematics also probe some radial range of this global profile, depending on the size of the measurement aperture relative to θ_{eff}, but unlike lensing information, kinematics are sensitive to λ, as they measure the true mass within a spherical shell. As such, the combination of many systems together can be thought of as probing the global profile at a range of radii with both lensing and kinematic data, where the measurements from one system can help support the unconstrained parts of another system. This combined information allows one to infer the populationscale λ behavior associated with the mapping of the global profile into a PL.
One can show that the MSTinvariant term extracted from a lensing measurement is:
where is the circularly averaged radial second derivative of the deflection evaluated at θ_{Ein}, and κ_{E} is the mean convergence of an annulus at θ_{Ein} (Sonnenfeld 2018; Kochanek 2020). From this point of view, the hierarchical inference takes a set of ξ constraints at different Einstein radii and uses the kinematic information to provide the conversion from the MSTagnostic description into the true mass distribution with an associated λ, which is allowed to be different for different Einstein radii.
The inference works by ascribing a distribution of MSTs to the galaxy population, where the λ_{int} distribution is described with Gaussian standard deviation σ(λ_{int}). The “int” denotation references the conceptual separation between an internal mass sheet transformation as described here from an external mass sheet, owing to convergence along the line of sight κ_{ext}. To account for differences in the radii at which the global mass profile is probed, a linear scaling with the effective radius relative to the Einstein radius is included as
with slope α_{λ} and intercept λ_{int, 0}. We note the distinction between the set of λ_{int} values each lens within a population individually undergoes and the intercept of our linear parameterization of the population, λ_{int, 0}, a descriptive parameter we recover. This scaling allows for lenses which are probed at different radii to have different MST amplitudes to match the powerlaw shape. The constraint on λ_{int} comes from kinematic information, for which even lenses without time delays can provide meaningful constraints. Lens population kinematics are included in a spherical Jeans framework parameterizing the distribution of anisotropy a_{ani} (further detailed in Sect. 3.4) as Gaussian with mean ⟨a_{ani}⟩ and standard deviation σ(a_{ani}). In total, the hierarchical framework recovers values for λ_{int, 0}, σ(λ_{int}), α_{λ}, ⟨a_{ani}⟩ and σ(a_{ani}) to describe the distribution of values within the population and place an informed constraint on H_{0}.
The hierarchical inference has been successfully tested on numerically simulated TDLMC profiles, which have artificial cores arising from numerical effects. However, it has not yet been confirmed to be effective for analytical profiles which do not have such cores. In addition to testing this, we sought to test the capacity of the method to incorporate an external nontimedelay lens population which is similar to the timedelay lens population, except with a different relative impact parameter, (i.e., different θ_{eff}/θ_{Ein}), as was done for real systems in TDC4.
To control for all other variables, we chose for the two lens populations to be identical, except for the source redshift, for which we selected one value for each population. Given the same lens surface mass density Σ(θ), a larger source redshift z_{s} decreases the critical lensing density Σ_{crit} and hence increase the convergence κ(θ):
This increased κ leads to a larger Einstein radius, changing which part of the profile is probed by lensing. If a profile has a changing shape with radius, probing different radii could change the resulting fit, resulting in a different MST and ultimately in a change in H_{0}. As such, the expected change in λ_{int} with respect to impact parameter encodes two effects: the change in redshift leading to a different probed radius of a given profile, and the variety of profile shapes arising from the set of twocomponent mass distributions at a given redshift. The hierarchical setup assumes that these effects collectively are welldescribed as linear with Gaussian scatter (Eq. (4)). If the combined effect is not welldescribed by this parameterization, there will be nonrandom departures from the linear description, resulting in an incorrect scatter estimation and a biased H_{0}. However if the linear parameterization is a good description of the combined effect, such that the variations between systems can be captured by the Gaussian scatter, the effect is accounted for and the correct H_{0} will be recovered. As such, an accurate H_{0} recovery in this experiment represents supporting evidence for the versatility of the hierarchical method. On the other hand, if the method were to fail in the case where the only difference between populations comes from source redshifts, the hierarchical framework in its current form would be unlikely to be reliable for the more complicated real population.
This experiment is not necessarily intended to perfectly emulate the data sets for either SLACS or TDCOSMO specifically. Instead, the intention is to perform a targeted test of the effect of probed radius within the hierarchical framework, using a set of lenses engineered to come from the same population.
3. Creating mocks
We constructed a set of twocomponent profiles consistent with the observed lens population. For each profile, we simulated two mock HSTquality images corresponding to the lens with two different source redshifts. Each lens image was modeled using a powerlaw profile for the lensing galaxy. Mock kinematic measurements were calculated and combined with the lensing posteriors in a hierarchical framework. This section details how these mass profiles were created (Sect. 3.1), how the population of these profiles compares to the observed lenses (Sect. 3.2), how mock images were created and fit (Sect. 3.3), and how kinematics were simulated (Sect. 3.4).
3.1. Mass distribution
We constructed lenses using a twocomponent composite mass model consisting of a Chameleon profile (Dutton et al. 2011; Suyu et al. 2014) representing the light distribution and an NFW (Navarro et al. 1996, 1997) component representing the dark matter. We assumed that the mass distribution of the baryon component corresponds to the light distribution of observed systems. Most observed light profiles of earlytype galaxies are modeled as Sérsic profiles, for which the corresponding surface mass density is given as
with axis ratio q where θ_{eff} is the effective radius enclosing half of the total light and n is the Sérsic index (Sérsic 1963; Ciotti & Bertin 1999). The constant b_{n} is set to ensure the halflight property of θ_{eff} given n. Unfortunately the elliptical Sérsic profile requires numerical integrals in the lens deflection calculation and as such is computationally expensive in this context. We therefore chose to represent the baryonic matter as a Chameleon profile (Dutton et al. 2011), which is a profile designed to closely match a Sérsic profile at lensing radii using a combination of two softened isothermal profiles (Kassiola & Kovner 1993),
with parameters w_{c} and w_{t} with w_{t} > w_{c}. The Chameleon profile is often used for this purpose in the context of lensing (e.g., Suyu et al. 2014; Birrer et al. 2019; Rusu et al. 2020; Shajib et al. 2020; Van de Vyvere et al. 2022a).
Given a Sérsic index and θ_{eff}, one can use the prescription of Dutton et al. (2011) to construct an analogous Chameleon profile. We note that this prescription claims that the Chameleon profile matches the enclosed mass of the Sérsic profile to within a few percent within the [0.5θ_{Ein}, 3θ_{Ein}] range, but actually requires a modification to reach this reported precision, which we have implemented and discussed in Appendix A. With this, we were able to construct a Chameleon profile which closely matches a provided Sérsic.
The NFW profile is described in 3D mass density as
with scale radius r_{s} and scale density ρ_{0}. Deflection angles and lens potentials for an elliptical NFW in projection can be calculated. Following Golse & Kneib (2002), ellipticity is added in the lens potential. The parameters for these Chameleon+NFW profiles must be selected to result in a population which is consistent with the observed SLACS and TDCOSMO lensing galaxies (Dutton & Treu 2014; Shajib et al. 2021).
3.2. Creating a lens population
To build up a realistic population of lens systems, we examined the existing population of lens galaxies which have been modeled as a composite profile. We considered the subset of 14 SLACS lenses used by TDC4 which have been modeled by Shajib et al. (2021) and therefore have fit constraints on their slope, γ. These lenses also have stellar masses measured via stellar population analysis using a Salpeter IMF (Auger et al. 2009). Combining these measurements, we examined the distribution of several observed properties (i.e., axis ratio of the light q_{L}, θ_{eff}, θ_{Ein}) as well as some recovered properties: axis ratio of the mass (q_{m}), mass within θ_{Ein} (M_{Ein}), total mass in baryons (M_{bar}), baryon fraction within θ_{Ein} (f_{bar}), and the MSTinvariant ξ = 2(γ − 2) for a power law. We chose to plot ξ as recovered from the lensing powerlaw fits rather than γ directly to better reflect that the fit slope is a property recovered from lensing and may not correspond to the true local slope at the Einstein radius (Kochanek 2020). Additionally, we compiled the same properties for the TDCOSMO set of 7 lenses, although we lack measurements for q, M_{bar} or f_{bar}, and plotted the two populations as histograms in Fig. 1.
Fig. 1. Histograms depict the distribution of observed SLACS (purple) and TDCOSMO (orange) lens properties. Our simulated populations (shown as KDEs of corresponding colors) were created by constructing a large number of profiles and selecting a subset which simultaneously match the 5 target parameter distributions, indicated by the green borders. The simulated SLACS and TDCOSMO populations are the same profiles, only differing by a changed source redshift. For parameters where the two sets are identical, a black line is used instead of orange or purple. 
We sought to create a synthetic population of twocomponent profiles which closely matches both observed populations. To create a twocomponent profile, 7 parameters (omitting centroid positions and position angles) had to be given values. For the baryon component we have Σ_{eff}, θ_{eff}, n, and q_{L}. For the dark matter halo component we have ρ_{0}, r_{s} and q_{h}. Furthermore, we also must select redshifts for the source of each lens.
The first simplifying assumption we made is that the two components have the same axis ratios (q = q_{L} = q_{h}). We also set the position angles and centroids to be aligned. In reality, these components need not be aligned (Gomer & Williams 2018; Shajib et al. 2019), and have been allowed to vary in lensing models (e.g., Rusu et al. 2020; Millon et al. 2020b). However, since the primary focus for this work is the effect of the change in Einstein radius regarding the radial profile, rather than the effects of the angular structure of the lens, we fixed the components to be aligned. This assumption is supported by the joint lensing+dynamics analysis of Shajib et al. (2021), who find that centroids of light and mass match within 218 pc (68% confidence) and position angles align within 10 deg for most SLACS galaxies, consistent with earlier results for quasar lenses (Keeton et al. 1998; Yoo et al. 2006; Sluse et al. 2012). Numerical results using the Illustris simulation also support this level of alignment (Xu et al. 2017). We therefore chose to adopt this simple angular structure for the input mocks, although the lens model we use to fit the systems does have additional azimuthal freedom through external shear.
The intention of our mock population is that it resemble the SLACS distributions for a particular source redshift, z_{s}, and the TDCOSMO distributions after a change to another z_{s}. In reality, both populations exist over a range of z_{s} and lens redshift z_{d}, shown in Fig. 2. Generally, the TDCOSMO set has larger z_{s} and z_{d}. For the targeted experiment in this work, we selected just one z_{d} to be the same for both sets to set aside the effects of a changing angular scale (which affects image resolution and kinematic aperture size) and a changing cosmological scale factor (which changes comoving values of effective radii). We selected z_{d} = 0.25, which is higher than the typical SLACS value and lower than the typical TDCOSMO value. Strictly speaking, this lens redshift is below the lowest TDCOSMO z_{d}, but serves as a nice middle ground which is capable of retaining a Σ_{crit} for each set which is consistent with each respective population. We selected just two z_{s}: one to represent each of our two populations. We selected a z_{s} for what we call the “SLACSlike” set to be 0.6 and a z_{s} for what we call the “TDCOSMOlike” set to be 2.0. These redshift selections result in Σ_{crit} values which are close to the medians of those of the respective SLACS and TDCOSMO populations. These calculations use a flat universe with H_{0} = 70 km s^{−1} Mpc^{−1} and Ω_{M} = 0.3.
Fig. 2. Redshift comparison with the observed lens population. Left: the distributions of source redshifts against lens redshifts for the populations of observed SLACS (purple) and TDCOSMO (orange) lenses. Our selections of source and lens redshifts are indicated as dashed lines. Right: these selections result in the dashedline Σ_{crit} values, which are consistent with the observed populations. 
For the baryonic mass component, we set the Sérsic index of the light component to be n = 4, equivalent to a de Vaucouleurs profile (Sérsic 1963). We then used the prescription described in Appendix A to create a Chameleon profile from an input Σ_{eff} and θ_{eff}.
Even once a model family is chosen, there still is considerable freedom in the handful of parameters required to make a lens. Some parameters (such as q) can be directly observed and inserted into a model, while others (such as NFW scale radius r_{s}) do not easily convert into the corresponding observable constraint (like θ_{Ein}, which indirectly is determined by r_{s}, in conjunction with other parameters). We wish to correctly probe the space of observed values, but must create our mocks using parameter values which are sometimes not observable. To solve this problem, we randomly sampled the input parameters over a range, resulting in a large set of profiles with wouldbe observed values. Comparison to the observed distribution provides a set of target parameters, from which one can select only the realizations which match the targets to our desired precision. This was done by defining a χ^{2} for each targeted quantity and selecting a subset of profiles below a threshold total χ^{2}. We selected a population which targets the observable quantities q, θ_{Ein}, and θ_{eff}^{2}^{3}.
After selecting a number of profiles which met our targets, we arrived at a population of 465 lens profiles, which are a good match to the observed lens populations^{4}, with distributions of parameter values approximated via Kernel Density Estimation (KDE) in Fig. 1. Aside from the intended match to the target parameters, we were also able to compare some nonobservable quantities which were not targeted directly, such as baryon fraction or total mass. These quantities also match well between the observed and simulated populations, with the possible exception that the recovered values of ξ were closer to isothermal than some of the observed SLACS lenses, which recover steeper slopes (recall that ξ = 0 corresponds to an isothermal slope). Nevertheless, our simulated distribution overlaps significantly with the real lenses. In particular we have reproduced the difference in impact parameter between SLACS and TDCOSMO lenses: Fig. 3 shows that TDCOSMOlike lenses have images forming at a larger distance (smaller θ_{eff}/θ_{Ein}) than SLACSlike systems. This agreement with the observed data is important for testing the relevance of the parameterization of the MST with radius (Eq. (4)) used in the hierarchical framework.
Fig. 3. Scaled radial quantity used by the hierarchical framework to parameterize radial dependence for MST, plotted for the observed population (histograms) and the simulated profiles used in this work (KDE curves). 
We plot the convergence profiles as a function of radius for our mock populations in Fig. 4. In comparison with the hydrodynamically simulated lenses from the TDLMC previously used to validate the hierarchical framework, our analytical profiles are less cored, with slight differences in shape at larger radii.
Fig. 4. Convergence profiles for 20 analytical Chameleon+NFW profiles used in this paper (purple: SLACSlike; orange: TDCOSMOlike) alongside 16 profiles from hydrosimulations used in the TDLMC (green). Colored bands correspond to the 16th and 84th percentiles. Vertical colored lines indicate the median Einstein radius of the corresponding population. The diagonal black line indicates an isothermal slope for comparison. 
We now have a population of mock lens profiles which are comparable to both observed populations, except that the only difference between our SLACSlike and TDCOSMOlike populations is that we have placed the sources at two specific source redshifts. By drawing a random subset from this population, we created a set of lenses to be modeled while still maintaining confidence that our set is realistic.
3.3. Creating mock images and fits
We followed the setup of Van de Vyvere et al. (2022a) to create HSTlike mock images using lenstronomy^{5} (Birrer & Amara 2018; Birrer et al. 2021). The images were created on a 151 × 151 pixel grid with a resolution of 0.08 arcsec and were convolved by the same PSF as used in the TDLMC (Ding et al. 2021a). We included noise using a zeropoint of 26 mag, with an exposure corresponding to 5400 s and sky brightness of 22.3 mag arcsec^{−2} with a read out noise of 21 e^{−}, characteristic of the HST WFC3 camera in the F160W filter (Dressel 2012).
For each profile, two images were created: one SLACSlike image with z_{s} = 0.6 and one TDCOSMOlike image with z_{s} = 2.0. In all cases, we restricted our analysis to quadruply imaged configurations. The source was placed randomly within the caustic area in the same location for both source redshifts, which required rescaling the source position by the same factor as the increased caustic size for the higherz_{s} TDCOSMOlike set. The source light was a circular Sérsic profile with n = 3, θ_{eff} = 0.1″ and a total unlensed magnitude of 20.5. Lacking an AGN in reality, the SLACSlike set was not given a point source or time delays. Meanwhile for the TDCOSMOlike set we included a central point source representing an AGN with an unlensed magnitude of 21. These magnitudes were selected to match the observed contrast between the quasar and host galaxy (Van de Vyvere et al. 2022a) and correspond to a system approximately 0.5 mag brighter than the mean intrinsic magnitude of the TDCOSMO lenses measured by Ding et al. (2021b). Time delays were included, centered on their true values, with an uncertainty of 2% or 1 day (whichever is larger) for likelihood calculation. For both sets of lenses, light from the lens galaxy was not included in an effort to focus the exploration on the changing Einstein radius rather than any possible effects of imperfect lens light subtraction. The end result for each lens was two mock highquality HST images with different Einstein radii probing different parts of the same lens profile. Ten example lenses are shown in Fig. 5.
Fig. 5. Images (log brightness) and residuals (in units of σ) of 10 systems at 2 source redshifts. Left: SLACS analog with z_{s} = 0.6 and no point sources, Right: TDCOSMO analog with z_{s} = 2.0. Lowerleft scale bars correspond to 1″. 
The images were fit via Particle Swarm Optimization (PSO, Kennedy & Eberhart 1995; Shi & Eberhart 1998) using lenstronomy with a Sérsic model for the source light and a Powerlaw Elliptical Mass Distribution (PEMD) + external shear model for the mass. Because external shear was included in the fit model but absent in the mock creation, the fit model has more azimuthal freedom than the input. This reduces the possibility for a bias to arise due to insufficient azimuthal freedom (Kochanek 2021). The central 0.4″ were masked in the fitting process to stand in for the sometimesuncertain nature of lens light subtraction in the centermost region, remaining agnostic to the possiblity of a fifth image (although our PEMD model does not produce a fifth image). We also masked all pixels with intensity below twice the background level, as in Van de Vyvere et al. (2022a). Because this work uses mocks where the source properties are known, we were able to initialize a preliminary fitting using the correct source, but we did not fix the source to the truth during the fitting procedure. We fixed the slope to isothermal for this preliminary fit. A secondary fit started from this result and allows all parameters to vary. From the second result, a Markov chain Monte Carlo (MCMC) via emcee (Goodman & Weare 2010; ForemanMackey et al. 2013) was used to sample the posterior distribution of the lens parameters and estimate uncertainties.
All systems were successfully fit with the PEMD+shear models, with reduced χ^{2} < 1.06 in all cases. Before accounting for the MSD, these fits recovered H_{0} with a median of 58.6 km s^{−1} Mpc^{−1}, which is biased by approximately 16% from the input value, with approximately 3% scatter. An example comparison between the input and fit profiles for a single pair of systems is given in Appendix B. The next steps were to add kinematic information to break the MSD and combine the posteriors together hierarchically.
3.4. Kinematics implementation
Stellar kinematic information is used in conjunction with lensing information to break the MSD. To do the same in this work, we had to create for each lens a predicted velocity dispersion corresponding to a mock aperture measurement.
The hierarchical framework uses spherical Jeans modeling as detailed in Binney & Mamon (1982), Birrer et al. (2020). This approximation is commonly employed in lens modeling (e.g., Treu & Koopmans 2002a; Koopmans et al. 2003; Suyu et al. 2010; Sonnenfeld et al. 2012; Wong et al. 2017; Birrer et al. 2019; Rusu et al. 2020) and calculates the velocity dispersion in the limit of a relaxed, spherically symmetric, nonrotating system. Radial () and tangential dispersions () are described by the spherical Jeans equation:
for a given stellar distribution ρ_{*}(r) and gravitational potential Φ(r) where the anisotropy β_{ani}(r) is defined as
The luminosityweighted projected velocity dispersion is given as
where R is the 2D projected radius and Σ_{*}(R) is the 2D projected stellar density. Finally, to produce an observable measure, one needs to integrate this quantity over an aperture 𝒜 and convolve by a PSF 𝒫:
All told, one can calculate the observed aperture velocity dispersion given a mass distribution, a light distribution, and a description of the orbital anisotropy. One common parameterization of the anisotropy is to use
which describes the orbits as isotropic at small radii and radial at large radii, with transition radius r_{ani} (Osipkov 1979; Merritt 1985). It is also useful for our purposes to describe this transition radius relative to the effective radius using the anisotropy parameter
The spherical Jeans framework has the advantage that it is tractable to calculate quickly for a large number of models, although this is a simplification of the true dynamics of galaxies. We initially experimented with using the more general Jeans anisotropic modeling (JAM; Cappellari 2008, 2020; Shajib 2019), which is starting to be explored for lens systems (e.g. using mock JWST kinematics, Yıldırım et al. 2020, 2021), to calculate an input velocity dispersion while keeping the hierarchical modeling framework as spherical Jeans. However, there are subtleties regarding the differences between the two methods which make interpretation difficult. The main challenge is that the anisotropy is parameterized differently in JAM (, instead with respect to the zcoordinate of a cylindrical coordinate system), which can change the aperturemeasured velocity dispersion. Therefore, the experiment would not be a test of the hierarchical framework’s multipleredshift capabilities, but instead a test of the accuracy of anisotropy assumptions. Determination of an optimal description of anisotropy for lensing work proves to be a significantly more involved task than the scope of this work, and requires a focused study of its own. In the future, it would be preferable to upgrade the hierarchical framework to work with JAM, but the challenge is that JAM is too computationally expensive to sample within an MCMC. Another project is currently underway to attempt to speed up this calculation through the use of a neural network to precalculate the posterior space for MCMC sampling (paper in prep.). In this present work, we simply used spherical Jeans modeling in both the input and in the hierarchical fitting.
Regarding previous work using Jeans modeling within the hierarchical framework, we note that the kinematics calculation of Birrer et al. (2020), Birrer & Treu (2021) approximated the light distribution as a Hernquist profile for the ρ_{*} tracer population in Eq. (9). Since our light distribution is a Chameleon profile, we updated lenstronomy so that it is capable of calculating kinematics for Chameleon profiles (see Appendix A for details). This way, we avoid any potential bias that could come from a mismatch between the true light distribution and the one used to calculate the kinematics.
The specific kinematic setup used in this work is the same throughout. For each lens we calculated mock velocity dispersions within 3 apertures with radii of 0.67″, 1.33″, and 2″, convolved with a seeing of 0.7″. Each dispersion was given an uncertainty of 10%. For real systems, IFU spectroscopy is capable of measuring 2D kinematic constraints (such as in MaNGA e.g., Bevacqua et al. 2022) which correspond to multiple aperture constraints (Shajib et al. 2020). Some of the existing lens set already has such constraints via MUSE, typically made with ≃5% precision (see Table E.1 of TDC4 for the SLACS set). Furthermore, resolved kinematics should be possible with JWST (Yıldırım et al. 2020) and adaptiveoptics IFU data such as MUSE (Bacon et al. 2010), KCWI (Morrissey et al. 2018) and OSIRIS (Larkin et al. 2006), providing even higher precision constraints for some subset of the lens sample. Indeed, IFU kinematics are anticipated within the TDCOSMO forecast using the hierarchical framework (Birrer & Treu 2021). As such, our kinematic setup uses a conservative estimate of the constraining power of real observations. For our mock kinematics, we must set an input anisotropy a_{ani}. For real systems a_{ani} must take on a range of values but the exact physical distribution is not wellknown. To focus the experiment on the effect of changing impact parameter rather than the parameterization of anisotropy, we set input a_{ani} = 1 for all systems. The limitations of this assumption are further discussed in Sect. 5.3.1. This describes the mock measured velocity dispersions calculated for the true profile.
We must also calculate the predicted velocity dispersions from the lensing model. This was done using the same spherical Jeans modeling for the PEMD fit profile, resulting in a measured velocity dispersion for each aperture. The posteriors of the lensing fit were sampled, with uncertainties in the fit slope and Einstein radius included. We also introduced a 10% uncertainty to the effective radius in this sampling. Considering the effect of the MST on velocity dispersion, we note that the infinite uniform sheet term in Eq. (2) does not affect kinematics. As such, the effect of the MST on the velocity dispersion is a modification by a constant factor of which follows from the application of the coefficient λ_{int} to the density in Eq. (11) (or from the simple scaling that v^{2} ∝ GM/R). The hierarchical framework seeks to match the predicted and measured velocity dispersions, sampling over the parameters which describe the populationlevel a_{ani} and λ_{int}.
A summary table of the parameters used for mock lens creation is made available. See Appendix C.
4. Results
We now present results for an ensemble of lens systems. We used the publicly available hierArc package^{6} (Birrer et al. 2020) to sample over the posteriors from the PEMD fits to each system. From the lensing posteriors and kinematic information, hierArc samples for H_{0}, using a linear parameterization of the populationscale MST (Eq. (4)) with intercept λ_{int, 0} and slope α_{λ}, with a spread of σ(λ_{int}). Additionally, the distribution of anisotropy is described as Gaussian with median ⟨a_{ani}⟩ and standard deviation σ(a_{ani}). This work uses uniform priors for all of these parameters with bounds listed in Appendix C. One can compare the results from including just the lenses from the one or the other source redshift, or a combination using the populations from both source redshifts.
4.1. 20 systems
The results for 20 Chameleon+NFW lenses are plotted in Fig. 6. A comparison of the SLACSonly population, the TDCOSMOonly population, and the combined populations indicates that the combination can increase the precision and accuracy on the final value of H_{0}, even though the systems span different ranges of θ_{eff}/θ_{Ein}.
Fig. 6. Results of the hierarchical framework using 20 Chameleon+NFW lenses for the population with SLACS source redshift and no time delays (purple), the population with TDCOSMO source redshift (orange) and the combination of both populations (blue). 
The TDCOSMOonly set recovered a median H_{0} which is higher than the fiducial value of H_{0} by approximately 6.2% (1.5σ). The broad width of the posterior is largely due to the degeneracy with ⟨a_{ani}⟩, which was poorly constrained from a singleredshift set of lenses.
The SLACSonly population had no point sources, therefore no time delays; as a result, it could not place a constraint on H_{0} alone, as evidenced by the wide and unconstrained purple contours in the first column of Fig. 6. Rather, the contribution of these systems is to inform the parameters other than H_{0}: those concerning λ_{int} and a_{ani}. This information is combined with the TDCOSMO systems (which themselves bring information about the H_{0} posterior) to break the MSD. This combination enabled the fit to recover an unbiased H_{0} value of km s^{−1} Mpc^{−1}. As a particularly nice visual example, the panel which shows α_{λ} and λ_{int, 0} (middle row, second column) illustrates how these parameters, while degenerate within the SLACS posteriors, were nonetheless able to inform the TDCOSMO set to recover an accurate description of λ_{int, 0}. As in TDC4, the combinedpopulation result was shifted downward, but still consistent with the TDCOSMOonly value of H_{0}.
4.2. 40 systems
While 20 systems is a reasonable shortterm goal of TDCOSMO, a larger number of lenses will soon become available with the advent of LSST, Euclid, and larger telescopes. Increased sample size should increase the precision on H_{0}, but intrinsic scatter in λ_{int} may inhibit this precision increase. We took this opportunity to test if an increased number of lens systems within the hierarchical framework can correctly quantify and account for the intrinsic scatter in λ_{int}. As such, we doubled the number of systems to 40 at each redshift.
With 40 systems at each redshift, we recovered the results in Fig. 7. As intended, the combined SLACS+TDCOSMO value errors did indeed shrink, zeroing in on the correct value of H_{0}, with a recovered value of km s^{−1} Mpc^{−1}.
Fig. 7. Effect of the number of lenses: the hierarchical result of a population of Chameleon+NFW lenses, combining the fits from the TDCOSMO and SLACS source redshifts. The blue is the same as that in Fig. 6, using 20 lenses from each redshift. The pink distribution corresponds to the same setup, but with 40 lenses each. 
4.3. Singleredshift test
We were curious about the importance of probing a range of redshifts, that is to say spread over a greater range of impact parameter. To discuss this, we must compare against the case in which the lenses come from the same redshift. In Appendix D, we show the results of a test where we set the systems to the same redshift and ran hierArc. With 20 time delay and 20 nontimedelay lenses all at the SLACS z_{s}, the result was km s^{−1} Mpc^{−1}.
5. Discussion
The goal of this project has been to test the systematics of the hierarchical framework. The framework has already been tested on a single sample of timedelay lens systems from hydrodynamical simulations; we have expanded this analysis to a different input model using twocomponent analytical profiles consistent with the observed lens population. Furthermore, we have expanded the analysis to test the effect of combining the sample with an analogous population of nontimedelay lenses with different source redshifts.
5.1. Prehierarchical result
To better understand the results of the hierarchical fitting, it can be useful to compare to the lensing results alone, before the MST is accounted for. Figure 8 shows the lensingonly recovered value of H_{0} for 20 systems at each z_{s}, expressed as H_{0,true}/H_{0,model} = λ_{int}, as a function of the radial scaling parameterization θ_{eff}/θ_{Ein} − 1. We note that since the SLACSlike systems have no point sources, they do not recover H_{0} values, and so we have instead plotted the results from the test where the same 20 lenses were given point sources in Sect. 4.3. The hierarchical parameters of λ_{int, 0}, σ(λ_{int}), and α_{λ} reflect the summary statistics of a linearfit description of this plot. We can see λ_{int} has a median of approximately 1.18 and scatter of 0.03, which should be approximately the truth values of λ_{int, 0} and σ(λ_{int}), assuming α_{λ} = 0 (which seems reasonable given the weak correlation, Pearson R = −0.18).
Fig. 8. Lensingonly values of H_{0} (i.e., λ_{int}) as a function of radial impact parameter defined in Eq. (4), which assumes this relationship is described as linear with Gaussian scatter. The SLACSz_{s} systems are indicated in red rather than purple because they are alternate versions of the same systems but with point sources and time delays (such that a value for H_{0} can be recovered and plotted). 
The lensing data is supplemented with kinematic data, and so it can be enlightening to see the velocity dispersion values from the model fits to the lensing data. Before any MST correction was applied, and setting a_{ani} to the truth value of 1, we found that all fits from both systems recovered model velocity dispersions (σ^{P}) which were approximately 7.5% lower than the true values. For individual systems this bias value is contained within the uncertainties of individual σ^{P} measurements (10% uncertainties), but at a population level systems require an average λ_{int} > 1 to reproduce the (observed) lens velocity dispersion compared to the powerlaw model. This illustrates the necessity to use information from the population of lenses rather than individually.
5.2. Performance of the hierarchical framework
We first consider the capacity of the hierarchical framework to reproduce H_{0} from a population at a single source redshift and therefore a narrow range of θ_{eff}/θ_{Ein}. Considering the case with 20 systems using only the TDCOSMOlike source redshift (orange set in Fig. 6), the spread in anisotropy radius ⟨a_{ani}⟩ propagates into 4.4% uncertainty in H_{0}, with a true value which lies within the 95% confidence interval of of the posterior. When expanded to 40 systems (not plotted), this resulted in km s^{−1} Mpc^{−1}, again with the true value within the 95% confidence interval, but now with an uncertainty of only 2.8%. It is interesting that the result lies approximately 1.5σ away from the truth in both cases, even as the value of σ narrows. This may indicate that discrepancies still could exist within the assumptions in the framework at the ≃1σ level, but 1.5σ is not severe enough to claim significance. Nevertheless, this precision represents a substantial improvement over the MSTbiased results from lens modeling alone, which had a median bias of 16% below the fiducial value. This illustrates the success of the hierarchical framework in breaking the MSD for a single population.
We further consider the capacity of the hierarchical framework to account for a change in the impact parameter for a set of otherwise similar lenses. The SLACSlike population is identical to the TDCOSMOlike except that they have no point sources and have a smaller z_{s} resulting in a larger θ_{eff}/θ_{Ein}. Here we compare the result from the singleredshift population to that of the combined population. The combinedpopulation constraints on ⟨a_{ani}⟩ are more informative than either individual set. The degeneracy between ⟨a_{ani}⟩ and H_{0}, visible in the TDCOSMOonly set of Fig. 6, shows how this improved anisotropy constraint leads to a more precise H_{0} (2.6% scatter for 20 systems). The accuracy is also improved: our combined set for 20 systems was centered just 0.3σ away from the correct value. When we increased the number of systems to 40, the precision and accuracy improved further, to a final value of km s^{−1} Mpc^{−1} (0.7% median offset, 1.8% scatter). The increase in accuracy comes primarily from tighter constraints on λ_{int, 0} as the large number of systems is better able to probe the systematic behavior. Nonzero σ(λ_{int}) was also recovered, meaning the framework was able to capture the presence of the intrinsic variance in the input population. Based on Fig. 8, this scatter is anticipated to be approximately 3%, which was wellmatched by the recovered value of σ(λ_{int}). Because the hierarchical framework explicitly captures this scatter through σ(λ_{int}), intrinsic scatter in the distribution of λ_{int} does not represent a systematic floor for H_{0}, that is to say it does not preclude the recovery of H_{0} at a better precision than the intrinsic scatter, and so one can expect the precision to improve as . Importantly, the correct distribution of a_{ani} was recovered, which forced H_{0} to the correct value. The power to constrain a_{ani} comes from having multipleaperture kinematic measurements, providing a probe of the behavior of velocity dispersion with radius. In summary, we confirm that the hierarchical framework is capable of incorporating lens populations from different source redshifts.
When considering the test where all lenses are placed at the same redshift, we find that the framework successfully recovered H_{0} in this test as well. We conclude that the method can result in an accurate value of H_{0} when using either a set of lenses at the same redshift or combining sets from different redshifts.
Birrer & Treu (2021) forecasted the improvement in precision achieved within the hierarchical framework by increasing the size of the lens population, assuming improved constraining power of the velocity dispersion from integral field unit (IFU) spectroscopic data. Unlike this work, Birrer & Treu (2021) did not create and fit mocks, but simply assumed Gaussian constraining power for each individual lens. With such a setup, they found that with 40 timedelay lenses and up to 200 nontimedelay lenses, the precision can be expected to increase to 1.2 − 1.5%, depending on the exact setup. Our result using simulated mocks from start to finish found 1.8% precision with 40 timedelay and 40 nontimedelay lenses, consistent with this forecast.
We note that the general behavior of having a lower value of H_{0} for the combined sample than the TDCOSMOonly sample qualitatively resembles the result from TDC4, where the same hierarchical framework was used on the observed lenses. TDC4 found km s^{−1} Mpc^{−1} when using 7 TDCOSMO lenses alone, but found that the value shifted to km s^{−1} Mpc^{−1} when incorporating 33 SLACS lenses. Speaking in terms of the downward uncertainty of the TDCOSMOonly set, the value shifted by 1.1σ, whereas for this work we saw a similar 1.3σ shift. While this work cannot confirm or refute the mutual similarity between the real lens populations, it is reassuring to see consistent results with the real population, as it affirms that the observed shift is not necessarily an indicator of dissimilar populations. Such mutual similarity between the populations is the key assumption required to make use of the SLACS lenses hierarchically; this work finds no evidence against this assumption.
5.3. Limitations
5.3.1. Mock data
The results of this project look optimistic for the hierarchical framework, but we would be remiss to not discuss the differences between this work and TDC4, as well as some of the present limitations of our method. Firstly, some choices made for the mock data quality may differ slightly from the real lens population. This work has a larger number of systems compared to TDC4, which at the time included 7 TDCOSMO lenses and 33 SLACS lenses of which only 9 had multipleaperture velocity dispersion measurements, with approximately 5% precision on σ^{P}. Here we used 20 (or 40) of each set, all of which have 3 kinematic apertures with 10% precision on σ^{P}. Image data quality in this work (e.g., source brightness, noise, exposure time, etc.) was modeled after the high S/N template from Van de Vyvere et al. (2022a), and as such are higherthanaverage quality compared with current singleband data. We also used transparent lenses, which assumes perfect lens light subtraction. On the other hand, we used singleband images while the observed multiband imaging add additional constraining power and can help deblend lens and source light (e.g., Shajib et al. 2020; Schmidt et al. 2022). This work is intended to be forwardlooking, and the sample of real lenses will grow over time to better match these settings in the future.
By design, this work uses only two source redshifts and one lens redshift to more accurately probe what is thought to be the key difference between the SLACS and TDCOSMO lenses. The real lens population spans a range of z_{s} and z_{d}, but since our populations already share similar distributions of θ_{eff}/θ_{Ein}, implementing this change would be unlikely to make a significant difference. However, if there are different selection effects between the observed SLACS and TDCOSMO samples, or changes to the global profile shape associated with galaxy evolution between the two sample lens redshifts, these effects are not captured in this work.
Regarding our choice of profile, the Chameleon+NFW profiles match many aspects of the observed systems, but may not perfectly describe the truth. Our SLACSlike ξ distribution has less scatter and is more isothermal than the observed set. Additionally, the fact that our PL lensing fits resulted in a different H_{0} value than the truth may indicate a different behavior than the real lens population, which find similar values of H_{0} whether a PL or composite model is used (Millon et al. 2020b). However, to call this a discrepancy assumes that fitting our systems with composite models would return the correct H_{0} (no MST), which we do not necessarily claim. At any rate, we can state that the hierarchical framework successfully accounts for the MSD using either our analytical mocks or the hydrosimulated mocks of the TDLMC, supporting the claim that it is reliable for the true population of (unknown) lens profiles.
A salient difference between our mocks and real galaxies is their elliptical shape, while real galaxies harbor more complex azimuthal structures. The level of impact these structures can have on H_{0} inference is debated. Earlytype galaxies can be triaxial (van de Ven et al. 2009; Chang et al. 2013; Weijmans et al. 2014), and can have 2D mass projections which can change in ellipticity or position angle with radius (Kormendy et al. 2009). We initially sought to create lenses from 3D triaxial profiles for this work, but challenges with the kinematics calculation forced us to simplify to 2D elliptical profiles. Furthermore, galaxies have been known to have “disky” or “boxy” isophotes, corresponding to an n = 4 multipole perturbation to their mass distributions (Hao et al. 2006). Finally, in the ΛCDM paradigm there are expected to be a large number of subhalos within any given galactic halo which add substructure to the mass distribution (Gilman et al. 2020). Additional angular complexity of a lens can lead to a bias in H_{0} when a model is not adequately azimuthally complex (Kochanek 2021). The effects of some of these angular structures in the context of lensing measurements of H_{0} are explored in other works (e.g., substructure by Gilman et al. 2020, multipole by Van de Vyvere et al. 2022a, isodensity twists by Van de Vyvere et al. 2022b). However, it is clear that any substructures are a subdominant factor to the elliptical shape of the mass distribution as evidenced by the success of cosmographygrade parametric lens models. As such, our simple elliptical lenses are likely a reasonable approximation for this work, for which the purpose is to probe the effect of a changing radial structure.
Along the same lines, our spherical Jeans framework used to calculate both the input and the fit velocity dispersion is a simplification. There is an inconsistency between using an elliptical mass distribution for the lens and a spherical mass distribution for the velocity. Yıldırım et al. (2020) found, by using JAM in the spherical limit, that this approximation did not bias the resulting D_{Δt} measurement for their single system, although a formal inclusion of λ_{int} is not incorporated. This spherical assumption is common in joint lens+kinematics modeling, in which kinematics are meant only to provide a single integrated measure of mass, but may become more important as spatially resolved kinematics become more widely attainable. Perhaps of greater concern is the description of anisotropy through the single value of a_{ani}. This parameterization is known to be a simplification which does not accurately match Nbody simulations (Mamon & Łokas 2005). Even so, it is unrealistic to expect each galaxy to have the same value of a_{ani} as has been set in this work as well as in Birrer & Treu (2021). Preliminary tests with changing the input distribution of a_{ani} have found that H_{0} is still recovered well when a representative prior is used, but it is possible for biases to arise if the prior is misinformed (similar to the findings of Birrer et al. 2016, 2020). There is also room to alter and explore the kinematics settings (e.g., number of apertures, radii, seeing, uncertainties) to learn if there is an optimal strategy for an ideal kinematics setup. A more thorough quantification of these effects is a target for a future paper. A future goal would be to incorporate the JAM framework into hierArc, which at present may be too computationally demanding.
5.3.2. MST parameterization
Here we describe an ancillary result with possible implications for the theoretical description of the MST. We sought to supplement the empirical success of the linear parameterization of λ_{int} with impact parameter with an analytical description. The intention was to plot a “truth” theoretical relation in Fig. 8 representing a mapping from Chameleon+NFW to PEMD profiles, and determine if a linear parameterization was appropriate. We performed some tests using the ξ values of each profile to predict the λ_{int} one would recover from a PL fit (since the PL fit should recover the same ξ as the true profile). However, we found that the predicted H_{0} using ξ differs from the actual value recovered from the mock fitting by approximately 3% (median). After considerable exploration, we find that approximately 60% of this discrepancy can be explained by effects of ellipticity: ξ is defined using a spherical lens for which the image radius is unambiguous. However, the remaining discrepancy is still not fully explained, and may be related to other degeneracies in the parameters, possibly dependent on image configuration. This discussion will be the subject of a separate publication (Gomer et al., in prep.). For this current work, this simply means that we cannot reliably use this ξ description to represent a theoretical target.
Finally, we mention in passing that the way the MST is parameterized physically in 3D kinematics is open to interpretation. In this work, the velocity dispersion is modified by a constant factor of , noting that the infinite sheet term in Eq. (2) does not affect kinematics. An alternative description by (Blum et al. 2020) describes the MST as a cored profile, which is wellbehaved at infinity and would leave a detectable imprint on kinematics (Yıldırım et al. 2021). Other descriptions are possible, as the deprojection of the MST is not unique. This work is confined to the first description, but future versions of hierArc could implement other parameterizations.
6. Conclusion
In this work, we conducted an experiment by creating a population of pairs of mocks of the same lens profile using two different source redshifts, representing the difference between the SLACS and TDCOSMO lens populations, controlling for any other possible differences. These images were constructed with twocomponent profiles, but fit as a power law, resulting in a ∼15% biased recovery of H_{0}. The posteriors from these fits were then combined together and supplemented with spherical Jeans velocity dispersion measurements using the hierarchical framework of Birrer et al. (2020). Our resulting value of H_{0} using a single population of TDCOSMOlike lenses was consistent with the fiducial H_{0} (using 20 lenses: 4.2% precision, 1.5σ median offset). When information from the SLACSlike population of lenses was used to supplement the recovery, H_{0} was recovered more precisely and more accurately (using 20 lenses from each set: 2.6% precision, 0.3σ median offset). This precision can be further improved as the number of systems grows (using 40 lenses from each set: 1.8% precision, 0.4σ median offset).
Even though the intrinsic profile was more complex than the powerlaw model used to fit it, we confirm that the hierarchical framework was able to combine systems from different redshifts and adjust the parameterization of λ_{int} as a function of θ_{eff}/θ_{Ein} in a way which results in an unbiased H_{0}. While there may be other possible parameterizations, the linear function with respect to scaled effective radius seems acceptable, provided the data quality is good and the number of lenses is large.
There are still a number of open questions for the hierarchical framework, such as the effect of the distribution of a_{ani} and the proper way to parameterize the anisotropy of galaxies, or the effect of different 3D physical interpretations of the MSD. These questions will be further explored in future work, but this work validates the utility of the method and serves as an important first step.
The Chameleon profile drops more quickly at large radii than the Sérsic it is meant to emulate. As such, the mass at infinity of the baryon component may not perfectly match a Sérsic. For θ_{eff} (and R_{eff}), the value of integrated light at infinity of the Sérsic profile was used, rather than the artificially truncated Chameleon. The other quantities which are not defined using behavior at infinity are expected to be sufficiently accurate: the Einstein radius changes by about 0.01″ comparing a Sérsic+NFW to a Chameleon+NFW. As such, for calculation of the quantities other than θ_{eff} (and R_{eff}), the Chameleon profile was used directly.
The only minor difference arises because of the z_{d} = 0.25 compromise, which resulted in the physical scale conversion being different for the TDCOSMO lenses than for our lenses, and so even though the Einstein radius and effective radius match well in terms of kpc, our synthetic values are slightly larger than the median TDCOSMO values in terms of arcseconds.
Acknowledgments
The authors thank Sherry Suyu, Akın Yıldırım, Martin Millon, Alessandro Sonnenfeld, and Matt Auger for their constructive discussion regarding this paper. This work uses the following Python packages: Python (Oliphant 2007; Millman & Aivazis 2011), Astropy (Astropy Collaboration 2013, 2018), Numpy (van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Pandas (McKinney 2010; The pandas development team 2020), and Seaborn (Waskom 2021). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 787886).
References
 Abbott, B. P., Abbott, R., & Abbott, T. D. 2017, Nature, 551, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, T. M. C., Abdalla, F. B., Annis, J., et al. 2018, MNRAS, 480, 3879 [Google Scholar]
 Addison, G. E., Watts, D. J., Bennett, C. L., et al. 2018, ApJ, 853, 119 [Google Scholar]
 Aiola, S., Calabrese, E., Maurin, L., et al. 2020, JCAP, 2020, 047 [Google Scholar]
 Anguita, T., Schechter, P. L., Kuropatkin, N., et al. 2018, MNRAS, 480, 5017 [NASA ADS] [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Auger, M. W., Treu, T., Bolton, A. S., et al. 2009, ApJ, 705, 1099 [Google Scholar]
 Auger, M. W., Treu, T., Bolton, A. S., et al. 2010, ApJ, 724, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Bacon, R., Accardo, M., & Adjali, L. 2010, in Groundbased and Airborne Instrumentation for Astronomy III, eds. I. S. McLean, S. K. Ramsay, & H. Takami, SPIE Conf. Ser., 7735, 773508 [NASA ADS] [Google Scholar]
 Barnabe, M., & Koopmans, L. V. E. 2007, ApJ, 666, 726 [NASA ADS] [CrossRef] [Google Scholar]
 Barnabè, M., Czoske, O., Koopmans, L. V. E., Treu, T., & Bolton, A. S. 2011, MNRAS, 415, 2215 [Google Scholar]
 Bevacqua, D., Cappellari, M., & Pellegrini, S. 2022, MNRAS, 511, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Mamon, G. A. 1982, MNRAS, 200, 361 [Google Scholar]
 Birrer, S. 2021, ApJ, 919, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Birrer, S., & Amara, A. 2018, Phys. Dark Univ., 22, 189 [Google Scholar]
 Birrer, S., & Treu, T. 2021, A&A, 649, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birrer, S., Amara, A., & Refregier, A. 2016, JCAP, 2016, 020 [CrossRef] [Google Scholar]
 Birrer, S., Treu, T., Rusu, C. E., et al. 2019, MNRAS, 484, 4726 [Google Scholar]
 Birrer, S., Shajib, A. J., Galan, A., et al. 2020, A&A, 643, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birrer, S., Shajib, A., Gilman, D., et al. 2021, J. Open Sour. Softw., 6, 3283 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, K., Castorina, E., & Simonović, M. 2020, ApJ, 892, L27 [Google Scholar]
 Bolton, A. S., Burles, S., Koopmans, L. V. E., Treu, T., & Moustakas, L. A. 2006, ApJ, 638, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Bolton, A. S., Burles, S., Koopmans, L. V. E., et al. 2008, ApJ, 682, 964 [Google Scholar]
 Bonvin, V., Millon, M., Chan, J. H. H., et al. 2019, A&A, 629, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cappellari, M. 2008, MNRAS, 390, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Cappellari, M. 2020, MNRAS, 494, 4819 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, Y.Y., van der Wel, A., Rix, H.W., et al. 2013, ApJ, 773, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, G. C. F., Fassnacht, C. D., Suyu, S. H., et al. 2019, MNRAS, 490, 1743 [NASA ADS] [CrossRef] [Google Scholar]
 Ciotti, L., & Bertin, G. 1999, A&A, 352, 447 [NASA ADS] [Google Scholar]
 Courbin, F., Bonvin, V., BuckleyGeer, E., et al. 2018, A&A, 609, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ding, X., Treu, T., Shajib, A. J., et al. 2018, ArXiv eprints [arXiv:1801.01506] [Google Scholar]
 Ding, X., Treu, T., Birrer, S., et al. 2021a, MNRAS, 503, 1096 [Google Scholar]
 Ding, X., Treu, T., Birrer, S., et al. 2021b, MNRAS, 501, 269 [Google Scholar]
 Dressel, L. 2012, Wide Field Camera, 3, 5 [NASA ADS] [Google Scholar]
 Dutton, A. A., & Treu, T. 2014, MNRAS, 438, 3594 [NASA ADS] [CrossRef] [Google Scholar]
 Dutton, A. A., Brewer, B. J., Marshall, P. J., et al. 2011, MNRAS, 417, 1621 [NASA ADS] [CrossRef] [Google Scholar]
 Falco, E. E., Gorenstein, M. V., & Shapiro, I. I. 1985, ApJ, 289, L1 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Freedman, W. L. 2021, ApJ, 919, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, W. L., Madore, B. F., Hoyt, T., et al. 2020, ApJ, 891, 57 [Google Scholar]
 Gavazzi, R., Treu, T., Rhodes, J. D., et al. 2007, ApJ, 667, 176 [Google Scholar]
 Gilman, D., Birrer, S., & Treu, T. 2020, A&A, 642, A194 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Golse, G., & Kneib, J. P. 2002, A&A, 390, 821 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gomer, M., & Williams, L. L. R. 2020, JCAP, 2020, 045 [CrossRef] [Google Scholar]
 Gomer, M. R., & Williams, L. L. R. 2018, MNRAS, 475, 1987 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
 Hao, C. N., Mao, S., Deng, Z. G., Xia, X. Y., & Wu, H. 2006, MNRAS, 370, 1339 [CrossRef] [Google Scholar]
 Huang, C. D., Riess, A. G., Yuan, W., et al. 2020, ApJ, 889, 5 [Google Scholar]
 Humphrey, P. J., & Buote, D. A. 2010, MNRAS, 403, 2143 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Kassiola, A., & Kovner, I. 1993, ApJ, 417, 450 [Google Scholar]
 Keeton, C. R., Kochanek, C. S., & Falco, E. E. 1998, ApJ, 509, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Kennedy, J., & Eberhart, R. 1995, Proceedings of ICNN’95 – International Conference on Neural Networks, 4, 1942 [CrossRef] [Google Scholar]
 Kenworthy, W. D., Scolnic, D., & Riess, A. 2019, ApJ, 875, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S. 2020, MNRAS, 493, 1725 [NASA ADS] [CrossRef] [Google Scholar]
 Kochanek, C. S. 2021, MNRAS, 501, 5021 [Google Scholar]
 Koopmans, L. V. E., Treu, T., Fassnacht, C. D., Blandford, R. D., & Surpi, G. 2003, ApJ, 599, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Koopmans, L. V. E., Bolton, A., Treu, T., et al. 2009, ApJ, 703, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J., Fisher, D. B., Cornell, M. E., & Bender, R. 2009, ApJS, 182, 216 [Google Scholar]
 Larkin, J., Barczys, M., Krabbe, A., et al. 2006, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, eds. I. S. McLean, & M. Iye, 6269, 62691A [Google Scholar]
 Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 363, 705 [Google Scholar]
 Merritt, D. 1985, AJ, 90, 1027 [Google Scholar]
 Millman, K. J., & Aivazis, M. 2011, Comput. Sci. Eng., 13, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Millon, M., Courbin, F., Bonvin, V., et al. 2020a, A&A, 640, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Millon, M., Galan, A., Courbin, F., et al. 2020b, A&A, 639, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morrissey, P., Matuszewski, M., Martin, D. C., et al. 2018, ApJ, 864, 93 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [Google Scholar]
 Oliphant, T. E. 2007, Comput. Sci. Eng., 9, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Osipkov, L. P. 1979, Pisma v Astronomicheskii Zhurnal, 5, 77 [Google Scholar]
 Pesce, D., Braatz, J., Reid, M., et al. 2020, H_{0}2020: Assessing Uncertainties in Hubble’s Constant Across the Universe, 20 [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poulin, V., Smith, T. L., Karwal, T., & Kamionkowski, M. 2019, Phys. Rev. Lett., 122, 221301 [Google Scholar]
 Riess, A. G., Yuan, W., Macri, L. M., et al. 2022, ApJ, 934, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Rusu, C. E., Wong, K. C., Bonvin, V., et al. 2020, MNRAS, 498, 1440 [Google Scholar]
 Saha, P., & Williams, L. L. R. 2006, ApJ, 653, 936 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, T., Treu, T., Birrer, S., et al. 2022, MNRAS, submitted [arXiv:2206.04696] [Google Scholar]
 Schneider, P., & Sluse, D. 2013, A&A, 559, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schneider, P., & Sluse, D. 2014, A&A, 564, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sérsic, J. L. 1963, Bol. Assoc. Argent. Astron. Plata Argent., 6, 41 [Google Scholar]
 Shajib, A. J. 2019, MNRAS, 488, 1387 [NASA ADS] [CrossRef] [Google Scholar]
 Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649 [Google Scholar]
 Shajib, A. J., Birrer, S., Treu, T., et al. 2020, MNRAS, 494, 6072 [Google Scholar]
 Shajib, A. J., Treu, T., Birrer, S., & Sonnenfeld, A. 2021, MNRAS, 503, 2380 [Google Scholar]
 Shi, Y., & Eberhart, R. 1998, 1998 IEEE International Conference on Evolutionary Computation Proceedings. IEEE World Congress on Computational Intelligence (Cat. No.98TH8360), 69 [CrossRef] [Google Scholar]
 Sluse, D., Chantry, V., Magain, P., Courbin, F., & Meylan, G. 2012, A&A, 538, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sonnenfeld, A. 2018, MNRAS, 474, 4648 [Google Scholar]
 Sonnenfeld, A., Treu, T., Gavazzi, R., et al. 2012, ApJ, 752, 163 [Google Scholar]
 Spingola, C., McKean, J. P., Auger, M. W., et al. 2018, MNRAS, 478, 4816 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Marshall, P. J., Blandford, R. D., et al. 2009, ApJ, 691, 277 [Google Scholar]
 Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201 [Google Scholar]
 Suyu, S. H., Treu, T., Hilbert, S., et al. 2014, ApJ, 788, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Suyu, S. H., Bonvin, V., Courbin, F., et al. 2017, MNRAS, 468, 2590 [Google Scholar]
 Treu, T., & Koopmans, L. V. E. 2002a, ApJ, 575, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Treu, T., & Koopmans, L. V. E. 2002b, MNRAS, 337, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Treu, T., Agnello, A., Baumer, M. A., et al. 2018, MNRAS, 481, 1041 [NASA ADS] [CrossRef] [Google Scholar]
 The pandas development team 2020, https://doi.org/10.5281/zenodo.3509134 [Google Scholar]
 Unruh, S., Schneider, P., & Sluse, D. 2017, A&A, 601, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van de Ven, G., Mandelbaum, R., & Keeton, C. R. 2009, MNRAS, 398, 607 [NASA ADS] [CrossRef] [Google Scholar]
 Van de Vyvere, L., Gomer, M. R., Sluse, D., et al. 2022a, A&A, 659, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van de Vyvere, L., Sluse, D., Gomer, M. R., & Mukherjee, S. 2022b, A&A, 663, A179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Meth., 17, 261 [Google Scholar]
 Waskom, M. L. 2021, J. Open Sour. Softw., 6, 3021 [NASA ADS] [CrossRef] [Google Scholar]
 Weijmans, A.M., de Zeeuw, P. T., Emsellem, E., et al. 2014, MNRAS, 444, 3340 [Google Scholar]
 Wertz, O., Orthen, B., & Schneider, P. 2018, A&A, 617, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, eds. S. van der Walt, & J. Millman, 56 [Google Scholar]
 Wong, K. C., Suyu, S. H., Auger, M. W., et al. 2017, MNRAS, 465, 4895 [NASA ADS] [CrossRef] [Google Scholar]
 Wong, K. C., Suyu, S. H., Chen, G. C. F., et al. 2020, MNRAS, 498, 1420 [Google Scholar]
 Xu, D., Sluse, D., Schneider, P., et al. 2016, MNRAS, 456, 739 [Google Scholar]
 Xu, D., Springel, V., Sluse, D., et al. 2017, MNRAS, 469, 1824 [NASA ADS] [CrossRef] [Google Scholar]
 Yıldırım, A., Suyu, S. H., & Halkola, A. 2020, MNRAS, 493, 4783 [Google Scholar]
 Yıldırım, A., Suyu, S. H., Chen, G. C. F., & Komatsu, E. 2021, A&A, submitted [arXiv:2109.14615] [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoo, J., Kochanek, C. S., Falco, E. E., & McLeod, B. A. 2006, ApJ, 642, 22 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Troubleshooting Chameleon profiles
A.1. Sérsicmatching modification
Initially we intended to use a Sérsic profile directly for the baryon mass component, matching the observed family of light distributions. Unfortunately, the Sérsic profile is computationally prohibitive to implement as a lens. Instead, we followed the practice of Dutton et al. (2011) and Suyu et al. (2014) to create a Chameleon profile which emulates the Sérsic profile by using a combination of two cored isothermal profiles. In the notation of Dutton & Treu (2014), Eq. 7 reads:
Dutton et al. (2011) provide a prescription to convert a given Sérsic profile into a Chameleon profile based on a polynomial fit. When we followed this prescription, we found we were unable to match our input Sérsic profile to the provided precision. Instead, we found that, at least for a Sérsic with n = 4, we had to manually modify the ratio of the two core radii, dividing α by a factor of 1.5. Only after this modification do our profiles match Fig. A.2 of Dutton et al. (2011). This effect is shown in Fig. A.1. We have not tested the general case. We suspect that there is a transcription error in the polynomial provided within the manuscript, but it is also possible that there is an error in our implementation, so we have made available a Python notebook which shows our calculations^{7}. At any rate, the Chameleon profiles we use in this work are constructed from our modified prescription and match Sérsic profiles to the precision shown in Fig. A.1.
Fig. A.1. Comparing the matching Chameleon profile to the input Sérsic profile. The modification we have made to the prescription better reproduces the corresponding figure from Dutton et al. (2011). 
A.2. Multiple Gaussian expansion
While troubleshooting the kinematics calculation within lenstronomy, we discovered a discrepancy can be introduced when the light profile is modeled as a Multiple Gaussian Expansion (MGE) to approximate the profile rather than the analytical profile. The MGE is used because it is simple to deproject, and a 3D mass distribution is required for kinematics. However, lenstronomy uses a 1D MGE for its kinematic calculation because it assumes spherical Jeans, which can lead to inaccuracies in deprojection when starting from an elliptical mass distribution. Namely, we found that the MGE introduces an excess of light at large radii in 3D when the lens is not circular. This can introduce a bias in the calculated σ_{v} at large radii, which can bias the apertureweighed v_{ap} by approximately 3% (see Fig. A.2). The effect worsens for larger radii and for more elliptical lenses (this example used q = 0.8). To avoid this, lenstronomy has been updated to include the 3D profile for the Chameleon profile, such that the MGE is no longer necessary, but this serves as a cautionary tale to beware of subtle effects when calculating kinematics of lens models.
Fig. A.2. Use of the 1D MGE introduces an excess of light in 3D at large radius for an elliptical lens. This causes the σ_{v} calculation to be biased by approximately 3% near the aperture radius of 1 arcsecond. 
Appendix B: Single fit example
As a demonstration of the MST, we plot the result from lensing fits to a single profile (created with a source placed at both z_{s}). Figure B.1 shows that the bestfit PEMD result is not equivalent to the convergence profile, but rather has been transformed by an MST. This results in a shallower slope than the truth. Kinematics are required to correct the effect through a determination of λ_{int}; the lensing data alone returns a biased value of H_{0}.
Fig. B.1. Example of a single profile, placed at each redshift and fit as a PEMD. Top: Convergence profiles as a function of radius; solid curves indicate input convergence profiles while dashed lines indicate convergence profiles corresponding to the PEMD fit. Bottom: residuals showing the differences between the input and PEMD fit profiles. Vertical lines correspond to the Einstein radii of the corresponding redshift. 
Appendix C: Fitting parameter bounds
The summarized parameters for lens construction and the uniform prior bounds used in this work are displayed in Table C.1.
Parameters for the mock creation and hierArc priors for this work.
Appendix D: Singleredshift test
Because the hierarchical framework is able to combine systems which span a range of impact parameters into an unbiased H_{0} result, we were curious if the range of spread was beneficial to the success of the method. We modified the set of 20 SLACSlike images from Sect. 4.1 by adding point sources with time delay information. For a set with no point sources, we use the remaining 20 SLACSlike images created in Sect. 4.2. As such, all lenses in this test have the same redshifts, with the case with no point sources simply corresponding to a different realization of profiles drawn from the same global distribution. We combine these sets together using the hierarchical framework and plot the results in Fig. D.1.
Fig. D.1. Hierarchical results from lens populations without point sources (purple) and with point sources (red) are plotted, as well as their combination (blue). Unlike Fig. 6, the red population has the same source redshift as the purple population. 
We find that the result is similar to Sect. 4.1, with the combination of populations successfully recovering H_{0} with 2.2% precision. One minor difference is that the intrinsic scatter σ(λ_{int}) is now consistent with zero; the framework captures the fact that this population of lenses has less intrinsic variation than the case where it is being probed at a wider range of different radii.
All Tables
All Figures
Fig. 1. Histograms depict the distribution of observed SLACS (purple) and TDCOSMO (orange) lens properties. Our simulated populations (shown as KDEs of corresponding colors) were created by constructing a large number of profiles and selecting a subset which simultaneously match the 5 target parameter distributions, indicated by the green borders. The simulated SLACS and TDCOSMO populations are the same profiles, only differing by a changed source redshift. For parameters where the two sets are identical, a black line is used instead of orange or purple. 

In the text 
Fig. 2. Redshift comparison with the observed lens population. Left: the distributions of source redshifts against lens redshifts for the populations of observed SLACS (purple) and TDCOSMO (orange) lenses. Our selections of source and lens redshifts are indicated as dashed lines. Right: these selections result in the dashedline Σ_{crit} values, which are consistent with the observed populations. 

In the text 
Fig. 3. Scaled radial quantity used by the hierarchical framework to parameterize radial dependence for MST, plotted for the observed population (histograms) and the simulated profiles used in this work (KDE curves). 

In the text 
Fig. 4. Convergence profiles for 20 analytical Chameleon+NFW profiles used in this paper (purple: SLACSlike; orange: TDCOSMOlike) alongside 16 profiles from hydrosimulations used in the TDLMC (green). Colored bands correspond to the 16th and 84th percentiles. Vertical colored lines indicate the median Einstein radius of the corresponding population. The diagonal black line indicates an isothermal slope for comparison. 

In the text 
Fig. 5. Images (log brightness) and residuals (in units of σ) of 10 systems at 2 source redshifts. Left: SLACS analog with z_{s} = 0.6 and no point sources, Right: TDCOSMO analog with z_{s} = 2.0. Lowerleft scale bars correspond to 1″. 

In the text 
Fig. 6. Results of the hierarchical framework using 20 Chameleon+NFW lenses for the population with SLACS source redshift and no time delays (purple), the population with TDCOSMO source redshift (orange) and the combination of both populations (blue). 

In the text 
Fig. 7. Effect of the number of lenses: the hierarchical result of a population of Chameleon+NFW lenses, combining the fits from the TDCOSMO and SLACS source redshifts. The blue is the same as that in Fig. 6, using 20 lenses from each redshift. The pink distribution corresponds to the same setup, but with 40 lenses each. 

In the text 
Fig. 8. Lensingonly values of H_{0} (i.e., λ_{int}) as a function of radial impact parameter defined in Eq. (4), which assumes this relationship is described as linear with Gaussian scatter. The SLACSz_{s} systems are indicated in red rather than purple because they are alternate versions of the same systems but with point sources and time delays (such that a value for H_{0} can be recovered and plotted). 

In the text 
Fig. A.1. Comparing the matching Chameleon profile to the input Sérsic profile. The modification we have made to the prescription better reproduces the corresponding figure from Dutton et al. (2011). 

In the text 
Fig. A.2. Use of the 1D MGE introduces an excess of light in 3D at large radius for an elliptical lens. This causes the σ_{v} calculation to be biased by approximately 3% near the aperture radius of 1 arcsecond. 

In the text 
Fig. B.1. Example of a single profile, placed at each redshift and fit as a PEMD. Top: Convergence profiles as a function of radius; solid curves indicate input convergence profiles while dashed lines indicate convergence profiles corresponding to the PEMD fit. Bottom: residuals showing the differences between the input and PEMD fit profiles. Vertical lines correspond to the Einstein radii of the corresponding redshift. 

In the text 
Fig. D.1. Hierarchical results from lens populations without point sources (purple) and with point sources (red) are plotted, as well as their combination (blue). Unlike Fig. 6, the red population has the same source redshift as the purple population. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.