Climbing out of the shadows: building the distance ladder with black hole images

In the era of precision cosmology it has became crucial to find new and competitive probes to estimate cosmological parameters, in an effort of finding answers to the current cosmological tensions/discrepancies. In this work, we show the possibility of using observations of Super Massive Black Hole (SMBH) shadows as an anchor for the distance ladder, substituting the sources usually exploited for such purpose, such as Cepheid variable stars. Compared to the standard approaches, the use of SMBH has the advantage of not needing to be anchored with distance calibrators outside the Hubble flow since the shadows physical size can be estimated knowing the mass of the SMBH. Furthermore, SMBH are supposed to inhabit the center of all galaxies which, in principle, means that we can measure the size of the shadows in any Supernova type Ia host galaxy. Under the assumption that the mass of the SMBH can be accurately and reliably estimated, we find that the Hubble constant can be constrained with a $\approx10\%$ precision even considering current experimental design of ground-based interferometers. By constructing a SMBH catalogue based on a specific choice of the SMBH Mass Function (BHMF), we forecast the constraints on the Hubble constant, finding that a precision of $\approx4\%$ may be within reach of future interferometers.


Introduction
In present day cosmology, one of the most puzzling results is the tension in the different available measurements of the current expansion rate of the Universe, H 0 , obtained through high and low redshift cosmological observations [1][2][3][4]. Taking the most recent local measurement H 0 = 73.04 ± 1.04 km s −1 Mpc −1 obtained by the Hubble Space Telescope (HST) and the Supernova H0 for the Equation of State (SH0ES) team [5], there is a ≈ 5σ tension with the constraint inferred by the Planck Collaboration H 0 = 67.4 ± 0.5 km s −1 Mpc −1 [6]. Notice that the latter is not a direct measurement of the expansion rate, but rather a measurement of the angular size of the sound horizon at recombination, which can be related to H 0 by assuming that the expansion of the Universe follows the cosmological constant-cold dark matter (ΛCDM) model [6].
The former result comes instead from the "distance ladder" method, which allows to infer the value of the Hubble constant from observations of Type Ia supernovae (SNIa) immersed in the Hubble flow. The idea behind the construction of the ladder is to calibrate the absolute magnitude of cosmological SNIa, by obtaining their distance through the observation of other sources in the same system for which this can be known. The local measurement quoted above is obtained using the luminosity-period relation of Cepheids [5,[7][8][9]] to obtain the distance of galaxies were SNIa are observed and, from this, their absolute magnitude, which can then be used with cosmological SNIa to obtain H 0 . Such a calibration is not the only possible, and new competitive methods have been used to anchor the distance ladder and measure the Hubble parameter, e.g. calibration using stars at the Tip of the Red Giant Branch (TRGB) from the Chicago Carnegie Hubble program (CCHP) [10,11], the Surface Brightness Fluctuations (SBF) method [12] or Mira variable stars [13,14]. The first one in particular has reached the same accuracy as the Cepheid method, but the result shows a moderate statistical tension (2σ) with the SH0ES inference [10,11]. Recent papers [15][16][17][18][19][20] have shown that such a discrepancy can be related to a difference in calibrating distances to common supernovae hosts, leading to an offset of 0.2 mag between Cepheids and TRGB SNIa calibration. Calibrations using the TRGB, independent from [10,11], show however very good agreement with both the SH0ES and CCHP results [21,22], supporting the possibility that the discrepancy in calibrating SNIa distances may be due to some internal systematics of the TRGB method, e.g. in the choice of the tip of the red giant branch in stellar catalogues [21]. In addition to these, other methods not relying on the distance ladder are also able to provide low redshift measurement of H 0 e.g., through observations of Strong Lensing Time Delay [23,24], Gravitational Waves [25,26], Quasars and Gamma-Ray Bursts [27][28][29][30], Cosmic Chronometers [31][32][33] and Fast Radio Bursts [34].
While no evidence that systematics can account for the tension in the H 0 measurements was found [5][6][7][8][35][36][37][38], it is helpful to find new ways of calibrating SNIa measurements, in order to minimise the risk of falling under systematic effects.
In this paper, we investigate the possibility of applying the distance ladder method with a completely new observable to anchor the cosmological SNIa. This relies on the possibility to use the shadow of a Super Massive Black Hole (SMBH) as a standard ruler, due to the relation existing between the size of the shadow and the mass of the SMBH [39][40][41][42]. Such a possibility has been opened by the recent observations of Event Horizon Telescope (EHT), a telescope array consisting of a wide network of radio telescope, that for the first time has observed the immediate environment of a SMBH [43][44][45].
The possible use of SMBH shadows as standard rulers was already explored in the literature [46,47], with the aim of constraining cosmological parameters. While our approach starts from the same idea, here we focus on the possibility of using such observations as a new calibrator for the distance ladder, thus as a viable option to obtain an estimate of H 0 with SNIa observations independently of the cosmological model.
The aim of this paper is therefore to determine the constraints achievable on the Hubble parameter by current and future interferometers dedicated to horizon-scale observations of SMBH. We also review the impact of the experimental uncertainties in the measurements of the size of the shadows and of the mass of the SMBH. We further assess how the number of observed SMBH will be affected by the angular resolution threshold of an ideal ground-based interferometer.
The paper is organised as follows. In section 2 we review the main steps of the distance ladder method to obtain H 0 from SNIa measurements, while in section 3 we explore the possibility of using measurement of the angular size of SMBH shadows as a standard ruler to calibrate the SNIa data set. We produce synthetic data, as expected to be observed from the present and next generation of ground-based interferometer, following the approach of section 4 and obtain the forecasted precision on H 0 achievable with this method in section 5. We finally summarise our conclusions in section 6.

The distance ladder
The distance to any astronomical object can be written in terms of the so-called distance modulus where m(z) and M are the apparent and absolute magnitude of the source while d L (z) its luminosity distance; The luminosity distance can be related to the line elements of the Universe given a form of the metric. Choosing a Friedmann-Lemaitre-Robertson-Walker (FLRW) line element i.e. ds 2 = −dt 2 + a 2 (t)dr 2 (assuming c = 1), the distance luminosity can be written in terms of the comoving distance χ as where we have introduced the Hubble rate, (1 + z)H(z) = −dz/dt in terms of the redshift z = a −1 (t) − 1. At low redshifts, Equation 1 can also be rewritten as where the quantity a B is the intercept of the magnitude-Hubble relation and it is defined, for a generic expansion rate at z > 0, with a cosmographic expansion [48] e a B +0.2m 0 = z 1 + 1 2 where m 0 = m(z ∼ 0), q 0 is the deceleration parameter and j 0 is the jerk. Thus, if for a given astrophysical object one has measurements of its apparent magnitude, and a way to obtain its absolute magnitude, it is possible to exploit these relations to obtain a measurement of H 0 .
Indeed, calibrating sources to obtain such a measurement is not always feasible and therefore the sources commonly used for such an approach are "standard candles", i.e. sources whose absolute magnitude is constant. The most common source used to infer cosmological distances are SNIa, sources that are extremely bright, and thus can be seen at quite high redshifts, and that behave as standardisable candles, i.e. there exists a relationship between SNIa peak luminosity and the shape of their light curve. While the peak luminosities are not exactly identical for all SNIa, due to the differences in mass and chemical composition of their progenitors, their intrinsic dispersion is small and therefore their absolute magnitude M can be taken as constant (see e.g [49][50][51][52][53] and reference therein for a detailed discussion about SNIa light curve reconstruction). Therefore, if one is able to calibrate the absolute magnitudes of SNIa at very low redshifts where their distance can be measured also by other means, such a calibration can be assumed to be valid also for SNIa that are immersed in the Hubble flow, and can therefore be used to infer H 0 .
One method to achieve such a measurement is to identify SNIa in galaxies where Cepheids stars are present. Thanks to their period-luminosity relation, which also needs to be calibrated through observations of Cepheids in nearby hosts, one is able to obtain the distance to the host galaxy, and therefore to use this to calibrate the magnitude M of the SNIa [5,[7][8][9].
Connecting the distance of nearby Cepheids, with the ones in SNIa hosts and then with the SNIa in the Hubble flow is what is typically called the distance ladder and each step of the ladder is called rung. Practically, the distance ladder relies on a three rungs approach: 1. The period-luminosity relation of Cepheids is calibrated using sources that are close enough for their distance to be measured via parallax. This allows to obtain the luminosity of more distant Cepheids; 2. The distance of SNIa whose host galaxies also contains a Cepheid can be obtain through the period-luminosity relation of the latter. With the distance to the SNIa measured, one can use Equation 1 to obtain a measurement of M ; 3. with such measurement of M , assuming the same absolute magnitude holds for all SNIa, Equation 3 can be used with the data of SNIa in the Hubble flow to obtain H 0 .
As we noted in section 1, this is not the only possible approach to calibrate the SNIa absolute magnitude. In particular, the second rung of the ladder can be performed by substituting Cepheids stars with any other methodology able to measure distance to galactic hosts: SBF, Mira variable and the TRGB are perfect examples of distance ladders independent of Cepheids distances. These alternative and independent methodologies have been proven successful in the inference of H 0 , showing consistency within the total error budget across the various methods. However the terrific increase in accuracy that these methodologies have experienced in the latest decades have shown that some small discrepancies exist between them [5,10,14]. We further stress that the common denominator to all these methods is the use of observations that need to be calibrated with nearby distances (or anchored as it is commonly said in literature). This constitutes the first step of the ladder and the same anchors are used across the various methods. While still within the region of being moderate statistical fluctuations, these discrepancies have led to question not only the validity of these approaches, but also the assumptions made in all of the rungs of the distance ladder [15][16][17][18][19][20][21]54].
It is therefore timely to find new methodologies that can overcome some or all these problems, in particular the need for anchoring in local hosts.
In this work, we investigate a different method to calibrate supernovae observations, and therefore to obtain the current rate of the expansion of the Universe H 0 , exploring the possibility of obtaining the distance to low redshift supernovae by observing the shadow of SMBHs in their host galaxies. As we elaborate later in this manuscript, SMBH shadows can be self-calibrated (if one has measurements of the mass of the BH casting the shadow) removing the need for anchoring with nearby hosts. Furthermore, SMBH (and BH generally speaking) are supposed to be present in any elliptical or spiral galaxies (and perhaps less massive and active BH with mass within 10 4 and 10 6 solar masses may be present in irregular and dwarf galaxies) making them very abundant in the Universe.

Black Hole shadows as standard rulers
The recent observations by EHT have highlighted how the measurement of the immediate environment of SMBHs is now within the reach of current facilities [44,[55][56][57][58]. Such observations can provide significant information since, when surrounded by an emission region, black holes are expected to exhibit a dark shadow caused by effects due to gravitational lensing of light rays and by the capture of photons at the event horizon [59,60]. The shadow is caused by the difference in luminosity between the accelerated particles falling into the BH and the BH itself (which is a perfect absorbing surface or else a perfect black body at temperature T ∼ 0 1 ). The apparent angular size of the shadow can be obtained from the size of the cone of photon escaping [39] where we assumed a Schwarzschild BH, defined the reduced mass as m = GM BH /c 2 and r is the distance between the observer and the centre of the BH. For a distant observer one can take the limit r 2m and promote the radial coordinate of the Schwarzschild metric into the angular diameter distance d A of the FRLW metric (see e.g. [40,41] for a full derivation). This leads to (assuming sin θ ≈ θ) where l BH = 3 √ 3m is the physical size of the shadow. BH shadows can be therefore used as standard rulers if one is able to obtain l BH , i.e. a measure of the BH 1 It is worth noting that the BH surface is expected to emit a thermal radiation by Hawking mechanism [61][62][63]; however such a radiation would have an energy so small that it would be negligible with respect to the luminosity of the photons orbiting around the BH. mass, M BH . This allows to infer the distance of a BH by measuring, independently, the apparent angular size of the shadow for a system at redshift z, θ BH (z), and the mass M BH of the SMBH.
It is important to stress here that until this point we assumed that the BHs under examination can be described by the Schwarzschild metric. However, such an assumption cannot hold for all BHs, since when the spin of the BH is non-vanishing, one needs to use metrics able to describe rotating BHs, such as the Kerr metric [64]. Throughout the rest of this paper, we keep relying on the Schwarzschild description; we assess the impact of considering Kerr BHs on our approach in Appendix B.
In addition to the correct description of the spacetime around the BHs, several other complications might affect the modelling that leads to Equation 5 [65]; evolving accretion mechanisms and changes in the plasma distribution of the accretion disk can indeed lead to modifications of this equation [66,67], and second order effects on the photons travelling to the observer, such as weak gravitational lensing, can also have an impact. In this work, we neglect these complications, taking advantage of the fact that we are only using BH shadows as calibrators for the distance ladder. This allows us to focus the analysis on systems at low redshift, where the details of the BH are better understood and modelled, while the effect of weak gravitational lensing is negligible, due to its nature of an integrated effect along the line of sight.
The recent observations by EHT highlighted how the relation between θ BH , l BH and M BH can exploited to measure with great precision the mass of the BH, having measurements of its distance and of the apparent angular size of the shadow [45,68].
In this work, we take a complementary approach, and investigate instead the possible use of such observations for cosmological inference provided one is able to measure the mass of the SMBH with other methods. We show that, under these assumptions, the observation of shadows of SMBHs is able to provide a measurement of the distance between these systems and the observer.
Knowing the black hole mass and the angular size of the shadow on the sky, it is indeed possible to infer the angular distance of the BH in a similar fashion to what is done for Baryon Acoustic Oscillations (BAO), which, similarly, rely on the observation of a standard ruler, i.e. the size of the sound horizon at recombination, to infer angular distances [69,70].
Compared to BAO, however, the BH shadows have the advantage of not being related to recombination physics, but only rely on the ability to measure the mass of the BH independently from the measurement of the size of the shadow , and on the validity of the assumed theory of gravity (General Relativity in this work). There is yet another advantage, the BAO peak in the galaxy 2D correlation function has amplitude that decreases as z → 0 so that in our local Universe the peak is suppressed by non-linear fluctuations produced by peculiar velocities or other bulk effects making difficult to measure the BAO angular scale below z 0.05 [69,70] . Conversely, the BH shadows angular size grows linearly with the BH mass which, in turn, is expected to be maximal today since the BHs have had time to grow accreting mass from the stars and gas that surround them or by merging with other BHs [71,72]. This peculiarity of the BH evolution allows us to measure bigger shadow sizes in our local Universe than at cosmological distance and open the possibility of using the BH shadows to perform distance measurements at low redshifts.
In the following section we describe in detail how measurements of d A (z) could be achieved by upcoming observations; however, in order to use such measurements as a calibrator for the distance ladder, we still need to convert the angular distance of the shadow (d BH A ) into a luminosity distance (d BH L ). Throughout the paper, we will assume that the mutual scaling of the two distances is fixed by the reciprocity theorem d A (1 + z) 2 = d L ; such a relation implicitly assumes that the Distance Duality Relation (DDR) holds, an assumption that could be violated in non-standard cosmological models [73][74][75][76]. However, departures from this relation are a cumulative effect, and in the low redshift regime, where we perform our investigation, the DDR holds at least at first approximation. This allows to write the luminosity distance to the BH shadows as:

Supernovae and Black Holes data sets
In order to forecast the results of the distance ladder calibrated with upcoming measurements of SMBH shadows, we need to simulate the data sets that will be available in the near future. We need therefore to create data for both SNIa and observations of SMBH shadows. We take • mock catalogues for future observations of SMBH shadows, where we impose a cut in the angular size of the shadow corresponding to the angular resolution of EHT and post-EHT experiments; • mock SNIa data for the Legacy Survey of Space and Time (LSST) of the Vera C. Rubin observatory [77].
We provide details on how these data sets are built in subsection 4.1 and subsection 4.2 respectively.

A catalogue for SMBH shadows
In this section we describe how to create our synthetic realisation of a SMBH shadows catalogue compatible with upcoming astrophysical observations. We start from the local SMBH distribution to get a catalogue of events with an associated mass and redshift. Given the redshift and assuming a cosmological model, we can obtain the angular diameter distance for each of these events, while from Equation 5 and Equation 6 we can obtain both the intrinsic and apparent size of the shadow for each of the events. We then cut out of the catalogue the events for which the shadow cannot be resolved with the chosen experimental set up, thus introducing a minimal observable size θ cut , and we associate an error on the observables (M BH and θ BH ) for each of the events.

The SMBH mass function
The first ingredient we need for our catalogue is the distribution in redshift and mass of the SMBHs that can be observed. This can be obtained starting from the local BH mass function (BHMF), which is the number of SMBH per unit volume and mass Obtaining this distribution from observations is not trivial, as astronomical surveys are inevitably incomplete. This is mostly due to selection bias in the counts of SMBH in any mass bin and because there is yet no consensus on a broadly applicable and precise technique to determine accurately the mass of the SMBH [78]. Thus many different variants of the BHMF exists nowadays. [79][80][81][82][83][84][85] As commonly done in literature [86][87][88], we choose here a Schecter-like BHMF with functional shape where φ is a normalisation constant and M is the cut-off mass of the BH population. Tuning the parameters of the BHMF we can represent very different BH populations and assess the impact of different parameterizations on the final population of SMBH. We choose here to use a representation of the BHMF as derived in [89,90]. The parameters of such BHMF are derived phenomenologically fitting observations, obtained inferring the SMBH masses from the bolometric luminosity of early and late-type galaxies with the functional form above. The fit is performed in the log-range of BH masses log 10 M BH ∈ [6,9] in the local Universe (z ∼ 0). This BHMF can be obtained fixing the parameters to α = −1.19 , log 10 M = 8.4, and log 10 φ = −3 [87]. Consequently, this BHMF is characterised by an exponential cut-off for masses above M ∼ 10 9 M . However it has been argued that ultra-massive BH with masses above 10 9 do exists in nature, but a cutoff may still exist preventing SMBH to have arbitrary high masses [72,83]. To account for this, we use a different value for M that we fix corresponding to the estimates of [91] to M = 3.5 × 10 10 M . Such a value is obtained calibrating the estimates in [83] for the upper limit of the mass of local SMBHs with the observational constraints on the mass of the SMBH in NGC600 [92]. This pushes the exponential cutoff of both mass functions to masses around 10 11 M . We will refer to this model as SH09MX.
Until now, we have discussed the local BHMF; however, the SMBHs we observe today have evolved from primordial seeds by (primarily) accreting matter from their surroundings and by merging. While the knowledge of accretion mechanisms have evolved significantly in the past decades [88,[93][94][95], a well motivativated physical model that can describe the evolution of SMBH throughout cosmic history is far from being conceived and one has typically to rely on phenomenological models to describe the mass evolution of SMBH (see e.g [96,97]).
It is therefore clear that the BHMF needs to evolve in time, thus in order to asses the number of SMBH of a given mass at different cosmological epochs one needs to model such an evolution. We choose here to follow a fully phenomenological approach that is based on few simple assumptions to construct a "visibility function", reproducing the evolution of the BHMF in a simplistic way. While we do not expect it to be as accurate as a model which takes into account the complicated mechanisms involved in SMBH ac-cretion, we will show that it gives realistic prediction given the uncertainties in the modelling of the SMBH mass function.
The assumptions that we use to construct our visibility function are the following: • the conditional probability to observe a SMBH with mass M BH at redshift z, P (M BH |z), is described by a log-Gaussian distribution, and, therefore, the distribution in log M will be Gaussian; • the higher mass tail of the BHMF has evolved from seeds with mass of the order 10 5 − 10 6 M , thus we assume this is the lower limit of the mass for the BH we consider; • we assume that the scaling in redshift of the mean and variance of the BHMF follows the scaling of distances in FLRW, Implementing these three conditions we can describe the redshift evolution of the BHMF by combining the local mass function with the conditional probability at a given redshift.
The first assumption allows us to write the conditional probability of measuring a mass M BH at redshift zero as a Gaussian distribution in x = log 10 M BH . The variance of this distribution is fixed requiring that at 2 standard deviation the maximum mass measurable is 10 11 M , in agreement with what we have discussed previously about the existence of a cutoff mass for ultra-massive BHs. This consideration can be translated in an equation to determine the standard deviation of the Gaussian σ M by imposing [75] This assumption would still allow for the presence of ultra-massive BH with M BH > 10 11 M with a probability of 5%. In other words, the high mass tail of the distribution would populate our catalogue with ultramassive BH with masses that can exceed significantly the threshold of 10 11 M . However, it must be taken into account that when combining with the BHMF, this would weight the higher mass tail of the distribution through an exponential cutoff which prevents SMBHs to have masses much bigger than M . Applying the second condition means that at high redshift a considerable number of SMBH still exist with mass of the order of 10 5 −10 6 M . This translate into an higher probability of measuring these masses compared to masses higher than 10 6 M . Combining this with the third condition, implies that only the variance of the distribution of P (x) will scale with redshift while its mean remains constant.
The final distribution for the mass M BH in logspace will be constructed normalising the peak of the distribution to unity at all redshift. The form of this distribution is as follows: This avoids that at z ∼ 0 the distribution deviates significantly from the local BHMF and at higher redshift the distribution would only shrink towards small masses rather than changing also in amplitude. In other words as we increase the value of z the distribution will scale self-consistently. Notice that, in order to exactly recover the local BHMF at z = 0, one would need to assume a uniform distribution i.e. P (x|z = 0) = 1, which would return the local BHMF when the two are combined. However, we find that our assumption of a Gaussian distribution leads to negligible deviations as it impacts only the higher mass tail of the BHMF. Multiplying this distribution with the local BHMF of Equation 9, φ BH , we obtain our final expression for the BHMF at given mass and redshift It is worth stressing that for the purposes of this work, using SMBH shadows to calibrate low redshift SNIa, we only need a catalogue of local SMBH, which are described by φ(x). However, the redshift distribution Φ(x, z) is necessary to compute the total number of systems, which will directly impact the number of systems that fall in the redshift range that we will use to calibrate SNIa. Indeed, a different choice of P (x|z) could impact the results and affect the final outcome of our methodology. While we rely on this simple approach throughout the rest of this work, we also assess the impact of such an assumption in Appendix A.

Simulated observations
With the redshift dependent BHMF now defined, we can generate our synthetic SMBHs catalogue, i.e. a set of mock observations of mass and redshfit. To do so, we apply the same procedure used to realize synthetic catalogues of gravitational waves observations (see e.g. [98,99]) and integrate Φ across the chosen range of mass and redshift. By construction, the integral of Φ across a range of mass and redshift gives the number of available SMBH in the comoving volume: where A is a normalisation constant and we have assumed a standard FLRW background to express the comoving volume in terms of the comoving distance χ(z) and the Hubble parameter H(z).
From the above equation, it is clear that we can easily normalise the BHMF to unity and use it as a probability density to extract the mass and redshift needed for our SMBH catalogue once the total number of SMBH observables in the Universe is fixed. The estimates of [91] show that the number of SMBHs exceeding the observational threshold of shadows with size θ > 10 µas and flux F > 10 −3 Jy and having optically thin disks (i.e. for experiment like EHT) is O(10). We found that fixing N tot = 10 6 our catalogue has only 6 shadows that exceed the threshold of 10 µas in agreement with the estimates of [91]. We therefore fix the total number of observables BHs (conservatively) to N tot = 10 6 . This will be our baseline setting for the total number of SMBHs. We stress however that the estimate of [91] falls short in predicting the number of M87-like shadows of a factor of ∼ 50 possibly making our estimates over pessimistic; we will assess in section 5 how changing N tot affects the final results.
Having fixed the normalisation through Equation 13 we can finally write the bidimensional distribution for mass and redshift as and the masses and redshifts of the SMBHs in our catalogue can then be obtained performing a MCMC sampling of P (x, z). Following this approach, we end up with a set of N pairs of {M BH , z BH }; we then assume a fiducial ΛCDM cosmology, with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.29, which allows to associate a distance d BH A to each redshift. Using Equation 6 we can then obtain the apparent angular size of the shadow for each mock SMBH. The catalogue is therefore transformed in a set of observations {M BH , θ BH } for each event drawn. Of the full catalogue we only keep the events at low redshift (z 0.01), i.e. the range in which we will calibrate the SNIa. Finally, to each event left we associate observational errors that we assume as σ θ /θ BH = 0.07 [43], and σ M /M BH = 0.07 [100][101][102]. This assumption on the relative uncertainties of the observables can be extremely optimistic, in particular for what concerns the measurement of the BH mass. While for the rest of our analysis we assume the values above, we also assess the impact of less optimistic assumptions on the final constraints in Appendix C.
Notice that the choice of the starting local BHMF can significantly impact the distribution of masses in the catalogues and determines how many of the simulated shadows will be actually observable.
We show the final baseline catalogue in Figure 1, and we recap the assumptions made to produce it in Table 1 It is worth noting that, the large range of masses for the SMBH population makes the angular size of the shadow span several orders of magnitude even at similar redshifts. However the state-of-the-art of BH shadows observations [43,91] only allows to observe angular sizes above a certain threshold size θ cut . For the EHT instrument θ cut = 10 µas [43,91] leaving most of the SMBH available in the comoving volume invisible. In order to account for this sensitivity threshold, we apply a hard cut to our simulated catalogue, only preserving events for which θ BH > θ cut . We will however assess in section 5 how changing this baseline simulation EHT Figure 1: The SMBH shadows baseline catalogue obtained from the SH09MX mass function. We report here with yellow data points only the SMBH appearing below z ≤ 0.01 and for θcut = 10µas. The red data point shows the event observed by EHT.
cut can affect the results.
The events shown in Figure 1 show those that survive both the angular and redshift cut (in red), alongside the data point from available EHT observations (yellow).

SNIa catalogue
As discussed above, a catalogue of SMBH shadow measurements will not be enough to constrain H 0 using the distance ladder. We also need SNIa observations that can be calibrated with SMBH measurements, together with a set of SNIa in the Hubble flow. The making of this additional data set requires again to estimate the redshift distribution in order to produce a mock catalogue. We start here from estimating the number distribution of SNIa from the Pantheon catalogue (which contains 1048 of such observations) [49]. This is done by reconstructing the redshift distribution of SNIa from the sources of the Pantheon catalogue and then applying an inverse transform sampling technique 2 to obtain a statistically equivalent version of the Pantheon data set in redshift space. 2 Inverse transform sampling is a technique that allows to translate a uniform distribution into a generic one using its cumulative distribution function (CDF). This can be done solving a integral equation of the form : We employ the sampled redshifts to create the set of mock observations needed for the inference of H 0 through the distance ladder. The values of the relative magnitudes at each redshift is then calculated using Equation 1 and assuming the same fiducial cosmology employed for the SMBH catalogues.
For this work, we fix the number of resampled SNIa to be 1000 in the range z ∈ [0.01, 10] and the value of the intercept a B = 0.715 corresponding to its ΛCDM value [7]. This range of redshift is used to ensure that the distribution of the SNIa in the data set is consistently reconstructed by the inverse sampling algorithm. However, the high-redshift tail of the distribution does not significantly exceed z ∼ 2 and basically no SNIa are found above z ∼ 3 as expected from current observations. As a further check, we note that out of 1000 resampled SNIa we found ∼ 270 below z ≤ 0.15 consistent with the real Pantheon data set [49]. This ensures that our synthetic SNIa data set has exactly the same statistical properties of the observed one.
To each SNIa relative magnitude we associate an error due to the brightness uncertainties following [103] where σ flux = 0.01 is the systematic uncertainties related to flux calibration , the intrinsic scatter of SNe at fixed colour is σ scat = 0.025, σ intr = 0.12 is the intrinsic distance scatter and we include an intrinsic dispersion in the distance modulus of the form: δµ(z s ) = e M z s with e M drawn from a Gaussian distribution N (0, 0.01) [103]. We use this error for all the SNIa in the second and third rungs of the ladder.

Constraints on the Hubble constant
In this section we provide details on our analysis pipeline, which exploits the method of section 2 using the observables of section 4. We then provide the constraint one can obtain on H 0 with both current and forecast data.

Analysis method
We start noting that the distance ladder performed with SMBH shadows is only composed of two rungs, precisely the second and third (see section 2), as SMBH shadows can provide calibrated distance measurements if the mass of the BH is known. To perform the second rung, a pair observations of a black hole shadow and of a SNIa in the same galactic host is required. If such configurations exist, one can estimate M from each pair as M = m SNe − 5 log 10 d BH L = 5 log 10 H 0 − 5a B , (16) where the above relation implicitly assumes d BH L = d SNe L . For each SMBH shadow in our catalogue, we assume that a measurement of a SNIa is available in the same galaxy, and we translate the SMBH distance, d BH L , into an estimate of M . We then combine all the inferred M into a joint likelihood that we employ to constrain cosmological parameters.
The second step of our ladder consists in fixing a B . This can be obtained fitting the low redshift part of the cosmological SNIa catalogue that we discussed in subsection 4.2 with the cosmographic expansion of Equation 4.
The full pipeline combines these two rungs in a simultaneous fit of the five free parameters {M, H 0 , a B , q 0 , j 0 }, where the absolute magnitude M is used to fit the magnitude measurements inferred from the SMBH catalogue, while the remaining four enter the expressions for the luminosity distance, Equation 1 and Equation 3, which are used to fit the low redshift SNIa. We employ a Gaussian χ 2 likelihood for a B and we combine this with the joint likelihood for M to obtain our final constraints on the five parameters. These parameters are sampled using the public sampler Cobaya [104], which employs a Metropolis-Hastings algorithms to sample the parameter space [105,106], assuming flat priors on the parameters set. The final results are obtained analysing the samples using GetDist [107].

Constraints using current and forecast data
We now apply the analysis method we presented to the baseline catalogues for SMBH and SNIa generated following subsection 4.1 and subsection 4.2. This will provide us with forecast results for upcoming observations of multiple SMBH images and estimate the precision that one will be able to achieve on measurements of the Hubble constant. With the baseline settings of Table 1, we find that our catalogue contains six observables shadows for z 0.01, i.e. observations available to anchor the distance ladder.
In addition to this, we also assess the precision that could be achieved with the currently available observations, i.e. a single shadow observed with the sensitivity of EHT. In order to do so, we use the EHT  Table 2: Mean values and 68% errors on the free parameter of our analysis, obtained using the EHT-like and baseline data sets ( Table 1).
The EHT collaboration used such a measurement alongside an independent estimate of the distance, to obtain a measurement of the SMBH mass. Here instead we assume that an independent estimate of the mass is available, with a measurement precision of 7%. Furthermore, we simulate the presence of an observed SNIa in the same host of the SMBH, whose magnitude is obtained from our fiducial cosmology and with an uncertainty following subsection 4.2. We name such a calibration data set, containing a single shadow and a single SNIa, EHT-like. We show the results obtained with EHT-like and baseline data sets in Figure 2, and we report the constraints on the free parameters in Table 2. From these results we can see that with observations equivalent to those currently available, one can already obtain a measurement on the Hubble constant, with H 0 constrained with a precision of ≈ 10%. In the baseline forecast case, with only a handful more images available, the precision on H 0 improves to 4%, i.e. a bound larger than those obtained with common methods to anchor the ladder, but still tight enough to be used to obtain independent measurements of H 0 .
Our results show how the bounds on a B and q 0 do not change switching from EHT-like to the baseline data set. This is due to the fact that the constraining power on a B and q 0 is mostly determined by the low 3 Let us note that in [43] it is reported a measure of the size of the observed diameter of the shadow of the SMBH in M87. Therefore in the main text we reported this result and the corresponding uncertainty divided by two. However one can instead multiply by two both side of Equation 5 to obtain the angular size of the shadow diameter. Note also that there is a small asymmetry in the photon ring due to the SMBH rotation as discussed in [108] redshift end of the SNIa in the Hubble flow, and they are therefore independent from the uncertainty of the calibrator distances. The uncertainty on H 0 is instead determined by the measurement error achievable on the absolute magnitude M , which is directly related to the errors on the observables, θ BH and M BH . In Table 2 and Figure 2, we can also notice how the jerk parameter j 0 is not constrained by the analysis. This is due to the fact that j 0 enters the cosmographic expansion at the second order in redshift; this means that its impact will be relevant at higher redshift with respect to q 0 . Since the redshift range of the cosmological SNIa we use for our analysis is limited to z < 0.15 (see subsection 4.2), we do not reach an high enough redshift for j 0 to be relevant and therefore we are not able to constrain it.
As mentioned in subsection 4.1, in our baseline settings, we assumed that the observational error on both the angular size of the BH (σ θ BH ) and on the its mass (σ M BH ) are 7% of the measured value for these observables. These assumptions on σ M BH and σ θ BH are quite optimistic given that the observation of SMBH shadows is still far from becoming a mature field for cosmological inference [65]. We therefore assess the impact of less optimistic assumptions on the observational errors in Appendix C, finding that this can significantly impact the final results of our approach.
Finally, we want to stress that a constant percentage error over a whole data set is highly unrealistic; each measurement of a SMBH shadow is independent from the others and the peculiarities of each system can make shadows very different across the catalogue. The examples proposed in this manuscript would therefore corresponds to a situation in which systematic uncertainties of the SMBH are subdominant compared to the experimental errors.

Angular threshold and total number of BHs
While we saw above that changing the errors on the observables with respect to our baseline choices can significantly degrade our constraints, we have two other assumptions that can instead impact the results: the angular size threshold for detection (θ cut ) and the total number of expected BHs (N tot ). Both these parameters will change (with respect to our baseline) the number of SMBH shadows that can be employed for the distance ladder, but in a slightly different way. The value of N tot will determine the available SMBH per bins of mass and redshift, i.e. the number of observable shadows. The value of θ cut will instead determine how many of those would be actually measured by an experiment with given angular resolution. We show the bounds obtained on H 0 changing these parameters in Figure 3, with the left panel referring to changes in θ cut and the right panel showing the changes due to N tot . As expected, an increase in θ cut leads to a lower fraction of the existing SMBH being observables (f obs ), which translate into looser bounds on H 0 . For the range of θ cut explored here, a significant change in the precision of this method is found, with our extreme cases being σ H 0 /H 0 (θ cut = 1 µas) ≈ 2% and σ H 0 /H 0 (θ cut = 20 µas) ≈ 9%. Similar results are found when changing N tot , which directly affect the number of BH available at the redshifts of interest. We find in this case that the precision on H 0 can vary between σ H 0 /H 0 (N tot = 10 5 ) ≈ 7% and σ H 0 /H 0 (N tot = 10 7 ) ≈ 2%.

Conclusion
One of the most pressing issue in modern cosmology is that of understanding the nature of cosmological tensions. A promising avenue is represented by the possibility that these discrepancies may be related to the failure of the standard cosmological model, thus unveiling the opportunity of discovering new physics beyond it. In this manuscript we have proposed a new way of performing the distance ladder with cosmological SNIa using horizon scale observations of SMBH as distance calibrator. Compared with the usage of Cepheid variable stars, SMBH shadows have the advantage of having absolute calibration (the physical size of the shadow) fixed by the theory of General Relativity (see also section 3).
This feature of SMBH shadows allows us to perform the distance ladder with only two rungs; in other words the BH shadows do not require being anchored with nearby calibrator like in the case of Cepheid stars. Furthermore SMBH are thought to inhabit any galaxy in the Universe thus, provided that the shadow and the SMBH mass can be measured independently, any SNIa measurement would have an associated SMBH calibrator.
As a proof of concept for the proposed pipeline, we applied it to mock data sets of SMBH and SNIa. First, see section 4, we generated a synthetic realisation of a SMBH shadow catalogue bounded on a specific choice of the BHMF that we produced using a Schecter-like functional form, also employing a toy-model for the conditional probability of finding a SMBH with a given mass at a given redshift. We then imposed a cut in redshift (z BH ) to keep only the SMBH present in the redshift range used to calibrate the distance ladder, and a cut θ BH > θ cut = 10 µas to remove from the catalogues all those SMBH shadows whose angular size is below the observation threshold. Second, we generated a mock data set for SNIa, ensuring that this is statistically equivalent to the observed Pantheon catalogue. With these two data sets in hand, we applied the distance ladder methodology we described in section 2 and section 3 to forecast the precision it could achieve on a measurement of H 0 . We found that our approach could provide a precision of ≈ 4% on H 0 with the catalogues obtained with our baseline settings. We also quantified the impact of other assumptions we made during the production of the synthetic data set, finding how more pessimistic assumptions on the errors on θ BH and M BH , or more stringent requirement on the threshold size θ cut , could worsen the expected precision.
Overall, our results seem to indicate that the approach presented in this work could be used as an alternative method to measure H 0 , and that it can therefore contribute to the investigation of the Hubble rate tension by providing an independent measure of this parameter. Nevertheless, the precision on H 0 that we estimated may be in the reach of upcoming observations if a reliable method to infer the mass of the SMBH is found. Currently, as discussed in section 3, there is no agreement in how to determine reliably the SMBH mass, which is critical to determine the distance of SMBHs directly from their resolved shadows.
The results presented in this work rely on the assumption that our choice of the metric describing the gravitational field around the SMBHs holds. Such an assumption is necessary to obtain the relation between the size of the shadow, the mass of the BH and its distance from the observer. Our baseline results are obtained assuming a Schwarzschild metric, and we investigated the effects on the results of a different assumption by switching to a Kerr metric (Appendix B); however, several other possibilities are available (see e.g. [42,109]) and, potentially, observations of BH shadows can be used to obtain constraints on the possible metrics describing the SMBH environment, an investigation that the EHT collaboration has started, finding consistency with the Kerr metric [110][111][112].
In concluding, while the observation of SMBH shadows is still in its infancy, and many systematic effects need to be better understood in order to obtain robust and reliable constraints, the future of these measurements looks very promising and may lead to a golden era for the study of the very foundations of gravity and of our Universe.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. In constructing our BHMF and its scaling with redshift we have assumed to know the form of the conditional probability P (x|z) appearing in Equation 12. However, as we noted already in the main text, this assumption is quite simplistic as it does not take into account all the complicated processes involved in SMBH accretion physics (see e.g [87,91]). We therefore review in this appendix the impact of this assumption over the behaviour of our catalogue for our baseline model.

Acknowledgements
Let us start discussing the effect of our assumption on the evolution of the BHMF by first looking at the behaviour of the marginalised distributions of mass and redshift of the BHMF, which affect the features of the SMBH in our catalogues. These functions can be obtained simply integrating over either the mass M BH or the redshift z, i.e.
For the conditional probability we choose instead a model with an additional parameter describing the index of the scaling with redshift i.e.
and we choose two values of p corresponding to our baseline model i.e. p = −1 and a value corresponding to p 1 = −0.8.
The corresponding marginalised distributions obtained with the two values of p are shown in Figure A.4. As we can appreciate from the plots, the change in the value of p significantly impact the behaviour of P (x) for SMBH masses 10 9 but it has a negligible impact for masses above this threshold. The reason is that the ultra massive BH in the tail of the distribution are exponentially less likely than those with smaller masses due to the behaviour of the BHMF we have chosen and this sums up to the imposed scaling of P (x|z) leading to a fast depauperation of the high mass tail of P (x) as the redshift increases (see also Figure A.5). However it must be taken into account that P (x) is integrated over redshift and, therefore, the difference between the two curves accumulates as the redshift increases. Looking at the full 2D distribution ( Figure A.5), it is clear that at low redshift the scaling of the conditional distribution has little impact, as its variance only decreases linearly with z. Thus, for the redshift of interest for the construction of the distance ladder, the choice of P (x|z) is negligible even though it has a significant impact on the full 2D distribution.
Furthermore, the impact on the redshift distribution is also negligible for the redshifts that are interesting for the distance ladder (i.e. z ∼ 0). Here the reason is that P (z) depends on the size of the comoving volume and not only on the scaling of P (x|z). For very small redshifts, P (x|z) is almost a constant in z, therefore for z ∼ 0 the scaling is determined solely by the scaling of the comoving volume, i.e. dV /dz ∼ χ 2 (z)/H(z). Therefore, there is no appreciable effect varying p on the SMBH redshift distribution at low redshift.
In conclusion the distribution of SMBH at low redshift in our catalogue is determined by the chosen BHMF and the assumed cosmological model. We note however that choosing an accurate scaling for P (x|z) and an accurate form of the BHMF is mandatory if one want to use SMBH shadows to reconstruct the distance-redshift relation at arbitrary redshift and not only to make the distance ladder.

Appendix B. Spinning SMBHs and their shadows
The analysis performed in this paper relies on the assumption that the observed SMBHs have a vanishing spin and can be described by a Schwarzschild metric, leading to the relation for the shadow size of Equation 5.
However SMBHs do rotate and the effect of their spin on the surrounding environment alters the shape and the size of their shadows, quantities that now need to be computed in the Kerr metric, describing rotating BHs [113][114][115]. This is the case e.g. for the two SMBHs for which the shadow sizes have been measured by the EHT collaboration [45,68,108,112].
It is therefore important to assess the impact of rotation on the constraints obtained in this work. Generally speaking, both the BH spin and the observer inclination angle will affect the size of the shadow, with departures from Equation 5 being maximal for an observer on the equatorial plane looking at a nearlymaximally spinning BH. In such a case, the shadow will be deformed of ∼ 7% along the horizontal axis while preserving the same size of a Schwarzschild shadow on the vertical axis (these axis are defined in the plane perpendicular to the line of sight) [116][117][118].
Therefore, rotating BHs will exhibit an asymmetric shadow, with a degree of asymmetry determined by the value of the BH spin and observer inclination [42]. Such an asymmetry, while complicating the description of the shadow size in terms of the BH parameters, allows to constrain the BH spin and observer inclination angle at the same time, independently of the mass of the BH [119,120].
In principle, one can use the asymmetry between the vertical and horizontal diameter of the shadow to improve constraints on the distance measure (see e.g. [42] and reference therein). However, while such quantities are well defined theoretically, on the observational side it is very complicated to infer these two lengths, as the centroid of the shadow cannot be determined accurately from observations. In order to overcome this issue, one can make use of quantities averaged over the full extension of the shadow, which are well connected to what EHT-like interferometers would measure. One example is the average radius of the shadow, defined as [117,118] R ≡ 1 2π We note that analytical solutions for the shadow size in a Kerr space-time can be found only for a handful of cases [42], with none available for a generic spin and inclinations. Therefore, one needs to rely on phenomenological relations, obtained from numerical simulations of light rays travelling around a Kerr BH to determine the average size of the shadow. With this approach, the average radius (for an observer at distance much bigger than the horizon size) can be obtained in the form [118]: with R 0 and R 1 being where a = J/(cM BH ), J being the angular momentum of the BH, and θ i being the angle between the observer line-of-sight and J.
To understand the effect of spin and inclination in the determination of the Hubble constant, we introduce an additional step in the analysis proposed in the main text. Instead of fitting directly to obtain H 0 from distances obtained with Equation 5, we now need to first determine the distance comparing Equation B.2 to the measured size of the shadow varying the mass, the spin, the inclination and the distance as free parameters; we can then convert the distance obtained this way into a constraint on H 0 . This procedure must be done singularly for each SMBH in our catalogue, as it is done e.g. to determine the luminosity distance from gravitational wave events [25]. For the purpose of showing the impact of spin and inclination on the constraints on H 0 we apply this exercise to the EHT-like case we proposed in the main text.
The results are presented in Figure B.6, where we show three limit cases for illustration: 1. a general situation where only the mass of the BH is known while the spin and inclination are free to vary. We choose broad priors for these parameters namely a ∈ [−1, 1] and We label this case Kerr BH (purple contours); 2. same as the previous case, but assuming the BH is not spinning (a = 0), which corresponds to our baseline case where the shadow size is determined using Equation 6. We do however still use the inclination as a free parameter. We label this case SW BH (yellow contours); 3. an optimal case where BH mass, BH spin and observer inclination are known (red contours). This case is exactly like the Kerr BH one, but with the inclusion of a prior on spin and inclination as determined in [120] from the photon dynamic of the shadow of M87. We label this case as Kerr BH + prior.
In all cases we assume the mass of the SMBH in M87 to be M BH = (6.6 ± 0.4) · 10 9 M [100].
To quantify the shift due to the inclusion of spin and inclination of the BH, we quantify the percentage error on the distance with the formula: whereD BH is the mean value of the posterior of the BH distance and σ +,− the right and left limit at 68% C.L. respectively . We find that σ rel ∼ 0.133 for case (1) while σ rel ∼ 0.10 for case (2) and (3). We further note that the spin and inclination parameters have a weak correlation with the determination of the BH distance when these parameters are free to vary. However, this correlation becomes more prominent for high spin, |a| > 0.5. This can be manifestly seen in the small shift of the distance posterior in the Kerr BH + prior case, see Figure B.7. The uncertainty on the distance however stays unmodified compared with the SW BH case. To determine the impact on the value of the Hubble constant we include this additional scatter in the EHT-like mock proposed in the main text. We found that the value of the Hubble constant is minimally affected leading to H 0 = 69.8 +6. 6 −8.1 km s −1 Mpc −1 , an increase of 1% of the uncertainty.
We further note that due to the (small) correlation between spin, distance and mass of the BH, an increase in the uncertainty of the mass determination would reduce the impact of spin and inclination on the final errors, as the uncertainty on the mass dominates. As an example, a factor two worsening of the uncertainty on M BH would lead to a 2% error on the distance instead of the 3.5% we find with our baseline M BH uncertainty. On the contrary, should the mass measurement improve significantly, the inclusion of spin and inclination in the analysis would become crucial.
However given the current observational uncertainties, the impact of the BH spin and observer inclination has almost no influence on the measurement of distances from SMBH shadows. These parameters, in fact, increase the uncertainty on the distance measure only by ∼ 3.5% ( 1% on H 0 ) if their values are completely unknown. Nevertheless, given that those parameters can be determined by the dynamic of the accretion disk around the BH [120,121], with some prior knowledge of the spin and inclination pa- rameters it is possible to achieve the same accuracy obtained with Equation 6 on D BH . Overall, our analysis hints to the fact rotating BHs do not pose significant limitation to the use of SMBH shadows as cosmological probe, at least with the current level of observational uncertainty. The dominant contribution to this measurement is still given by the determination of the mass (as well as the shadow size) even in the case of Kerr BHs.

Appendix C. On the mass determination of cosmological SMBHs
Throughout the analysis done in this paper, even when varying the observational specifications such as the value of θ cut , we have kept fixed the uncertainties on the measurements of θ BH and M BH to the level of 7%. This was motivated by existing observations (see section 4), but such uncertainties are indeed optimistic when dealing with observations of shadows at distances higher than those currently observed.
In particular, there are no reliable and robust methods that can be applied to every system to infer the mass of SMBH, a significant (and arguably the most prominent) limitation to the use of these objects as distance tracers. The main issue is that the avail- able methods are precise, but not necessarily accurate, an indication that astrophysical systematics are still poorly understood. As an example, M87 has two mass determinations to date: stellar kinematics gives M BH = (6.6 ± 0.4) · 10 9 M [100] while gas dynamics gives M BH = (3.5 +0. 9 −0.7 ) · 10 9 M [102]. The two determinations have an accuracy of ∼ 7% and ∼ 22% respectively, but they exhibit a ∼ 3.5σ tension showing that systematic uncertainties are still far from being under control. The EHT collaboration has also determined the mass from the shadow size of M87 to be M BH = 6.5 ± 0.2| stat ± 0.7| sys · 10 9 M [45] showing good agreement with the measurement of [100] and a 3σ discrepancy with that of [102].
In order to avoid this issue, one could rely on more conservative estimates of the observed quantities (θ BH and M BH ), in order to avoid inaccuracies due to systematic effects. This however increases the observational uncertainties on the distance to the BH, and will impact the final constraints on H 0 . In order to estimate how the final constraints degrade with increased uncertainties, we performed again our baseline analysis, but this time we vary, separately, the relative uncertainties σ θ /θ BH and σ M /M BH in the range [7%, 80%].
Our results are shown in Figure C.8, yellow and red lines showing how changes in σ θ /θ BH and σ M /M BH , respectively, translate into changes in the relative error on H 0 . We find that indeed the precision with which the observables are measured have a significant impact, with the final constraint on H 0 ranging from ≈ 5% to ≈ 40% moving from the most optimistic to the most pessimistic case.