The Hubble constant from strongly lensed supernovae with standardizable magnifications

The dominant uncertainty in the current measurement of the Hubble constant ($H_0$) with strong gravitational lensing time delays is attributed to uncertainties in the mass profiles of the main deflector galaxies. Strongly lensed supernovae (glSNe) can provide, in addition to measurable time delays, lensing magnification constraints when knowledge about the unlensed apparent brightness of the explosion is imposed. We present a hierarchical Bayesian framework to combine a dataset of SNe that are not strongly lensed and a dataset of strongly lensed SNe with measured time delays. We jointly constrain (i) $H_0$ using the time delays as an absolute distance indicator, (ii) the lens model profiles using the magnification ratio of lensed and unlensed fluxes on the population level and (iii) the unlensed apparent magnitude distribution of the SNe population and the redshift-luminosity relation of the relative expansion history of the Universe. We apply our joint inference framework on a future expected data set of glSNe, and forecast that a sample of 144 glSNe of Type~Ia with well measured time series and imaging data will measure $H_0$ to 1.5\%. We discuss strategies to mitigate systematics associated with using absolute flux measurements of glSNe to constrain the mass density profiles. Using the magnification of SNe images is a promising and complementary alternative to using stellar kinematics. Future surveys, such as the Rubin and \textit{Roman} observatories, will be able to discover the necessary number of glSNe, and with additional follow-up observations this methodology will provide precise constraints on mass profiles and $H_0$.


INTRODUCTION
The current expansion rate of the Universe, the Hubble Constant H 0 , anchors the scale and the age of the Universe. There is an ongoing debate about the precise value of H 0 , where some local distance ladder measurements based on calibration using Cepheids (e.g., Riess et al. 2021) are in significant statistical disagreement with measurements extrapolated from the cosmic microwave background (CMB) (e.g., Planck Collaboration et al. 2020; Aiola et al. 2020). Another distance ladder analysis based on calibration using the tip of the red giant branch (TRGB) stars results in a consistent measurement with the CMB (Freedman et al. 2020;Freedman 2021). This discrepancy either indicates unaccounted systematics in one or multiple measurements (e.g., Efstathiou 2020; Mortsell et al. 2021), or new physics beyond the stan-et al. 2020; Millon et al. 2020b). These measurements assumed particular forms of the mass density profiles of the deflector galaxies. The mass-sheet degeneracy (see Falco et al. 1985;Schneider & Sluse 2013, hereafter, MSD), an inherent transform leaving the lensing observables invariant while changing the time-delay prediction, poses limits in the precision of H 0 measurements in the absence of additional data. Birrer et al. (2020) introduced an additional degree of freedom to the mass density profiles to avoid constraining the lens model based on the specific form of the mass profiles previously chosen. Birrer et al. (2020) constrained the MSD solely by stellar kinematics observations of the deflector galaxy hierarchically on the deflector population level mitigating covariances among the assumptions of individual lenses. For the achieved 5% precision measurement of H 0 , Birrer et al. (2020) combined the 7 TDCOSMO lenses with 33 galaxy-galaxy lenses from the Sloan Lens ACS (SLACS) survey (Bolton et al. 2008;Shajib et al. 2021). The interpretation of the kinematics measurements are impacted by the mass-anisotropy degeneracy (Binney & Mamon 1982;Dejonghe & Merritt 1992) and mitigating this degeneracy requires assumption on the stellar anisotropy distribution or spatially resolved kinematics measurements (e.e., Cappellari 2008; Barnabè et al. 2011;Yıldırım et al. 2020). A forecast for future constraints using kinematics observations in breaking the MSD within the assumptions of the Birrer et al. (2020) analysis is provided by Birrer & Treu (2021).
An alternative to lensed quasars, as in fact anticipated in the original work by Refsdal (1964), are multiplyresolved gravitationally lensed supernovae (glSNe). glSNe are exquisite laboratories for fundamental physics, as well as astrophysical properties of the host and lens galaxies (see Oguri 2019, for a review of strong lensing of SNe and other explosive transients). We refer to Goobar et al. (2002), for example, for early explorations of cosmological parameter forecasts with hundreds of glSNe.
Although strongly lensed galaxies and quasars (QSOs) are more common than currently discovered occurrences of glSNe, glSNe have notable advantages, particularly if they are of Type Ia (SNe Ia). The luminosities of SNe Ia have a small dispersion after correcting for the relations with the lightcurve shape, observed color, and properties of their host galaxies, making them a "standardizable candle" (e.g., see Phillips 1993;Guy et al. 2007;Scolnic et al. 2018). Knowledge about the apparent magnitude, at the source redshift of the SNe, in the absence of any lensing effect allows us to directly measure the lensing magnification factor at the locations relevant to predict the time delays, breaking the MSD (see also e.g., Kolatt & Bartelmann 1998;Oguri & Kawano 2003;Foxley-Marrable et al. 2018). Thus, breaking the MSD does not require an a priori knowledge of the SN Ia absolute magnitude, and thus keeping the inference from time-delay cosmography independent from the local distance ladder calibration. Additionally, SNe Ia have a well-studied family of light curves with a well-defined maximum at ∼18 days from explosion (e.g. Yao et al. 2019;Miller et al. 2020) and hence, can be used for an accurate measurement of time delays, with significantly fewer follow-up observations than quasars. Recent simulations of lensed SNe Ia also find that the timedelay measurement is not impacted significantly from microlensing Huber et al. 2021). Moreover, since SNe fade away, we can obtain post-explosion imaging to validate the lens model (see e.g., Ding et al. 2021).
glSNe Ia are also complementary to lensed QSOs in the strategy for discovering these system. Owing to their small luminosity scatter, glSNe Ia can be discovered due to the lensing magnification increasing their brightness. This would not need highly spatially resolved observations, as is the case for lensed QSOs, important for testing potential biases from selecting high angular separation events. This was demonstrated in the discovery of the first resolved strongly lensed SN Ia, iPTF16geu (Goobar et al. 2017). At z = 0.409, the SNe was found to be 30 standard deviations too bright compared with the SNe Ia population, prompting space-based and laser guided star-adaptive optics (LGS-AO) follow-up. While iPTF16geu had a short time delay of ∼ 1 day (Goobar et al. 2017;More et al. 2017), hence, not ideal for measuring H 0 , the system could uniquely be used for a direct inference of the lensing magnification (Dhawan et al. 2020). Simulations of wide-field surveys like the Zwicky Transient Factory (ZTF) suggest a median time delay of ∼ 1 − 5 days (Goldstein et al. 2019;Wojtak et al. 2019) based on a magnification discovery channel. Upcoming deeper surveys, such as the Vera Rubin Observatory Legacy Survey of Space and Time (LSST), are expected to discover glSNe based on image multiplicity at fainter magnitudes, shifting the median time delay to ∼ 10 days (Wojtak et al. 2019), making glSNe Ia compelling probes of H 0 . glSNe also offer a unique opportunity to obtain time-delays from resolved spectroscopy, a method which requires very few epochs of observations Bayer et al. 2021).
With the advent of transient astrophysics and the anticipated discovery of more glSNe from current and future timedomain facilities, glSNe can play a major and complementary role in time-delay cosmography and beyond. In particular, the complementarity in constraining the MSD with magnification measurements in addition to stellar kinematics measurements allows one to rigorously check for systematics inherent in either of the two approaches as well as gain further statistical precision in the most limiting domain of time-delay cosmography to date.
Additionally, strong gravitational lensing systems are powerful probe of elliptical galaxy properties and evolution (e.g., Treu & Koopmans 2002;Auger et al. 2010;Shajib et al. 2021). Non-imaging data -such as the time delays, the stellar kinematics, or the image magnifications -provide additional constraints on the gravitational potential, or equivalently the mass distribution. In such studies, the adopted values of cosmological parameters can indeed have significant physical outcomes. For example, Blum et al. (2020) demonstrated that adopting the CMB-based H 0 value for the 7 TD-COSMO systems leads to galaxy mass distributions with a cored component in the dark matter profile. However, Shajib et al. (2021) combined only the stellar kinematics with the lens imaging data (without any time-delay measurement) to find that the deviation from the power-law profile in elliptical galaxies can also be caused by a higher normalization in the dark matter profile instead of having a cored component. The hierarchical analysis of Birrer et al. (2020) simultaneously constrained both the mass distribution in galaxies and H 0 for the first time from the combination of stellar kinematics and lensing information. glSN Ia will similarly provide simultaneous constraints on the galaxy mass distribution and cosmological parameters. Furthermore, a sufficiently large glSN Ia sample spanning a wide redshift range can provide direct insights into the evolution of massive elliptical galaxies.
In this paper, we aim to exploit the uniform behaviour of the population of SNe Ia to reduce uncertainties in the lens modeling arising from the MSD. We extend the hierarchical inference framework by Birrer et al. (2020) and incorporate SNe Ia apparent magnitudes, for both lensed and global unlensed populations, on the likelihood level in the cosmographic inferences. We perform forecasts at the same level of complexity as presented by Birrer & Treu (2021), now replacing the kinematics observables with SNe Ia brightness measurements, for different scenarios and highlight the key ingredients required to achieve an H 0 measurement with precision below 2%.
The paper is structured as follows. Section 2 provides a general review on the key concepts of time-delay cosmography, with a special focus on the MSD and approaches to constrain it. Section 3 defines the methodology, model parameterization, likelihood, and sampling approach used for the forecasts. Section 4 presents different forecast scenarios in regard to H 0 inferences. We discuss the key components and implications of this work in Section 5 and conclude in Section 6.
The formalism and inference schemes presented in this work are implemented in the open-source software HIER-ARC 1 and the scripts to reproduce the presented work is 1 https://github.com/sibirrer/hierArc publicly available 2 . Lensing calculations are performed with LENSTRONOMY 3 (Birrer & Amara 2018;Birrer et al. 2021).

TIME-DELAY COSMOGRAPHY WITH STRONGLY LENSED SNE
In this section, we review the principles of time-delay cosmography for lensing and time delays (Section 2.1). We then emphasize how an MSD affects the observables and thus the inference of cosmographic quantities, and specifically discuss the ability of glSNe in breaking the MSD with absolute lensing magnifications (Section 2.2).

Cosmography with strong lenses
The phenomena of gravitational lensing can be described by the lens equation, which maps the source plane coordinate β to the image plane θ as where α is the angular shift on the sky between the original unlensed position and the lensed observed position of an object. For a single deflector plane, the lens equation can be expressed in terms of the physical deflection angleα as where D s and D ds are the angular diameter distances from the observer to the source and from the deflector to the source, respectively. In the single lens plane regime, we can introduce the lensing potential ψ such that and the lensing convergence as The relative arrival time ∆t AB between two images θ A and θ B originated from the same source is where c is the speed of light, is the Fermat potential (Schneider 1985;Blandford & Narayan 1986), and is the time-delay distance (Refsdal 1964;Schneider et al. 1992;Suyu et al. 2010); D d is the angular diameter distance from the observer to the deflector. In the last line of Equation 5 we chose the notation ∆τ AB to describe the relative Fermat potential between two images. Constraints on the Fermat potential difference ∆τ AB and a measured time delay ∆t AB allow us to constrain the timedelay distance D ∆t . This absolute physical distance anchors the scale in the Universe within the redshifts involved in the lensing configuration. The Hubble constant is inversely proportional to the absolute scales of the Universe and thus scales with D ∆t as mildly dependent on the relative expansion history from current time (z = 0) to the redshifts of the deflector and the source.
2.2. The MST and the ability of lensing magnifications in breaking it

MST impact on time delays and imaging data
The mass-sheet transform (MST) is a multiplicative transform of the lens equation (Eqn. 1) given by which preserves image positions (and any higher order relative differentials of the lens equation) under a linear source displacement β → λβ (Falco et al. 1985). The term (1−λ)θ in Equation 9 above describes an infinite sheet of convergence (or mass), and hence the name mass-sheet transform.
Only observables related to the unlensed apparent source size, to the unlensed apparent brightness, or to the lensing potential are able to break this degeneracy. The convergence field transforms according to Thus, the same relative lensing observables can result if the mass profile is scaled by the factor λ with the addition of a sheet of convergence (or mass) of κ(θ) = (1 − λ). The different observables described in Section 2.1 relevant for time-delay cosmography transform by an MST term λ as follows: the image positions remain invariant the source position scales with λ as the Fermat potential scales with λ as and so does the time delay as When transforming a deflector profile with an MST, the inference of the time-delay distance (Eqn. 7) from a measured time delay and inferred Fermat potential transforms as Thus, the Hubble constant, when inferred from the timedelay distance D ∆t , transforms as (from Eqn. 8) Achieving precise and accurate constraints on the radial density profile required to measure H 0 necessitates external data and puts high demand on the precision and accuracy of those measurements and priors. We refer the reader to Section 2 of Birrer et al. (2020) for a discussion on interpretations of an MST in regard to a parameterized profile and physical limits of it.
There are two promising observables that have the ability to break the MST independent of the time delays: the stellar velocity dispersion measurements of the deflector galaxy and the absolute magnification measurement from knowledge of the apparent unlensed brightness of a source component.
For the remainder of this paper, we chose the convention of λ to be the mapping from a model prediction ignoring MST effects to the target prediction of the correct answer.

Stellar velocity dispersion
The stellar velocity dispersion of the main deflector galaxy is directly sensitive to the deflector potential. Joint lensing and kinematics measurements have been used to constrain the mass profiles of massive elliptical galaxies  and is the sole constraining anchor on the MST in the H 0 measurement by Birrer et al. (2020). The observed stellar velocity dispersion σ v scales with an MST as A fractional uncertainty in the velocity dispersion measurement σ v,obs , or model prediction σ v,model propagates to a fractional uncertainty in the MST as where we identified the target truth (measured) velocity dispersion with σ v,λ = σ v,obs and the model without the MST correction with σ v = σ v,model of Equation 17. Thus, an achievable 5% uncertainty in the measurement of σ v propagates to a 10% uncertainty in λ. Beyond the measurement uncertainty in σ v , projection uncertainties and degeneracies are present in the interpretation of the measurement, the model uncertainty. In particular, the mass-anisotropy degeneracy limits the precision, so only spatially resolved kinematics observations are able to break this secondary, but relevant, degeneracy (e.g., Yıldırım et al. 2020;Birrer & Treu 2021). Constraints on the radial extent of the mass profile with spatially resolved kinematics are possible. Equation 18 is, however, applicable for the covariant uncertainties among multiple measurements or integral-field unit spectroscopy. A forecast utilizing kinematic measurements of ground-and spacebased facilities on a larger sample of lenses within the same assumptions as Birrer et al. (2020) is presented by Birrer & Treu (2021).

Absolute lensing magnifications
The alternative to kinematics, and key element in the exploration in this work, are absolute magnification constraints (Kolatt & Bartelmann 1998;Foxley-Marrable et al. 2018). Absolute lensing magnifications, µ, change under an MST by A fractional uncertainty in the lensing magnification propagates to a fractional uncertainty in the MST as The observed magnification µ λ is the ratio where F obs is the observed flux of an image, and F unl is the unlensed apparent brightness of the object in the absence of the lensing effect. While measuring the observed flux of a lensed object is achieved to sub-percent precision on a regular basis, a lensing-independent measurement of µ λ requires, in addition, knowledge of the unlensed apparent brightness of the object in the same observational band as the measurement. We stress that the measurement of the lensing magnification does not require knowledge or calibration of the absolute luminosity, which is a key requirement in measuring H 0 with SNe Ia (e.g., Riess et al. 2019;Freedman et al. 2019). Only the probability distribution function of the apparent magnitude of the source at the redshift of the source is required.
The estimation of the MST scaling for a given lens model relevant for the time-delay prediction, and thus the measurement of H 0 , requires, in addition, an accurate lensing magnification prediction in accordance to the Fermat potential prediction. While the time-delay prediction is less susceptible to small scale model inaccuracies 4 as it relies only on an accurate lensing potential, the local magnification is impacted more significantly, as it relies on the second-order differentials of the potential. In addition to small scale dark matter structure, both along the line of sight (LOS) and within the main deflector as substructure, stellar microlensing is an additional source of lensing magnifications for sources of the size of exploding SNe (e.g., Dobler & Keeton 2006;Foxley-Marrable et al. 2018;Suyu et al. 2020).
We can approximately separate the different components entering the local magnification prediction into a smooth macro-model component µ macro and an additional perturbation by dark matter structure on milli-arcsecond scales, ∆µ milli , and stellar microlensing, ∆µ micro , as 5 µ local ≈ µ macro + ∆µ milli + ∆µ micro .
Milli-lensing depends on the halo substructure and on the line-of-sight abundances of small field halos. Stellar microlensing depends on the local projected stellar surface density and can vary significantly from lens to lens and from image position to image position.
The relative difference in λ for either an infinitesimal change in the apparent unlensed magnitude, δF unl , a change in the lensed observed flux, δF obs , the model predicted magnification, δµ model , or the physical cause of local milli (micro) lensing, δµ milli (δµ micro ), while keeping all other quantities fixed, can be expressed as (Eqn. 20,21,22) In words, while keeping all other parameters fixed, an increase in F unl leads to an increase in λ, an increase in F obs leads to a decrease in λ, an increase in the lensing magnifications µ macro , ∆µ milli and ∆µ milli leads to an increase in λ. On the other hand, errors in the measurement or estimation of these quantities result in shifts of λ in the opposite direction.
The intrinsically small scatter of Type Ia supernovae is a well suited population to constrain the MST. An intrinsic scatter of 10% in the peak brightness after lightcurve width and color corrections  allows one, at least in principle, to constrain the MST to 5% in the absence of other uncertainties. Thus, glSNe are not only able to provide precise time-delay measurements due to their transient and well characterized nature, but at the same time Type Ia or any other standardizable form of SNe allows one to constrain the currently dominating error budget of time-delay cosmography, the MST.
The constraints on the MST rely on precise and accurate determinations of all the parameters listed in Equation 23. The uncertainty in apparent unlensed brightness of SNe F unl , milli-lensing ∆µ milli contribution and the microlensing effect ∆µ micro are the dominant uncertainty components in constraining the MST. Systematic limitations in the usage of glSNe relate to dust extinction impacting the flux measurement F obs , and selection effects related to milli-and microlensing. We will review limitations and systematics of glSNE in breaking the MST in Section 5.

METHODOLOGY
In this section, we describe the methodology to measure H 0 from a set of glSNe by constraining the MST with the apparent magnitude distribution of an unlensed SNe sample. We layout the model assumptions and define the hyper-parameters governing the cosmological expansion, SNe brightness distribution, and the mass profiles of the lensing galaxies. Furthermore, we detail the implementation of the likelihood for the different observations that allow us to efficiently perform a joint hierarchical sampling of the posteriors. We describe the SNe Ia population assumptions and analysis in Section 3.1, and the glSNe population assumptions and analysis in Section 3.2. Separately, we discuss the impact and the treatment of LOS structure in Section 3.3. In Section 3.4, we state the joint hierarchical inference problem based on the previous parts of this section. In Section 3.5, we provide an approximate analytical error propagation. The methodology presented here, in terms of parameterization and likelihood calculation, is implemented in the open-source software HIERARC.

Supernovae of Type Ia population
We focus in this work on the SNe Ia population. Here we describe the SNe Ia magnitude-redshift relation and the likelihood for the unlensed SN Ia sample in Section 3.1.1. In Section 3.1.2, we then state the specific model and likelihood assumptions for the combined lensed and unlensed SNe Ia samples that we implement in this work.

Unlensed SN Ia sample
SNe Ia can be standardized as precise (relative) distance indicators. This involves well-known corrections for their luminosity-width and luminosity-color relations (Guy et al. 2007;Guy et al. 2010). In addition, the SN Ia inferred luminosity needs to be corrected for its dependence on the host galaxy properties (e.g. stellar mass) (refer to Scolnic et al. 2018, for details of the standardization procedure and bias corrections). The distance modulus, µ B , relates the standardized apparent magnitude of an unlensed SN Ia, denoted as m * B , with the absolute magnitude M B as where D L (z) is the luminosity distance from the observer to the redshift of the SNe, and µ dist is the distance modulus.
There are several large samples of the cosmological SNe Ia in the literature, e.g. Pantheon , Joint Light-curve Analysis (JLA; Betoule et al. 2014), or the Dark Energy Survey (DES; Abbott et al. 2019). For our analyses we use the largest, up-to-date compilation, i.e. the Pantheon SN Ia dataset ) For these SN Ia samples, covariances in the calibration parameters ξ sys and their evolutionary trends need to be taken into account.
Lensing magnifications, be it weak lensing from largescale structures, strong lensing from massive deflectors, or micro-lensing from stars, change the flux or apparent magnitude by where µ is the unsigned absolute magnification. For a single SN, we can formally write down the likelihood of an observed peak brightness F obs given a luminosity distance D L , lensing magnification µ, and absolute magnitude M B as the likelihood of the data given a flux prediction F model while marginalizing over calibration and other uncertainties, such as dust extinction, or uncertainty in the peak time, as For a sample of SNe, denoted as D SNe , different procedures and models have been employed for the parameterization, calibration, and marginalization of systematic errors on the population level and we refer to the relevant work for details (e.g., Betoule et al. 2014;Scolnic et al. 2018;Abbott et al. 2019). For this work, when combining such a supernovae sample with glSNe, we require the likelihood L(D SNe | π, ξ SNe ) of the global data set given the cosmological prediction of the luminosity distances with parameters π and the intrinsic brightness distribution of the SNe population ξ SNe .
To facilitate the evaluation of the likelihood, for example, Scolnic et al. (2018) compressed the marginalization in a Gaussian covariance matrix across all the measured apparent magnitudes in the SNe sample as where ∆m is the difference in the observed and predicted apparent magnitude of a non-evolving intrinsic mean brightness of the SNe population, n sn is the length of the data vector, and Σ cov is the error covariance matrix when marginalized over the systematics variables in Gaussian form.

Model parameterization and assumptions with pivot magnitude
We assume, for simplicity of this work, that the intrinsic peak brightness distribution is redshift independent. This is typically assumed in cosmological analyses with SNe Ia, based on comparisons of spectroscopic and photometric properties of local and high-z SNe Ia (e.g. Petrushevska et al. 2017, and other studies of high signal-to-noise data of high-z SNe Ia).
To obtain the absolute luminosity of an unlensed SN Ia at the redshift of the glSNe in our sample would require an independent calibration of the absolute luminosity distance, e.g. as done for the distance ladder. However, since we want to derive a relative magnification at the lensed source redshift, we can use the apparent magnitude of the unlensed SNe Ia from the cosmological sample and infer it at the redshift of the lensed SN. For this, we replace using an M B term with a apparent magnitude at a specific redshift, z pivot , m p m sn (m p , z) = m p + 5 [log 10 D L (z) − log 10 D L (z pivot )] .
(28) This parameterization results in a likelihood which is only dependent on relative distance ratios without the need of external data or constraints on top of a population of observed peak brightness of SNe. We describe the intrinsic distribution of apparent peak brightness at the pivot redshift p(m p ) by a Gaussian in astronomical magnitude space with a mean m p and width σ(m p ).
With these simplifications, we can write the likelihood for a single SNe as where F (m p ) is the shorthand form of the model predicted flux given an apparent magnitude m sn calculated by Equation 28 from m p and the luminosity distance ratio, and then turned into flux units of the observations while considering the lensing magnification (Eqn. 25). σ obs is the Gaussian error in the flux measurements and the term p m p | m p , σ(m p ) in the equation above describes the likelihood of a specific pivot magnitude to be drawn from the Gaussian distribution N (m p , σ(m p )).
For an ensemble of SNe characterized with the likelihood of Equation 27, the marginalization over a Gaussian distribution in m p is analytic and can directly be folded in the error covariance matrix as and the marginalized likelihood is given by

Deflector population
We first discuss general considerations about the deflector parameterization and the necessary degrees of freedom to allow for an accurate recovery of the time-delay prediction (Section 3.2.1). We then formulate the inference problem and the general form of the joint likelihood of the imaging data, time delays, and lensed SNe peak brightness observations provided by a glSNe (Section 3.2.2). Lastly, we provide a Gaussian approximation of the likelihood for a fast marginalization and efficient evaluation (Section 3.2.3).

Deflector parameterization and measurements
The uncertainty in the deflector mass distribution dominates the current error budget in the H 0 inference, a statement directly reflecting the MST. A popular model describing strong gravitational lensing imaging data on galaxyscale lenses is the power-law elliptical mass distribution (PEMD; Barkana 1998;Tessore & Metcalf 2015) combined with an external shear component. The popularity of the PEMD+shear model is a consequence of its ability to describe the data sufficiently well while keeping the degrees of freedom in the deflector model to a computationally affordable number.
Although considered simplistic, the PEMD+shear model's degrees of freedom can describe the primary azimuthal and radial observables. However, the observable in the radial direction are related to the third order differential of the lensing potential, while the parameterization of the PEMD profile explicitly assumes a one-to-one connection between the observable invariant quantity and the mass density at the position of the Einstein ring, leading to over-constrained mass profiles and potentially biased inferences in the radial density profile, and subsequently H 0 (see e.g. Kochanek 2002;Sonnenfeld 2018;Kochanek 2020Kochanek , 2021Birrer et al. 2020;Birrer 2021).
To mitigate possible over-constraints on the internal mass density profile, Birrer et al. (2020) added an additional degree of freedom with an MST on top of the PEMD+shear profiles of the TDCOSMO sample. A PEMD+shear+MST profile has the adequate degrees of freedom at and around the Einstein ring, where the multiple images appear, to estimate the relative Fermat potential. Higher-order differentials are subdominant in the effect on the predicted time delays (e.g., Sonnenfeld 2018). In addition, since the constraints on the MST from the lensing magnification are directly derived at the region relevant for the time-delay prediction, potential inadequacies of the PEMD+shear+MST profile further outwards or towards the center of the deflector do not impact the accuracy in the inferred time-delay prediction, and thus H 0 measurement.
In this work, we assume that the population of lenses can be described by a PEMD+shear+MST profile. Furthermore, we assume that the PEMD+shear parameters can be measured accurately for each lens individually from the imaging data without population covariances, and that the MST parameter λ transforming the internal density profile of the main deflector, denoted as λ int , follows a Gaussian distribution with mean λ int and sigma σ(λ int ). We refer to Section 3.3 for a discussion on different MST components. We highlight that there can be physical covariances between the PEMD parameters and λ int , as well as among the physical projected scale and λ int . A possible physical projection dependence has been accounted for by Birrer et al. (2020) with an explicit parameterization of λ int as a function of the ratio of deflector half-light radius relative to the Einstein radius. In this work, for the purpose of providing a forecast, we do not include secondary dependencies and covariances of λ int , and instead refer to Birrer et al. (2020) for the radial dependence as well as to Wagner-Carena et al. (2021) for a general treatment of lens model hyper-parameter inferences within a hierarchical framework.

glSNe inference
From the imaging data, I, we can measure the lens model parameters within our model assumptions, ξ pl , which in turn provide the Fermat potential differences between the multiple images ∆τ pl and the lensing magnifications µ pl at the position of the appearances of the glSN. With measured relative time delays ∆t and a model providing values for ∆τ pl , λ, and D ∆t , we can predict the time delay and evaluate the time-delay likelihood of the data given the model. From the same lens model, we can compute the likelihood of the mea-sured glSNe brightness F given the model prediction of µ pl , λ, and m sn .
The joint likelihood L(I, ∆t, F | D ∆t , m sn , λ) of the imaging, time delay, and flux measurements -given the relevant parameters of the hierarchical inference, D ∆t , m sn , and λ -can be written as product of the likelihoods of the different independent data sets where we explicitly marginalized over the magnification and Fermat potential parameters µ pl and τ pl . To describe the imaging data I and to compute the likelihood at the pixel level, we require a lens model ξ pl and a model of all the light components ξ light . We can describe the imaging likelihood and prior product on ∆τ pl and µ pl as Here, ∆τ (ξ pl ) and µ(ξ pl ) are unique functions of ξ pl , and the ∂(∆τ pl , µ pl )/∂(ξ pl , ξ light ) is the Jacobian determinant. This means that the likelihood and prior product of Equation 33 can be computed by sampling ξ pl from the posterior L(I | ξ pl , ξ light )p(ξ pl , ξ light ), and evaluating for the posterior sample of the quantities ∆τ pl (ξ pl ) and µ pl (ξ pl ). For the modeling choices and posterior sampling when marginalizing over complex source structure, we refer to previous work (e.g., Suyu et al. 2009;Birrer et al. 2015).

Gaussian likelihood approximation
Until this point in this subsection, we did not make any assumption on the form of the likelihood (Eqn. 32) nor on the shape of the imaging modeling posteriors (Eqn. 33). To facilitate the calculation of the likelihood in Equation 32, we approximate the likelihood in Gaussian form. In particular, we write the imaging likelihood from Equation 33 as where ∆ ∆τ µ ≡ (∆τ −∆τ 0 , µ−µ 0 ) with (∆τ 0 , µ 0 ) being the maximum likelihood estimator, n ∆τ µ is the length of the vector ∆ ∆τ µ , and Σ ∆τ µ is the error covariance matrix describing the Gaussian uncertainties in the measurement from the imaging data. The Gaussian form of the time-delay likelihood L(∆t | D ∆t , λ, ∆τ pl ) of Equation 32 reads where n ∆t are the number of relative time delay measurements, Σ ∆t is the relative time-delay measurement error covariance matrix, and ∆ ∆t is the difference between the predicted time delay (Eqn. 5 including an MST term of Eqn. 14) and measured time delay ∆t. The Gaussian form of the flux amplitude likelihood where n F is the number of flux measurements, Σ F is the glSNe flux measurement error covariance matrix, and ∆ F is the difference between the predicted peak flux of the SNe and measured flux F . The marginalized likelihood of Equation 32 with the Gaussian approximations for the individual likelihood components (Eqn. 34, 35, 36) is a Gaussian integral.
We can join the data vector of the time delays ∆t and fluxes F as d ∆tF ≡ (∆t, F ) and write the joint measurement covariance matrix as In this forecast, we assume no covariant measurement uncertainties between the time delays and the micro-lensing impact on the magnification. It has been shown that the microlensing effect on time-delay measurements can be mitigated with early phase multi-color light curves when the microlensing effect is achromatic ). At the same time, we can transform the covariance matrix of the imaging posteriors Σ ∆τ µ into the data vector space, resulting in where 1 n∆τ and 1 nµ are vectors of size n ∆τ and n µ , respectively, with 1 at each element. The joint likelihood of Equation 32 is then given by with n ∆tF ≡ n ∆t + n F , and the total error covariance matrix being the sum of the measurement covariance and the marginalized uncertainty in the model and ∆ ∆tF being the difference of the data vector d ∆tF and the model prediction.

LOS mass distribution
Mass over-or under-densities along the LOS of the strong lensing system cause, to first order, shear and convergence perturbations. Reduced shear distortions do have a measurable imprint on the azimuthal structure of the strong lensing system (see e.g., Birrer 2021) while the convergence component of the LOS, denoted as κ ext is equivalent to an MST, and thus not directly measurable from imaging data. The total MST, the relevant transform to constrain for an accurate cosmography and H 0 measurement, is the product of the internal and external MST (e.g., Schneider & Sluse 2013;Birrer et al. 2016Birrer et al. , 2020 λ The lensing kernel impacting the linear distortions, both shear and κ ext is different from the standard weak lensing kernel (McCully et al. 2014(McCully et al. , 2017Birrer et al. 2017Birrer et al. , 2020Fleury et al. 2021). The lensing kernel can be described as the product of three different angular diameter distances entering D ∆t in Equation 7 Fleury et al. 2020), and thus κ ext can be described as the product of the individual kernels entering Equation 7 as where κ d is the weak lensing effect from the observer to the deflector, κ s from the observer to the source, and κ ds from the deflector to the source, respectively .
Alternatively, but equivalently, the kernel can be described in the multi-plane formalism with the main deflector included, while keeping the Born approximation in between (e.g., Birrer et al. 2017;Fleury et al. 2021).
The LOS lensing contribution can be estimated by tracers of the large-scale structure, either using galaxy number counts (e.g., Greene et al. 2013;Rusu et al. 2017), or weaklensing measurements (Tihhonova et al. 2018). These measurements, paired with a cosmological model including a galaxy-halo connection are able to constrain the probability distribution of κ ext to few per cent per LOS.
For an accurate measurement of H 0 , the combined internal and external MST of Equation 42 is required. Since glSNe magnification is directly probing the combined λ, the LOS contribution effectively only adds a scatter in the inference and an accurate overall population selection function is not required (see Birrer et al. 2020, for the same argument using kinematics to break the MST). The overall lensing selection function is only relevant when demanding a physical interpretation of the internal and external contributions separately.
In this work, for practical simplicity but without impact on expected biases or uncertainty budget, we assume a Gaussian scatter in κ ext of 0.03 with a population mean at zero along the LOS's of the lenses.

Hierarchical analysis and sampling
Our goal is to jointly infer and marginalize over population hyper-parameters in the SNe distribution, lensing deflector profiles, and cosmological parameters, given the joint data set of lensed and unlensed SNe, and MST-invariant lensing quantities from imaging data. We follow the same approach as Birrer et al. (2020), except that we add the SNe magnification likelihood instead of the stellar kinematic one, and as an external data set we are using a sample of unlensed SNe instead of a sample of galaxy-galaxy lenses with measured kinematics.
We want to calculate the probability of the cosmological parameters, π, given the joint data set, p(π | {D i L } N , D SNe ), where D i L is the data set of an individual strong lens (including imaging data, time-delay measurements, SNe flux measurement, and LOS properties), N is the total number of lenses in the sample, and D SNe is a SNe data set.
In addition to π, we introduce ξ pop that incorporates all the additional population-level model parameters not yet marginalized over the individual data sets including their covariant impact on the likelihoods of individual lenses. Using Bayes' rule and considering that the data of each individual lens D i is independent, we can write: Table 1 summarizes the hyper-parameters describing the cosmological parameters, the SNe brightness distribution as well as the lens population that we are sampling hierarchically. We also state the parameter priors we employ in the forecast. We refer to Birrer et al. (2020) for the formal approximation we are making in the Bayesian analysis while treating other lens model parameters independently among the different lenses and to Wagner-Carena et al. (2021) to a hierarchical analysis inferring a wider range of lens model hyper-parameters.
The likelihood of an individual lens for a given set of hyper-parameters, L(D i L | π, ξ pop ), is given by the integral of the individual parameters according to the specified distribution of the hyper-parameters where p(ξ | ξ pop ) is the distribution function of the individual parameters ξ as specified by the population parameters ξ pop , and L(D i L | π, ξ) is the likelihood specified by Equation 32 and its Gaussian form (Eqn. 39) when stating the angular diameter distances as a function of the cosmological parameters π. The same statement as for the lens likelihood (Eqn. 44) applies for the SNe sample likelihood L(D SNe | π, ξ pop ). The marginalization in L(D SNe | π, ξ pop ) goes over the supernovae brightness distribution hyper-parameters m p and σ(m p ). We note that the SNe distribution parameters are shared for both the SNe population likelihood and the individual lens likelihoods, as well as the cosmological parameters relevant to describe the relative expansion history. The absolute scales of the Universe, stated in the form of H 0 , only enter explicitly in the time-delay likelihood.

Analytic error propagation
Before we present the forecast and results with the full hierarchical sampling and propagating of the covariances in the model described in Section 3.4, we also provide an analytic, simplified, approximate error propagation. This calculation is easily accessible, fast to compute, and provides valuable insights in the relative importance of different uncertainty components impacting the final H 0 constraints.
To first order, the relative H 0 uncertainty, σ H0 /H 0 , comprises of the uncertainty in the population mean of the MST parameter 6 , λ, and the uncertainty when performing an uncorrelated error propagation when fixing λ = const as 6 including internal and external MST effects In the following, we approximate the uncertainty budget for the distinct terms in Equation 45. For simplicity of this analysis, we assume that for all lenses, and all images, the uncertainty terms are identical. In practice, and in the full inference, inverse uncertainty weighting must be considered.

Uncertainty terms in the MST
The first term on the right-hand side of Equation 45 above can be determined with absolute lensing magnifications. The population level uncertainty in λ can, to first order, be expressed as the uncertainty in the population mean of the apparent unlensed brightness m p 7 , which is covariant among all lenses, the uncertainty in the relative expansion history translating the apparent magnitude of the distribution of the external SNe sample to the glSNe source redshift, and uncorrelated measurement uncertainties for each individual lens (Eqn. 20 and 23) as where N lens is the number of lens systems. Furthermore, for simplicity, we assumed equal precision in the individual relative magnification measurements in Equation 46 above for each lens. The relative magnification uncertainty per lens with fixed source population mean m p , can be written following Equation 23 as 7 The differential in logarithmic astronomical magnitude m in regards to relative linear flux I is I∂m/∂I = −2.5 log 10 (e) ≈ −1.086. Thus small scatter described in astronomical magnitudes are approximately the same scatter in relative flux.
The first term on the right hand side of the equation above is the intrinsic scatter in the standardizable source, the second is the flux measurement uncertainty, and the following ones are the different scales of the lensing effect. The factor 1/4 comes from the fact that we consider quadruply lensed quasars as this approximation assumes the random errors in the milli-and micro-lensing effects to be uncorrelated among the different images. The macro-model magnification uncertainties are covariant and thus we omit the factor 1/4 in the approximation.

Time-delay and Fermat potential uncertainties
The second term on the right hand side of Equation 45 encompasses all other sources of uncertainties not related to global inference shifts due to the MST. In particular, this involves uncertainties in the time-delay measurements, the Fermat potential uncertainty for a specified mass profile family (in our case PEMD+shear) from high-resolution imaging data, and the random uncertainties in the LOS convergence estimates and the internal MST. In addition, we include in this second term uncertainties in the relative expansion history that translate the angular diameter distance measurements to the lensing system, D ∆t , relative to to the scales at current time, and thus H 0 , which we denote as σ(H 0 /D z=SL ).
In terms of distance measurements, we can approximately write where the relative time-delay distance measurement uncertainty can be estimated by the relative Fermat potential uncertainties from imaging modeling and the relative time-delay uncertainties and the scatter and random uncertainty in λ coming from the internal and external scatter, which can be approximated as The time-delay distance uncertainty per lens (Eqn. 49) is, to first order, a weighted product of all the different images. The random uncertainty in the MST acts as a noise term for the individual distance measurements for each lens.

FORECAST
Having formulated the methodology and parameterization in the previous sections, we perform different forecast scenarios based on predicted number of glSNe, quality of measurements and systematics effect. In Section 4.1 we state the expected number of glSNe and time-delay measurements and our assumptions on milli-and micro-lensing effects in the magnification. In Section 4.2 we state the lens model, source configuration and uncertainties expected from imaging data on the Fermat potential and magnifications. In Section 4.3 we present the scenario for current and future unlensed SNe data sets. Finally, in Section 4.4 we present the inference results for the different forecast scenarios.

Lens population, time-delay and magnification uncertainties
In this work, we focus on the discoveries expected by the Vera Rubin Observatory Legacy Survey of Space and Time (LSST). We do not perform an independent forecast and derive our fiducial forecast scenario based on previous work in the literature. (2017); Goldstein et al. (2018) estimated, based on the catalogue by Oguri & Marshall (2010), the number of glSNe Ia to be up to 500-900 in 10 years of LSST with unresolved photometric magnification detection where the brightest SN image reaches a peak apparent i-band magnitude of 22.15 or brighter. Wojtak et al. (2019) compared two different discovery techniques, by magnifaction and resolved image multiplicity and estimated the annual discovery rate with LSST to be 61 with magnification, 44 with resolved image multiplicity and 89 in hybrid discovery scheme.

Goldstein & Nugent
It has been noted that lensed supernovae found via image multiplicity exhibit longer time delays and larger image separations making them more suitable for cosmological constraints than their counterparts found via magnification (Wojtak et al. 2019;Huber et al. 2019). Huber et al. (2019) finds, when restricting the expected time-delay measurement to a minimum precision of < 5% and an accuracy of < 1% (if based solely on LSST observations) would reduce the number of lensed type Ia supernovae to about 1 per year. This rate can be increased by a factor of 2-16 by employing other instruments for follow-up observations. Beyond LSST, for example, Pierel et al. (2021) predicts that the Roman observatory will discovery ∼ 11 glSNe Ia. With follow up efforts in measuring the time delays of the sub-sample restricted on the most promising time-delay measurements (Huber et al. 2019), LSST+follow up is able to provide < 1% overall statistical precision on H 0 form the time-delay uncertainties of 20 glSNe Ia (e.g., Suyu et al. 2020).

Milli-and micro-lensing
Milli-and micro-lensing effects on the magnification of the glSNe can significantly impact the ability of glSNe to be used as standardizable candles. Milli-lensing, an effect caused by dark subhaloes of the main deflector or along the line of sight (Dalal & Kochanek 2002;Gilman et al. 2020a;Hsueh et al. 2020, e.g.,), or baryonic effects (e.g., Hsueh et al. 2016;Gilman et al. 2017). Flux ratio anomalies at the ∼ 10% level have been studied and used to constrain dark matter properties with quardruply lensed quasar flux ratio anomalies (e.g., Gilman et al. 2020a;Hsueh et al. 2020). For physical source size of SNe, Kelly et al. (in prep) estimated for SN Refsdal (Kelly et al. 2015) about a ∼ 10% scatter from milli-lensing based on the forward modeling methodology by Gilman et al. (2019Gilman et al. ( , 2020a. Microlensing caused by stars or other compact objects in the foreground lens-ing galaxy or along the line of sight can be a more significant limit to the standardization of glSNe. Microlensing can independently magnify or de-magnify individual images of the background source (Dobler & Keeton 2006;Bagherpour et al. 2006), introducing scatter into the shape and amplitude of the resulting light curves. The effect of microlensing on each lensed image depends on the local smooth lensing properties (convergence κ, shear γ) and the stellar (or compact) projected mass fraction κ * /κ. For example, Schechter & Wambsganss (2002) investigated stellar micro-lensing effects on lensed quasars at image magnifications of µ ∼ 10 with moderate compact object mass fractions and showed that for such scenarios, the expected micro-lensing scatter can result in more than an astronomical magnitude.
Foxley-Marrable et al. (2018), with the aim of assessing glSNe Ia to be standardizable in the same spirit as this work, evaluated the effect of microlensing on glSNe Ia for various image configurations. They found that there are regions of parameter space where the effect of microlensing is suppressed enough for the glSN Ia to be standardizable. Specifically, regions of low κ, γ and high s are subject to microlensing scatter of σ ML ∼ 0.15 in astronomical magnitude, particularly at early times. Physically, this corresponds to asymmetric configurations with at least one image located far outside the Einstein radius, which will experience the least amount of microlensing.
When Foxley-Marrable et al. (2018) combined their microlensing models with the glSNe Ia catalogue from Goldstein & Nugent (2017), they predicted that ∼ 22% of the ∼ 930 glSNe Ia to be discovered by LSST will be standardizable (σ ML ∼ 0.15 or below for at least one image). The standardizabe sample has a median maximal time delay of 44 days and consists of 5:1 ratio of doubles vs quads. Foxley-Marrable et al. (2018) further concluded that from their sample of 650 glSNe Ia, of which accurate time delays can be measured, the MSD can be broken at the 0.5% level when considering microlensing and intrinsic scatter of the SNe as the source of statistical uncertainties.

Specific numbers and uncertainties of this forecast
Overall, restricting the follow-up effort to a considerably smaller number than the overall expected discoveries optimized to derived time-delay precision and accuracy, LSST is expected to provide sufficient statistical precision on time delays with a sub-percent error budget on final H 0 constraints. However, using glSNe for standardizable magnification constraints may require a larger and potentially different subset of the glSNe Ia population to be further investigated with follow-up efforts. Given the mass profile uncertainties are at the 10% level for individual lenses, we consider in this forecast a scenario with an extended sample of glSNe Ia beyond the subset of Huber et al. (2019); Suyu et al. (2020) with lower precision time-delay measurements, including both glSNe with shorter time delays as well as fainter images.
In this forecast, we design a scenario where time-delay precision and standardizable nature of glSNe Ia can be utilized.
We stress that time-delay measurement and flux standardization do not necessary need to come from the same lenses 8 .
We chose a lens population roughly following Foxley-Marrable et al. (2018). In total, we perform our forecast with 144 glSNe, among which 24 are quads and 120 are doubles. For the quad population, we split the sample in 8 crosses, 8 cusps and 8 fold configurations 9 . The doubles we split into three different configurations each consisting of 40 systems.
For the redshift distribution, we assume a uniform distribution of the deflector redshift, z lens , between z = 0.1 and z = 0.5, and for the source redshift, z source , a uniform distribution in U[z lens + 0.2, 1.], similar as the distribution by Huber et al. (2019); Suyu et al. (2020) restricting to the brighter population for both accurate time-delay and flux measurements. We stress the importance of rapid spectroscopic follow-up to confirm the SNe type and we assume that the follow-up has been acquired for the SN sample and the SNe have been robustly typed.
For the time-delay measurement, we assume that the light curves can be resolved in follow up observations and the relative time delays can be measured with a precision of 2 day per image pair 10 . Along with spectroscopy obtained for the typing, these cadenced observations provide further evidence to distinguish the normal SNe Ia from peculiar subtypes (e.g., see Taubenberger 2017, for review), since fast-declining and super-Chandra subtypes do not show a second maximum in the NIR, unlike normal SNe Ia. The presence of an NIR second maximum was further confirmation that iPTF16geu is a normal SN Ia (Dhawan et al. 2020). In addition to precise time delays, obtaining resolved photometry, in multiple wavebands, is crucial for constraining the extinction properties. We assume that similar to the case for iPTF16geu, there are cadenced observations in multiple optical and NIR filters to constrain the extinction in the host galaxy and the individual lines of sight in the lens for each image (Dhawan et al. 2020). Accounting for extinction correction, the magnification is inferred robustly with small uncertainties, making it a subdominant contribution to other sources.
For the flux uncertainty at peak brightness for the individual images, we use an effective relative magnitude uncertainty σ eff (m) that includes possible uncertainties from small scale milli-and micro-lensing effects This is a practically convenient noise definition when assuming Gaussian uncorrelated error in terms of the uncertainty relevant to constraining the MST. Beyond the intrinsic scatter in the SNe population, σ(m p ), and the uncertainty in the macro-model magnification, σ(µ macro ), the term in Equation 51 above can play a dominant role in the uncertainty budget and is by itself uncertain given the current rare discoveries and follow-up data of glSNe systems. We separate the σ eff (m) term for the different images into one image denoted as the standardizable one, σ eff,std (m), and all other images denotedas the microlensing dominated ones, σ eff,ML (m). For σ eff,ML (m) we assume a scatter of one magnitude, σ eff,ML (m) = 1.0, making most images of glSNe inefficient probes of the mass profile.
For the 'standardizable image' we perform three different scenarios for σ eff,std (m). The first scenario, denoted as IDEAL, sets σ eff,std (m) = 0 for all measurements, assuming no milli-and micro-lensing effects and perfect flux measurements. The IDEAL scenario is meant to assess the error budget and the precision floor of any other uncertainty component. The second scenario, denoted as REALISTIC, sets σ eff,std (m) = 0.2 for all measurements. The REALISTIC scenario represents a likely scenario for the uncertainty terms contained in σ eff,std (m). A specific split among its constituents is not required but is motivated by a < 10% flux measurement uncertainty, a ∼ 10% milli-lensing uncertainty, and a ∼ 15% micro-lensing uncertainty. The third scenario, denoted as EXTREME, sets σ eff,std (m) = 1, a scenario where the magnification of every single image of a glSNe is dominated by small scale micro-lensing magnification. We highlight, that these uncertainty terms should be interpreted as statistical averages for the population of glSNe. In particular, the micro-lensing component is expected to vary from image to image substantially depending on the stellar surface brightness. Table 2 summarizes our choices for the forecasts presented in this work. We emphasize that our forecast scenario and numbers operate under the assumption of imminent and complete follow-up observation after a discovery or promising candidate. The total number and numbers per year may be lower when the dedicated follow-up and we provide an extended forecast prediction as a function of glSNe in Section 4.5.

Deflector model
The model parameters for the PEMD+shear model are described in Table 3. We chose the same lens model for all glSNe systems for simplicity, but with general application of the error propagation and uncertainties. The different source position of the glSNe for the cusp, cross, and fold configurations are also provided in Table 3. We mimic high-resolution imaging data constraints on the lens model parameters with Gaussian errors on the lens model parameters, also stated in Table 3. For the image positions of the multiply lensed SNe, we assume an astrometric precision of ±0.005 arcseconds, achievable with high-resolution imaging around SNe peak brightness.  highlighted the importance and requirements on the astrometric precision of the images of the time-variable sources. Our chosen precision meets the requirement not to be the dominant uncertainty in our inference.
We sample the posterior of the imaging data (Eqn. 33) with the Gaussian likelihood in the lens model and image position parameters while demanding the image positions originating from the same source position for the proposed lens model as a solution of the lens equation. We then transform the posteriors into the relative Fermat potential and absolute magnifications at the predicted image positions (Eqn. 33).
The joint relative Fermat potential and magnification posteriors for the cusp configuration are illustrated in Figure 1. Similar posterior products are derived for the cross and fold configurations and are presented in the Appendix A (Figs. 5,  6).
The effective macro-model magnification uncertainty is ∼ 5% per image.
The effective relative Fermat potential uncertainty is ∼ 4% per image pair. The uncertainties are comparable for the three different image configurations chosen in this forecast and compatible with uncertainties obtained from the analysis of real data by the H0LiCOW/SHARP/STRIDES/TDCOSMO collaborations (Suyu et al. 2010(Suyu et al. , 2013Wong et al. 2017;Chen et al. 2019;Rusu et al. 2020;Shajib et al. 2020). The posteriors in Fermat potential and magnification for our chosen configurations and uncertainties are well approximated by multivariate Gaussians, which justifies the use the Gaussian likeihood of Equation 34 with the covariance matrix Σ ∆τ µ .

Unlensed field SNe data set
The data set of unlensed (field) SNe fulfils two purposes. First, it anchors the apparent unlensed population of SNe, m p and σ(m p ), and their uncertainties. The parameter m p directly translates to λ, and thus to H 0 . Second, the relative luminosity distances of SNe constrain the relative expansion history of the Universe, and thus Ω m in flat ΛCDM. The uncertainty on the relative expansion history can have two ways to impact the resulting H 0 uncertainty: (i) the translation of  Table 3. The configuration of the image position (diamonds), inner caustic (green) and critical curve (red) are illustrated in the top right figure. The posteriors for the cross and fold configurations are presented in Appendix A. Table 2. glSNe forecast scenarios in terms of numbers of glSNe, their redshift distribution and measurement uncertainties. The parameters of the macro model, and their uncertainties for the forecast, are presented in Table 3. The effective magnitude precision (Eqn. 51) is split between one image that is less affected by microlensing (σ eff,std ) and to the other images more strongly affected by microlensing (σ eff,ML ). the distance measurement corresponding to the glSNe systems at intermediate redshifts to the local distance constraints for a given MST parameter λ (Eqn. 48), similar to an inverse distance ladder; (ii) the translation of the apparent magnitudes from the distribution of unlsensed (mostly lower redshifts) to the glSNe source redshifts (mostly higher redshifts) (Eqn. 46). To assess current and future uncertainties coming from field SNe data sets, we set up two scenarios. First, we utilize the Pantheon data set . In particular, we are using the full covariance matrix product as described by Scolnic et al. (2018). The covariance matrix includes the intrinsic scatter in the SN Ia distribution as well as covariant systematic uncertainties. Second, we mimic a future SNe data set with an anticipated increase in the sample and lowering of systematics over the coming 10 years with the onset of the Roman Space Telescope. We use the forecast covariance matrix by Hounsell et al. (2018). The comparison between the hierarchical glSNe inference with the current Pantheon sample and the future SNe sample allows us to emphasize the importance of the field SNe sample in the next decade to utilizing glSNe to their full potential. Table 4 provides the one dimensional marginal constraints on Ω m and m p derived from the two samples. Table 4. Summary of constraints provided by the two field SNe samples used in the forecast, the Pantheon sample by Scolnic et al. (2018), and a forecast for the Roman Space Telescope by Hounsell et al. (2018).

Scenario
Ωm We perform the hierarchical analysis of the parameters and their priors presented in Table 1. We make use of the Gaussian likelihoods of individual glSNe system as presented in Section 3 with the numbers of glSNe and uncertainties presented in Tables 2 and 3. We specified three different uncertainty scenarios for σ eff,std (m) (Eqn. 51, Section 4.1.1), IDEAL (0.0), REALISTIC (0.2), and EXTREME (1.0). We also specified two different unlensed SNe scenarios, PANTHEON and FUTURE (Table 4). Any combination of SNe sample and σ eff,std (m) uncertainties results in six forecast scenarios. Figure 2 shows the posterior inference with the scenarios of the PANTHEON sample. Figure 3 shows the same inferences with the ROMAN sample.
In addition to these six inferences with a fully covariant MST component in the deflector model, we perform, as a reference for the time-delay and PEMD+shear lens model uncertainties, the forecast also without a covariant MST component by fixing λ int = 1 for both SNe scenarios. The scenarios without the MST do not depend on the error budget of the lensing magnifications σ eff (m) and the difference in the unlensed SNe sample and the glSNe sample only impacts the translation of the distance measurements into H 0 . Table 5 summarizes the results in regard of the relative precision on H 0 for the eight different scenarios considered in this work.
First, ignoring the MST, the mock data of measured time delays and Fermat potential allow one to constrain H 0 to 0.5% precision with both, PANTHEON and ROMAN sample. This set of forecast serves as a statistical reference and do not require standardizable magnifications to add information.
Once the MST is let free and only constrained by the magnifications, both the impact of the uncertainties of σ eff m and the external SNe sample significantly impact the resulting constraints. The difference between the constraining power of the PANTHEON and ROMAN sample can be seen prominently when comparing the scenarios with σ eff,ml m = 0, the IDEAL case without micro-lensing. The PANTHEON inference results in a precision of 0.8% while the increased constraining power of the ROMAN sample results in a 0.6% precision on H 0 . The error budget of the PANTHEON IDEAL scenario is dominated by uncertainties in the unlensed SNe population whereas the ROMAN IDEAL achieves almost the same precision as a scenario without an MST uncertainty.
When including REALISTIC or even EXTREME microlensing uncertainties in our forecast, the uncertainties in λ int start dominating the constraining power on H 0 as expected from the constraining power on the magnification constraints (Eqns 46, 47). Overall, we highlight our fiducial future scenario, ROMAN REALISTIC, which provides a 0.9% precision measurement on H 0 with a full 10-years LSST survey paired with a ROMAN supernovae sample.

Generalized forecast and expected timeline
Overall, the results of the full hierarchical inference performed in Section 4.4 can be well approximated with the analytical error propagation terms of Section 3.5. In this section, we make use of the analytic error propagation and generalize the forecast results of Section 4.4 for a range in the number of glSNe. Figure 4 shows the expected relative precision on H 0 as a function of the number of glSNe to be expected in the future for the three different micro-lensing scenarios and the two different external SNe samples considered in this work. In about 2 years of the LSST survey when expecting ∼ 28 glSNe, we forecast for the REALISTIC scenario a ∼ 3% precision on H 0 . With ∼ 150 glSNe for the ROMAN REALISTIC scenario we expect a 1% precision on H 0 . The precision of the external SNe sample substantially impacts the total error budget on H 0 for > 50 glSNe in the REALISTIC scenario. These numbers in terms of years of LSST survey assume an optimal follow-up effort of the discovery candidates.

DISCUSSION
The forecast results presented in Section 4 did only cover a limited range of possible systematics and opportunities regarding studying glSNe and measuring H 0 . In this section, we discuss key systematics, other windows of opportunities, and we give some general recommendations driving the design requirements in future studies of glSNe to achieve a sub 2% precision and accuracy of an H 0 measurement.

Selection effects
Brightness selection effects in the discovery and follow-up analysis of glSNe systems may pose significant limitations in the standardizable magnification methodology. Bright glSNe are easier to discover and to follow-up. Such a selection can impact unlensed brightness selection as well as local lensing magnification selection.
In our forecast and methodology, we assume an identical unlensed peak SNe brightness distribution for the unlensed field sample and for the glSNe population (m p ). Unaccounted differences between the unlensed field sample, m p,field , and the glSNe sample, m p,glSNe , results in a dif-ferential shift in H 0 by Thus, an unaccounted relative selection effect of the field SNe and glSNE of 2% results in a 1% bias in H 0 . Or in terms of an error budget, an uncertainty in the relative magnitude selection effect of 2% results in an additional error term of 1% on H 0 on top of the presented forecast results in Section 4. Local lensing magnification, a combination of micro-, milli-, and macro-lensing effects, may overall dominate the brightness selection. In particular, large (up-)scatter in brightness for rare micro-lensing events could significantly impact the selection function. It is thus crucial to understand the micro-lensing selection effect. Macro-lensing selection biases are less of an issue when performing the cosmographic analysis with time delays obtained by the identical selection function. However, when applying inferred mass profile constraints to lenses with different selection criteria, such as lensed quasars, the relative selection function comes into play.

SNe dependence with redshift and host galaxies
Beyond the glSNe systems and the required understanding of their selection function, breaking the MST and measuring H 0 also relies on an accurate and precise relative luminosity distance and intrinsic SNe distribution derived by an unlensed SNe data set. Such data sets are also used as a standalone cosmological probe or as a key component of a combined cosmological probe analysis, and their requirements and precision impact a glSNe+SNe analysis, as presented in this work.
For example, strong ∼0.1-0.2 magnitude dependence on the local host-galaxy UV surface brightness, as reported by Rigault et al. (2015), needs to be understood when making inferences from high-redshift SNe Ia. However, if there are reliable apparent magnitudes for unlensed field SNe available at the same redshifts as the glSNe, this can circumvent systematics limiting an SNe sample in measuring the late-time relative expansion history of the Universe.
We also note that with increased distance (higher redshifts) lensing effects also increasingly affect the apparent magnitudes of the field SNe sample as well. Relative selection effects (see Section 5.1.1) do also need to consider lensing selection effects in the field SNe sample.
We note that it is well known that the dust properties of SN Ia hosts, parametrised by the total-to-selective absorption ratio, R V , are very diverse and differ from the canonical value of the Milky Way of R V = 3.1 (e.g., see Brout & Scolnic 2021;Thorp et al. 2021;Johansson et al. 2021, for recent studies). Therefore, we require multiband data for each glSN  Table 1 (see also Tables 2 and 3 for details on the uncertainties) with the Pantheon unlensed SNe sample. We specified three different uncertainty scenarios for σ eff,std (m) (Eqn. 51, Section 4.1.1), IDEAL (blue; 0.0), REALISTIC (orange; 0.2), and EXTREME (violet; 1.0). Figure 3 presents the same forecast with a Roman unlensed SNe sample.  Table 5. Summary of the achieved precision on H0 for the six forecast scenarios of this work, and the two scenarios when keeping λ fixed. We specified three different uncertainty scenarios for the standardizable image σ eff,std (m) (Eqn. 51, Section 4.1.1), IDEAL (0.0), REALISTIC (0.2), and EXTREME (1.0). We also specified two different unlensed SNe scenarios, PANTHEON and FUTURE (Table 4). Any combination of SNe sample and σ eff (m) uncertainties results in six forecast scenarios. The resulting posterior inference on H0 are given in the last row. The posteriors are also presented in Figure 2 for the PANTHEON and Figure 3 for the FUTURE supernova sample, respectively.  Table 1 (see also Tables 2 and 3 for details on the uncertainties) with a Roman unlensed SNe sample (Table 4). We specified three different uncertainty scenarios for σ eff,std (m) (Eqn. 51, Section 4.1.1), IDEAL (blue; 0.0), REALISTIC (orange; 0.2), and EXTREME (violet; 1.0). Forecast Figure 2 presents the same forecast with the current Pantheon unlensed SNe sample. in our sample to constrain the R V and color excesses in the host and lens galaxies. This is important since unresolved photometry alone has been shown to underestimate the inferred magnification, as seen for iPTF16geu (Goobar et al. 2017;Dhawan et al. 2020), mandating the need for optical and NIR coverage for each image of the glSNe.

Gaussian uncertainty approximations
In the forecast of this work, we assumed Gaussian uncertainties in the measurements (linear flux units), log-normal scatter in the intrinsic SNe peak brightness distribution, as well as Gaussian scatter in the milli-and micro-lensing magnifications. The tails of the distributions need to be accurately captured to guarantee an unbiased joint inference 11 . In the current forecast, we explicitly distinguish between logarithmic and linear units and Gaussian likelihoods in either magnitude or flux units. This is not meant to be accurate for any specific scenario but primarily to emphasize the importance of accurately describing a likelihood or a posterior product. ing the probability density function (PDF) of the different components of the lensing magnifications. Specifically, non-Gaussian tails in the distributions, when combining a large set of glSNe, may significantly impact the resulting posterior PDF. The hierarchical sampling and marginalization over population distributions further poses challenges in the accuracy of the likelihood evaluation and computational requirements. Gaussian or multivariate Gaussian distributions have the advantage of analytic solutions for marginalizations and likelihood evaluations, but the assumptions of Gaussian PDF's need to be tested to the requirements of the combined posterior densities.

Opportunities
Aside from additional potential systematics considerations, there are also opportunities and circumstances that might increase the resulting precision on H 0 from glSNe relative to our fiducial forecast scenario. This section lists and briefly discusses a few of those opportunities.

glSNe without a time delay
The expected number of glSNe derived by Huber et al. (2019) that we adopt in our forecast is, in part, based on the requirement to achieve a time-delay measurement. There are potentially many more glSNe expected to be discovered (see e.g. Goldstein et al. 2019) where a precise time-delay measurement might not be expected. However, the availability of measured time delays is not the dominant source of uncertainty in our forecast. The primary information requirement to improve constraints on H 0 is foremost a precise absolute magnification measurement.

Galaxy-SNe lensing
There is also a set of "semi-strongly" lensed SNe expected with a single magnified image available that is lensing through the outskirts of a lensing galaxy. An absolute magnification measurement remains possible in the absence of multiple images and such an enhanced sample might provide significant information of the more extended galaxy density profile and thus also constraining the physically plausible MST components (see e.g., Rodney et al. 2015, for such an analysis with a singly-lensed SNe in a cluster environment). Such a probe is conceptually similar to galaxy-galaxy lensing and can possibly enhance the signal-to-noise in the very inner-most scales of galaxies where galaxy-shape information is less accessible and non-linear perturbations may arise on the distortion of the shapes (see e.g., Coupon et al. 2013, for work using magnifications of galaxies for such type of analysis).

Other type of standardizable sources
Our forecast has focused on SNe Ia, in terms of the expected numbers, intrinsic scatter and light-curve properties to measure a peak brightness and a time delay. There are other transient sources that can be standardizable. Different studies succeeded in constructing a Type II SN Hubble diagram with a dispersion of ∼ 10 − 14% in distance (e.g., Nugent et al. 2006;Poznanski et al. 2010;de Jaeger et al. 2015). The more abundant Type II SNe may provide a valuable addition. Though the light curves of Type II SNe are not as suited for time-delay measurements as with SNe Ia, there might be advantages in measuring an absolute magnification effect with Type II SNe.
Beyond SNe, there are also gravitational waves (GW) that can be standardized remarkably well and thus may open-up opportunities beyond the capabilities of SNe. Repeated fast radio bursts (FRB's) may also provide the possibility for a standardization. For GW and FRB's, one challenge will be the required astrometric precision to precisely determine the Fermat potential and macro-model magnification (see e.g., Birrer & Treu 2019).

Constraints from stellar kinematics
In our forecast, we left out anticipated constraints from stellar dynamics measurements on density profiles and breaking the MST. In part because there is a larger literature on stellar kinematics in breaking the MST and existing recent literature providing a forecast for this methodology for the decade to come (Birrer & Treu 2021). Another reason is to assess a kinematic-independent methodology in breaking the MST and thus constraints on the MST can be combined, provided both kinematics and standardizable magnifications are consistent. We highlight that stellar kinematic measurements can be performed on the glSNe lenses once the glSNe have faded away and thus might provide similar, but independent, constraining power per glSNe. Given that both methodologies are expected to provide about 1.5% precision on H 0 in the next decade, this can result in stringent consistency checks, mitigation of currently non-anticipated systematics effects and establish a precise direct distance anchor of the Universe.

Recommendations
Based on our forecast and the discussion of possible systematics and opportunities, we provide here some recommendations for the community to help guiding successful future strategies in providing both accurate and precise measurements of H 0 with glSNe. We focus on some aspects that either emerged directly from this work or deserves special emphasis.
1. Perform follow-up observations for standardizable glSNe candidates regardless of the expected time-delay precision. The precision on the mass profiles and hence H 0 relies on the ability of standardizable magnifications. Among the glSNe Ia discoveries, those systems with low expected micro-lensing events are the most valuable in breaking the MST. A significant number of glSNe Ia where at least one image is at lower magnification and lower projected stellar density are necessary, regardless of the time-delay precision (see also Foxley-Marrable et al. 2018). It is thus important to allocate significant follow-up efforts for those glSNe to be able to perform the analysis as forecasted in this work.
2. Integrate weak and strong lensing SNe analysis. To some extent, the division of the field SNe sample and the glSNe sample is an artificial cut in an underlying population of SNe that get lensed. Most lensing is weak with few percent magnification while the tails in the lensing magnification are effectively leading to glSNe. It is important to characterize the lensing effects across the entire spectrum to accurately describe the relative selection effects. With a more distant SNe sample, lensing effects may inevitably become more prominent also for the field SNe sample.
3. glSNe discovery strategy must provide a reproducible selection function. Relative selection effects are possibly a dominant source of uncertainty or unaccounted systematics. Making use of the standardizable magnification effect to break the MST, it is crucial to understand and reproduce the relative selection effect to the percent level. A survey and discovery strategy must account for the feasibility to reproduce the selection function it contains. Known selection effects can then be mitigated by e.g. large-scale simulations (see e.g., Scolnic & Kessler 2016;Kessler & Scolnic 2017, for the use for field SNe samples).

4.
Extension of the hierarchical analysis to incorporate the astrophysics of micro-lensing. The microlensing event statistics is by itself a phenomena that can probe the compact matter composition and fraction(e.g., Schechter & Wambsganss 2002;Kochanek 2004). Correlations between stellar surface brightness and (microlensing) magnification events allows one to distinguish and measure the stellar initial mass function (IMF) and other forms of compact objects, such as primordial black holes (PBH).

CONCLUSIONS
Strongly lensed supernovae (glSNe) can provide, in addition to measurable time delays, lensing magnification constraints when knowledge about the unlensed apparent brightness of the explosion is imposed. In this paper, we discussed the theoretical aspects that allow absolute lensing magnifications to constrain a key property of the lensing mass profile that is insufficiently constraint with lensing-only data due to the mass-sheet degeneracy. We then presented a hierarchical Bayesian analysis framework to combine a data set of SNe that are not strongly lensed and a data set of strongly lensed SNe with measured relative time delays. We jointly constrain (i) the unlensed apparent magnitude distribution of the population of SNe, (ii) the lens model profiles with the magnification ratio of lensed and unlensed fluxes on the population level, (iii) the relative expansion history of the Universe with the relative brightness of SNe with redshift, and (iv) H 0 with the time delays as an absolute distance indicator.
We applied our joint inference framework on a future expected data set of glSNe from 10 years of the Rubin Observatory LSST in combination with a future unlensed SNe sample from the Roman Space Telescope. We forecast that a sample of 144 glSNe with well measured time series and imaging data have the statistical power to measure H 0 to 1.0% in the next decade.
We discuss further expected covariant systematic uncertainties due to relative selection effects, dust extinction, and SNe redshift evolution. We discussed strategies to miti-gate systematics associated with using absolute flux measurements of glSNe to constrain the mass density profiles. Among the key systematic effect are relative selection biases in the discovery and usage of the glSNe and the unlensed SNe population due to micro-lensing magnification effects. We emphasize that for a 1% precision on H 0 , a 2% overall accuracy in the standardization of apparent brightness distributions between SNe population in the field and the glSNe population needs to be achieved. With an additional 1% systematic uncertainty we forecast an overall precision on H 0 of 1.5%.
The methodology presented in this work is implemented in the public software HIERARC and compatible with the hierarchical analysis by Birrer et al. (2020). The implementation allows one to combine different observational constraints self-consistently and can be adopted to the uncertain predictions of the expected glSNe depending on survey and follow-up strategies.
Using SNe is a promising and complementary alternative to using stellar kinematics observations to constrain the radial mass density profiles of strong lensing deflectors and can achieve comparable precision with independent assumptions and systematics. Future surveys, such as the Rubin and Roman observatories, will be able to discover the necessary number of glSNe, and with dedicated additional followup observations this methodology will provide precise constraints on mass density profiles and H 0 . These constraints will be key to understand the source of the current Hubble tension, and will additionally provide insights into the formation and evolution of massive elliptical galaxies.
We thank Ariel Goobar, Justin Pierel, and Sherry Suyu for useful comments on an eariler version of the manuscript. This research was supported by the U.S. Department of Energy (DOE) Office of Science Distinguished Scientist Fellow Program.

DATA AVAILABILITY
The formalism and inference schemes presented in this work are implemented in HIERARC 12 and the scripts to reproduce the presented work is publicly available 13 . Lensing calculations are performed with LENSTRONOMY 14 . In this Appendix, we provide the posteriors of the Fermat potential differences and lensing magnification for the three quad and three double configuration lensing configuration based on the lens model and source position parameters and uncertainties of Table 3. We present the quadruply lensed configurations of the cross in Figure 5, and the fold configuration in Figure 6. The cusp configuration is presented in the main body of the text in Figue 1. The posteriors for the three double configurations are provided in Figure 7.  Table 3. The configuration of the image position (diamonds), inner caustic (green) and critical curve (red) are illustrated in the top right figure.