Generalised model-independent characterisation of strong gravitational lenses V: reconstructing the lensing distance ratio by supernovae for a general Friedmann universe

TL;DR

This paper presents a model-independent method to determine the lensing distance ratio using supernovae data, enabling precise lens property measurements without assuming a specific cosmological model.

Contribution

It introduces a novel, data-driven approach to reconstruct the lensing distance ratio and lens properties in a general Friedmann universe without relying on predefined cosmological parameters.

Findings

Lensing distance ratio $D$ can be determined with 1.7% relative imprecision.

The method allows for model-independent lens characterization.

Distance measures are valid for $z\le2.3$ in any homogeneous, isotropic universe.

Abstract

We determine the cosmic expansion rate from supernovae of type Ia to set up a data-based distance measure that does not make assumptions about the constituents of the universe, i.e. about a specific parametrisation of a Friedmann cosmological model. The scale, determined by the Hubble constant $H_0$, is the only free cosmological parameter left in the gravitational lensing formalism. We investigate to which accuracy and precision the lensing distance ratio $D$ is determined from the Pantheon sample. Inserting $D$ and its uncertainty into the lensing equations for given $H_0$, esp. the time-delay equation between a pair of multiple images, allows to determine lens properties, esp. differences in the lensing potential ($\Delta \phi$), without specifying a cosmological model. We expand the luminosity distances into an analytic orthonormal basis, determine the maximum-likelihood weights for the basis functions by a globally optimal $\chi^2$-parameter estimation, and derive confidence bounds by Monte-Carlo simulations. For typical strong lensing configurations between $z=0.5$ and $z=1.0$, $\Delta \phi$ can be determined with a relative imprecision of 1.7%, assuming imprecisions of the time delay and the redshift of the lens on the order of 1%. With only a small, tolerable loss in precision, the model-independent lens characterisation developed in this paper series can be generalised by dropping the specific Friedmann model to determine $D$ in favour of a data-based distance ratio. Moreover, for any astrophysical application, the approach presented here, provides distance measures for $z\le2.3$ that are valid in any homogeneous, isotropic universe with general relativity as theory of gravity.

Tables

Table 1. Table 1: Parameter-free, analytic, orthonormal basis function sets used in this work.

Basis	Name	$u_{\alpha }\left(a\right)$	$u_0\left(a\right)$	$u_1\left(a\right)$	$u_2\left(a\right)$	$u_3\left(a\right)$
1	EdS	$a^{\alpha /2-1}$	$1/a$	$1/\sqrt{a}$	$1$	$\sqrt{a}$
2	$\Lambda$	$a^{\alpha -2}$	$1/a^2$	$1/a$	$1$	$a$
3	comp	$a^{\alpha /2-2}$	$1/a^2$	$1/\sqrt{a^3}$	$1/a$	$1/\sqrt{a}$
4	mod	$a^{\left(\alpha -1\right)/2}$	$1/\sqrt{a}$	$1$	$\sqrt{a}$	$a$

Table 2. Table 2: Λ Λ \Lambda CDM parametrisations of Planck Collaboration et al. ( 2016 ) (first row) and Scolnic et al. ( 2018 ) (second row).

$\Lambda$CDM	${\Omega }_{\mathrm{r0}}$	${\Omega }_{\mathrm{m0}}$	${\Omega }_K$	${\Omega }_{\Lambda }$	$H_0$
model					[km/s/Mpc]
(Planck)	0.0	0.3089	0.0	0.6911	67.74
(Scolnic)	0.0	0.298	0.0	0.702	73.52(1)

Table 3. Table 3: Quality assessment of the bases of Table  1 according to Section  3.4 .

B	$N_B$	${\chi }_{\nu }^2$	P>A	${\hat{c}}_0$	${\hat{c}}_1$	${\hat{c}}_2$	${\hat{c}}_3$	${\hat{c}}_4$	${\hat{c}}_5$	${\sigma }_{c_0}$	${\sigma }_{c_1}$	${\sigma }_{c_2}$	${\sigma }_{c_3}$	${\sigma }_{c_4}$	${\sigma }_{c_5}$
1	2	1.5534	✗	0.5658	-0.2402					0.0013	0.0006
1	3	0.9065	✗	0.6805	-0.3229	0.0264				0.0045	0.0032	0.0010
1	4	0.9063	✓	0.6717	-0.3158	0.0217	0.0013			0.0094	0.0083	0.0046	0.0012
1	5	0.9069	✓	0.6649	-0.3087	0.0162	0.0041	-0.0007		0.0146	0.0144	0.0102	0.0048	0.0012
1	6	0.9052	✓	0.6515	-0.2926	-0.0001	0.0169	-0.0076	0.0019	0.0163	0.0169	0.0136	0.0087	0.0042	0.0011
2	2	1.7080	✗	1.0009	-0.2296					0.0024	0.0006
2	3	0.9439	✗	0.8121	-0.1480	-0.0238				0.0069	0.0028	0.0008
2	4	0.9072	✗	0.7397	-0.1064	-0.0464	0.0061			0.0132	0.0071	0.0036	0.0010
2	5	0.9070	✓	0.7246	-0.0961	-0.0538	0.0098	-0.0010		0.0191	0.0118	0.0077	0.0035	0.0009
2	6	0.9050	✓	0.7060	-0.0810	-0.0679	0.0204	-0.0066	0.0015	0.0212	0.0141	0.0106	0.0066	0.0031	0.0008
3	2	2.9099	✗	1.1351	-0.2941					0.0028	0.0008
3	3	0.9144	✗	0.7629	-0.1150	-0.0511				0.0083	0.0039	0.0011
3	4	0.9055	✓	0.7179	-0.0858	-0.0668	0.0041			0.0161	0.0098	0.0049	0.0013
3	5	0.9064	✓	0.7168	-0.0849	-0.0675	0.0045	-0.0001		0.0209	0.0144	0.0094	0.0044	0.0011
3	6	0.9064	✓	0.7099	-0.0773	-0.0759	0.0116	-0.0041	0.0011	0.0219	0.0162	0.0126	0.0083	0.0042	0.0011
4	2	1.5841	✗	0.7229	-0.4075					0.0017	0.0010
4	3	0.9510	✗	0.9193	-0.5750	0.0419				0.0076	0.0064	0.0016
4	4	0.9062	✓	0.8163	-0.4693	-0.0050	0.0105			0.0164	0.0163	0.0068	0.0015
4	5	0.9064	✓	0.8026	-0.4536	-0.0146	0.0148	-0.0010		0.0229	0.0246	0.0133	0.0053	0.0012
4	6	0.9062	✓	0.7925	-0.4393	-0.0273	0.0242	-0.0060	0.0013	0.0244	0.0274	0.0171	0.0098	0.0046	0.0012

Table 4. Table 4: Confidence intervals of Λ Λ \Lambda CDM parametrisations of Planck Collaboration et al. ( 2016 ) (first row) and Scolnic et al. ( 2018 ) (second row) to determine confidence intervals on model-based E ​ ( a ) 𝐸 𝑎 E(a) , D L ​ ( a ) subscript 𝐷 L 𝑎 D_{\mathrm{L}}(a) , and D ​ ( z l , z s ) 𝐷 subscript 𝑧 l subscript 𝑧 s D(z_{\mathrm{l}},z_{\mathrm{s}}) .

$\Lambda$CDM	$\Delta {\Omega }_{\mathrm{r0}}$	$\Delta {\Omega }_{\mathrm{m0}}$	$\Delta {\Omega }_K$	$\Delta {\Omega }_{\Lambda }$	$\Delta H_0$
model					[km/s/Mpc]
(Planck)	0.0	0.0062	0.0	0.0062	0.46
(Scolnic)	0.0	0.022	0.0	0.020	1.62(1)

Equations

$$D(z_{\mathrm{l}},z_{\mathrm{s}})=\frac{D_{\mathrm{l}}{D}_{\mathrm{s}}}{D_{\mathrm{ls}}},$$

$${\tau }_{ij}=D\frac{(1+z_{\mathrm{l}})}{c}\left(\varphi (\mathbf{y},{\mathbf{x}}_i)-\varphi (\mathbf{y},{\mathbf{x}}_j)\right)\equiv D\frac{(1+z_{\mathrm{l}})}{c}\Delta \varphi .$$

$$\varphi (\mathbf{x},\mathbf{y})=\frac{1}{2}{\left(\mathbf{x}-\mathbf{y}\right)}^2-\psi (\mathbf{x}),$$

$$D_{\mathrm{L}}(a)$$

$$H(t{)}^2\equiv {\left(\frac{\dot{a}(t)}{a(t)}\right)}^2=\frac{8\pi G}{3}\rho -\frac{K{c}^2}{a^2}+\frac{\Lambda c^2}{3},$$

$$H(a)\equiv H_{0}E(a),$$

$$H(a)=H_0\sqrt{{\Omega }_{r0}{a}^{-4}+{\Omega }_{m0}{a}^{-3}+{\Omega }_K{a}^{-2}+{\Omega }_{\Lambda }},$$

$$D_{\mathrm{L}}(a_1,a_2)$$

$$D_{\mathrm{A}}(a_1,a_2)$$

$$f_K(r)=\{\begin{array}{cc}\frac{\sinh (\sqrt{|K|}r)}{\sqrt{|K|}}& \text{for}K<0\\ r& \text{for}K=0\\ \frac{\sin (\sqrt{K}r)}{\sqrt{K}}& \text{for}K>0\end{array}.$$

$$D_{\mathrm{L}}=10^{\frac{\mu }{5}+1}.$$

$$D_{\mathrm{L}}(a,\mathbf{c})=\sum _{\alpha =0}^{N_{\mathrm{B}}-1}{c}_{\alpha }{\varphi }_{\alpha }(a)=\mathbf{c}\circ \Phi$$

$$D_{\mathrm{L},i}=10^{\frac{{\mu }_i}{5}+1}=10^{\frac{d_i-M}{5}+1}\equiv 10^{-\frac{M}{5}}{\tilde{D}}_{\mathrm{L},i},\forall i=1,...,N_{\mathrm{D}}.$$

$$D_{\mathrm{L}}(a,\mathbf{c})=10^{-\frac{M}{5}}{\tilde{D}}_{\mathrm{L}}(a,\tilde{\mathbf{c}})=10^{-\frac{M}{5}}\tilde{\mathbf{c}}\circ \Phi .$$

$${\Sigma }_{ij}$$

$$\equiv 10^{-\frac{2M}{5}}{\tilde{\Sigma }}_{ij},\forall i,j=1,...,N_{\mathrm{D}}$$

$$k_{\Sigma }=10^{\frac{{\Sigma }_{\mu ,ij}}{5\sqrt{{\Sigma }_{\mu ,ii}}}+\frac{{\Sigma }_{\mu ,ij}}{5\sqrt{{\Sigma }_{\mu ,jj}}}}-10^{\frac{{\Sigma }_{\mu ,ij}}{5\sqrt{{\Sigma }_{\mu ,ii}}}}-10^{\frac{{\Sigma }_{\mu ,ij}}{5\sqrt{{\Sigma }_{\mu ,jj}}}}+1.$$

$$\arg \underset{\mathbf{c}}{\min }{\chi }^2,$$

$${\chi }^2={\left({\mathbf{D}}_{\mathrm{L}}-{\mathbf{D}}_{\mathrm{L}}(a,\mathbf{c})\right)}^{\top }{\Sigma }^{-1}\left({\mathbf{D}}_{\mathrm{L}}-{\mathbf{D}}_{\mathrm{L}}(a,\mathbf{c})\right),$$

$${\chi }^2$$

$$\hat{\mathbf{c}}$$

$$\underset{a}{\overset{1}{\int }}\frac{\mathrm{d}x}{x^{2}E(x)}=f_K^{-1}\left(\frac{H_0}{c}a{D}_{\mathrm{L}}(a)\right),$$

$$-\frac{1}{a^{2}E(a)}=\frac{H_0}{c}\frac{\partial f_K^{-1}\left(\frac{H_0}{c}a{D}_{\mathrm{L}}(a)\right)}{\partial \left(\frac{H_0}{c}a{D}_{\mathrm{L}}(a)\right)}\left(D_{\mathrm{L}}(a)+a\frac{\mathrm{d}{D}_{\mathrm{L}}(a)}{\mathrm{d}a}\right)$$

$$E(a)=-{\left[a^2\frac{H_0}{c}\frac{\partial f_K^{-1}\left(\frac{H_0}{c}a{D}_{\mathrm{L}}(a)\right)}{\partial \left(\frac{H_0}{c}a{D}_{\mathrm{L}}(a)\right)}\left(D_{\mathrm{L}}(a)+a\frac{\mathrm{d}{D}_{\mathrm{L}}(a)}{\mathrm{d}a}\right)\right]}^{-1}.$$

$$E(a)=\frac{\tilde{E}(a)}{\tilde{E}(a_{\max })}.$$

$$E(a)=-{\left[a^2\tilde{E}(a_{\max })\hat{\mathbf{c}}\circ \left(\Phi +a\frac{\mathrm{d}\Phi }{\mathrm{d}a}\right)\right]}^{-1}.$$

$$H(a)=H(a_{\max })E(a).$$

$$u_{\alpha }(a)=a^{\frac{\alpha }{2}-1},\alpha =0,...,N_{\mathrm{B}}-1,$$

$$D_{\mathrm{L}}(a)=\frac{c}{H_0}\frac{1}{a}\underset{a}{\overset{1}{\int }}\frac{\mathrm{d}x}{x^2\sqrt{{\Omega }_{\Lambda }}}$$

$$u_{\alpha }(a)=a^{\alpha -2},\alpha =0,...,N_{\mathrm{B}}-1$$

$$u_{\alpha }(a)=a^{\frac{\alpha }{2}-2},\alpha =0,...,N_{\mathrm{B}}-1.$$

Full text

††pubyear: 2019††pagerange: Generalised model-independent characterisation of strong gravitational lenses V: reconstructing the lensing distance ratio by supernovae for a general Friedmann universe–D

1 Introduction

1.1 Motivation from previous works

There is hardly any astrophysical research question that does not involve the distance to an object of interest. Distance measurements in our galactic neighbourhood can be performed to high precision and accuracy, Gaia Collaboration et al. (2018). Compared to that, extra-galactic distance measurements are more difficult, Cuesta et al. (2015); Muraveva et al. (2014); Tegmark (2002). Thus, with the development of an observation-based cosmic distance ladder still ongoing, cosmic distances are usually inferred from a measured (spectroscopic or photometric) redshift in combination with a cosmological model that assigns the redshift to a cosmic distance. So far, the cosmological standard model, as most precisely measured by Planck Collaboration et al. (2018), is inserted into these distance measures.

As one example, in the gravitational lensing formalism, angular diameter distances between the observer and the lens $D_{\mathrm{l}}$, the observer and the source $D_{\mathrm{s}}$, and the distance between the lens and the source $D_{\mathrm{ls}}$ appear in the lensing equations to scale the (multiple) images, the source, and the deflection potential with respect to each other. The so-called lensing distance ratio

[TABLE]

appears, for instance, in the time delay $\tau_{ij}$ between two multiple images $i$ and $j$ of the same background galaxy located at angular position $\boldsymbol{y}$ in the source plane, which is given by

[TABLE]

$\boldsymbol{x}_{i}$ are the measured angular positions of the two images in the lens plane and $\phi(\boldsymbol{y},\boldsymbol{x})$ is the lensing potential; $c$ denotes the speed of light and $z_{\mathrm{l}}$ the redshift of the lens. The lensing potential is given by

[TABLE]

in which $\psi(\boldsymbol{x})$ denotes the projected, two-dimensional gravitational deflection potential of the lens in the lens plane. A detailed introduction of the gravitational lensing formalism can be found, e.g. in Petters et al. (2001); Schneider et al. (1992).

In the previous four papers of this series, Wagner (2017); Wagner & Tessore (2018); Wagner (2018a, b), we investigated gravitational lensing from a model-independent perspective. We derived equations that determine local lens properties for different configurations of multiple images solely from observed properties of these images without assuming a specific model for the lens and determined invariance transformations of these equations to derive the degeneracies in our approach. As the next step of generalisation of this ansatz, we now replace the distance measures defined in Equations (8) and (9) based on a cosmological model by data-based distance measures that do not involve a particular parametrisation of a Friedmann-Robertson-Lemaître-Walker cosmological model. We call the derived distances therefore Friedmann-parameter-free in the following. This enables us to determine $\Delta\phi$ without the need to specify in which way the total energy density of our universe today is distributed among radiation, matter, curvature, a cosmological constant, or dark energy. Since most lenses and their background sources are located at redshifts between $z=0$ to $z=2.0$, reconstructing their distances by a supernova sample, like the Pantheon sample, Scolnic et al. (2018), is possible. Hence, we can obtain data-based distances for most lensing configurations without the need to calibrate several probes of the cosmic expansion with respect to each other.

The paper is organised as follows: in the remainder of this section, we discuss related work on the usage of supernovae to reconstruct the cosmic expansion history and highlight the differences between model fits and data-driven reconstructions of the cosmic expansion history. Then, Section 2 gives a brief introduction into cosmological model building based on Friedmann universes and discusses the observational prerequisites and basics of the standardisation of supernovae that influence the reconstructions of the cosmic expansion and distance measures. The Pantheon sample (contrary to its antecessor, the UnionSample 2.1, Suzuki et al. (2012)) is not calibrated by an overall scale. Therefore, we derive the equations to reconstruct the expansion function from a scale-free sample by normalising the expansion function. Subsequently, we insert this normalised expansion function into the definitions of distance measures. In Section 3, details about the implementation of the approach are given, before it is applied to the Pantheon sample and a Pantheon-like simulation in Section 4. Since $H_{0}$ is not determined from the scale-free supernovae data, it remains a free parameter that has to be set independently. We insert $H_{0}$ as derived from cosmic-microwave background measurements or as measured in our local neighbourhood into the distance measures. Together with the expansion function from the supernovae, we obtain data-based luminosity and angular diameter distances. We compare both choices for $H_{0}$ and show their deviations in the resulting distances. Independent of $H_{0}$, we evaluate the relative accuracy and precision of our implementation to reconstruct the expansion function and the distances up to this overall scaling by $H_{0}$. In addition, we compare our Friedmann-parameter-free reconstructions of the expansion function, the luminosity distances, and the lensing distance ratio to the ones of $\Lambda$ cold dark matter ($\Lambda$CDM) models, as parametrised by Planck Collaboration et al. (2016)111being the most recent data release at the time of performing our experiments and Scolnic et al. (2018). As a consistency check, we also perform a fit of the data from the Pantheon sample to a flat $\Lambda$CDM cosmology with our $\chi^{2}$-parameter-estimation function to show that our mathematical framework leads to a Friedmann $\Omega_{\mathrm{m}0}$-parameter value that is in agreement with the value for $\Omega_{\mathrm{m}0}$ as obtained by Scolnic et al. (2018) with a different parameter-estimation approach on the same dataset and that is in agreement with the value obtained by Planck Collaboration et al. (2016). In Section 5, we compare the precision of the lensing distance ratio to the precision of the other observables entering the time-delay equation of two multiple images of the same background source to estimate the loss of precision when dropping the parametrisation of the Friedmann model in favour of a data-based distance measure. Section 6 summarises the results, discusses the advantages and disadvantages of Friedmann-parameter-free distance measures and gives an outlook for further applications. This work focuses on the development of the Friedmann-parameter-free distance measures to determine $\Delta\phi$ for a fixed $H_{0}$. Determining $H_{0}$ for a given $\Delta\phi$ will be the subject of the next part of the paper series.

1.2 Methodology and related work of supernovae of type Ia as cosmological probes

The idea to reconstruct the evolution of the cosmic density (perturbations) from luminosity distances of standardisable candles goes back to Starobinsky (1998) and has been employed in many ways and variants since then. The basic ansatz solves

[TABLE]

for $H(a)$ having measured the luminosity distance $D_{L}(a)$ from the observer to supernovae at the left-hand side and assuming that the universe is flat. Tegmark (2002) extended and detailed the original idea and suggested combinations of different probes of the cosmic density. Huterer & Turner (1999) provided first feasibility tests for different dark energy potentials by Monte-Carlo simulations.

Depending on their way, how to invert Equation (4), several approaches are distinguished. Shafieloo et al. (2006) contains a detailed list of approaches, therefore, we focus on more recent ones here: Direct methods reconstruct $H(a)$ from the smoothed or averaged version of the original data set, e.g. as done in Shafieloo et al. (2006); Shafieloo (2007); Wang & Tegmark (2005). Basis-function methods reconstruct $H(a)$ by expanding $D_{\mathrm{L}}(a)$, the distance modulus $\mu$, $H(a)$, or a parametrisation thereof, into a set of basis functions, e.g. as done in Gómez-Valent & Amendola (2018); Ishida & de Souza (2011); Mignone & Bartelmann (2008). Bayesian methods reconstruct $H(a)$ either employing Gaussian processes or other priors, e.g. as done in Gómez-Valent & Amendola (2018); Porqueres et al. (2017); Seikel et al. (2012), and the most recently by Lemos et al. (2018) (and references therein) and Capozziello et al. (2019).

In addition, other ansatzes exist that use supernovae of type Ia to establish a model-independent hypothesis test to reject the hypothesis that the universe is not expanding, e.g. like Seikel & Schwarz (2008, 2009). Another hypothesis test can be found in Zhao et al. (2017). It is based on the Kullback-Leibler divergence and uses supernovae of type Ia in combination with other cosmological probes to find tensions between the data sets and in the current cosmological standard model.

Reconstructing $H(a)$ by fitting parametric models, e.g. a $\Lambda$CDM model, to the data, yields the most likely parameter values together with their confidence bounds. Due to the limited amount of parameters, parametric models can become unlikely for an increasing amount of data or with increasing measurement precision. For instance, the flat $\Lambda$CDM model can explain the cosmic evolution as measured by the Pantheon data set but may require extensions for future supernovae data of higher measurement precision or for increasing redshift. Extending a parametric model like the flat $\Lambda$CDM model is not unique and introducing further parameters is motivated by adding physical assumptions. As a disadvantage, the resulting model may not be the only feasible solution and may have intrinsic degeneracies among the parameters (as shown e.g. in Planck Collaboration et al. (2016)). Parametric models thus yield a description for $H(a)$ that can be directly interpreted in physical terms and that tests the compatibility of our assumptions with observations. They are the most suitable ansatz to reconstruct $H(a)$ when only sparse observational data is available.

Complementary to parametric model fitting, $H(a)$ can be reconstructed from standardisable observational data. This ansatz yields a flexible representation of $H(a)$ that is uniquely extended for an increasing amount of data. The disadvantage of data-driven reconstructions is often the data representation itself that, like in case of basis functions, lacks a direct physical interpretation. The width of the confidence bounds to which the cosmic evolution can be constrained by the data is a measure of the constraining power of these data. The wider the confidence bounds of the reconstruction, the less constraining the data are. Usually, parametric model fits have much smaller confidence bounds due to the additional assumptions they make. With an increasing amount of data becoming available, data-driven methods will overcome their disadvantage of being less precise than model-based methods and observations will replace model assumptions. Considering the supernovae of type Ia available in Suzuki et al. (2012) and in the Pantheon sample of Scolnic et al. (2018), the amount of data has almost doubled within six years and the measurement precision and the redshift range have also been increased.

In our approach to determine distance measures from a set of supernovae of type Ia, we choose a data-driven reconstruction based on statistical sampling theory as an alternative to the most commonly used parametric models. As data representation, we employ the orthonomal set of basis functions as developed in Mignone & Bartelmann (2008), which is physically motivated. This set of basis functions has been further investigated in Benitez-Herrera et al. (2012, 2013). From all distance measures mentioned above, it seems to be most suitable for the reconstruction of the lensing distance ratio $D$ in our approach to gravitational lensing, as will be further detailed in Section 2 and 3, and tested in Section 4.

2 Theoretical derivations

2.1 Cosmological model prerequisites

Assuming that our universe is spatially homogeneous and isotropic on large scales and general relativity is the theory of gravitation, the universe can be described by a Friedmann-Lemaître-Robertson-Walker metric (FLRW metric), which fulfils the first order Friedmann equation,

[TABLE]

in which $G$ denotes the gravitational constant, $\rho$ is the matter-radiation energy density, $K$ is the constant spatial curvature, and $\Lambda$ the cosmological constant. The dot denotes the derivative with respect to cosmic time $t$. Without loss of generality, we choose the scale factor today to be one, $a_{0}\equiv a(t_{0})=1$, in order to determine $a$ uniquely.

Being agnostic about the constituents and the state of the universe (i.e. not knowing anything about $\rho,K$, or $\Lambda$), we can generally define

[TABLE]

calling $H(a)$ the Hubble function of the universe with the expansion function $E(a)$ and today’s Hubble constant $H_{0}=H(a_{0})$. Thus, $E(a)$ is normalised, such that $E(a_{0})=1$. This is the most general functional form of $H(a)$ to describe the evolution of a spatially homogeneous and isotropic metric fulfilling Einstein’s field equations. Splitting the constituents into radiation, matter, curvature and a $\Lambda$-term, we arrive at the usual parametrisation of the Hubble function

[TABLE]

with $\Omega_{i}$ being today’s density divided by the critical density.

With the help of Equation (6) and the distance duality relation of Etherington (1933), we can calculate the luminosity and the angular diameter distances between two scale factors $a_{1}$ and $a_{2}$ with $a_{2}<a_{1}$, as

[TABLE]

with

[TABLE]

2.2 Observational prerequisites

In order to employ Equations (8) and (9) to reconstruct $H(a)$, observations from standardisable candles, rulers, or sirens can be used. For supernovae of type Ia, determining $D_{\mathrm{L}}$ by fitting their light curves to standardised templates has become a routine (see e.g. Amanullah et al. (2010); Betoule et al. (2014); Burns et al. (2011); Scolnic et al. (2018); Suzuki et al. (2012)). Light curve fitters, like SALT2 (Guy et al. (2010)), determine the distance modulus $\mu$, which is a function of the observed light curves (usually measured as magnitudes in the $b$-band). The light curves are parametrised by the peak magnitude, the time-stretching of the supernova, its colour at maximum brightness, the absolute magnitude $M$ for a standardised supernova of type Ia, and parameters of the light curve fitter to correct for biases, e.g. in the distance modulus due to the host-galaxy mass of the individual supernovae. The impact of these biases can be calibrated by simulations that may depend on an underlying Friedmann model with a special parametrisation. Thorough analyses have been performed that investigate the dependence of the inferred quantities, e.g. $\mu$ or cosmological parameters, on the trained and calibrated light-curve model, see Mosher et al. (2014) and Hauret et al. (2018). Enhanced simulations that employ several parametrised Friedmann models for the calibration and blind the light curve fitting with respect to a reference cosmology, have also been established, Kessler & Scolnic (2017). In the case of the Pantheon sample, the absolute magnitude $M$ of a standard supernova is left as a free parameter. Given $\mu$, we can determine the luminosity distance

[TABLE]

Setting Equation (11) equal to (8) and inserting Equation (7) into the latter, cosmological parameter values can be retrieved. Usually, nuisance parameters like the absolute magnitude $M$ are fitted together with the cosmological parameters222Determing $M$ for different, redshift-dependent subsets, as performed in Suzuki et al. (2012), revealed that all $M$ are in good agreement with each other. Thus, a single scale $M$ for the entire data set is sufficient.. $M$ and $H_{0}$ both define an overall scale for the distances and are thus not independent of each other. Since the light-curve standardisation is performed without fixing $M$, observations of supernovae only determine $E(a)$ without constraining the overall distance scale. This has to be set by one of the following options:

• •

scaling with an independently determined $H_{0}$, e.g. as derived from Planck Collaboration et al. (2016) or from Riess et al. (2018), or

• •

a measurement of $M$ for supernovae in our local neighbourhood, Richardson et al. (2014), or

• •

employing Bayesian statistics to infer cosmological parameter values after having marginalised over the prior probability of all nuisance parameters, e.g. assuming a non-informative prior probability distribution for $M$, see e.g. Conley et al. (2011) and Scolnic et al. (2018) (which is not pursued in this work).

In our work, we use the most recent Pantheon sample from Scolnic et al. (2018). It provides $\mu+M$, i.e. the distance to the supernovae up to an overall scaling factor, from a compilation of $N_{\mathrm{D}}=1048$ supernovae. These were observed between $a_{\mathrm{min}}=0.307$ ($z=2.26$) and $a_{\mathrm{max}}=0.990$ ($z=0.01$), which set the limits in $a$ (and $z$) in which our distance measure can be used. Basing on the frequentist framework of statistics for the remainder of this work, we reconstruct the normalised expansion function, $E(a)$, from these scale-free supernovae and subsequently scale with $H_{0}$ from Planck Collaboration et al. (2016). Scaling by $H_{0}$ instead of using a local value for $M$ seems to be the most consistent approach because observations of local supernovae properties might be subject to small-scale anisotropy biases that are not included in the FLRW metric from which we derive $E(a)$ and the distance measures. See e.g. Planck Collaboration et al. (2016) for a detailed discussion about the tension of currently available measurements of $H_{0}$ and Marra et al. (2013), Bolejko (2018), and Macpherson et al. (2018) for recent advances to reconcile the measurements from the cosmic microwave background and the local neighbourhood.

2.3 Scale-free series expansion of $D_{\mathrm{L}}$

In a first step, we expand the scale-free distances to the supernovae into a set of basis functions: We define the luminosity distance function $D_{\mathrm{L}}(a,\boldsymbol{c})$ as an expansion into orthonormal basis functions $\phi_{\alpha}(a)$

[TABLE]

for $a\in\left[a_{\mathrm{min}},1\right]$, in which the $c_{\alpha}$ denote the weights of the basis functions and each entry in $\boldsymbol{c}\in\mathbb{R}^{N_{\mathrm{B}}}$ is multiplied by the respective column in $\Phi\in\mathbb{R}^{N_{\mathrm{D}}\times N_{\mathrm{B}}}$ and these terms are summed up in the short-hand notation of the right hand side. This ansatz has also been pursued in Benitez-Herrera et al. (2012, 2013); Mignone & Bartelmann (2008). Since, apart from some minor restrictions detailed in Appendix A, $D_{\mathrm{L}}(a,\boldsymbol{c})$ can be exactly represented in any basis (for $N_{\mathrm{B}}\rightarrow\infty$), we do not insert a specific one until Section 3.1 and establish a quality measure to compare different bases with respect to their practical usefulness in Section 3.4.

Denoting the data as provided by Scolnic et al. (2018) as $d_{i}\equiv d(a_{i})=\mu_{i}+M$, we factor out the unknown overall scale in Equation (11) as

[TABLE]

In order to link Equations (12) and (13), we define

[TABLE]

Given the covariance matrix between the $d_{i}$ (including systematic correlations), $\Sigma_{\mu}\in\mathbb{R}^{N_{\mathrm{D}}\times N_{\mathrm{D}}}$, we obtain the entries of the scale-free covariance matrix for $D_{\mathrm{L}}$, $\tilde{\Sigma}$, by calculating

[TABLE]

with

[TABLE]

A derivation of $\Sigma_{ij}$ can be found in Appendix B.

To obtain the $\boldsymbol{c}$ in Equation (12) from Equation (13), we set up a generalised linear-least-squares parameter estimation as

[TABLE]

with

[TABLE]

in which $\boldsymbol{D}_{\mathrm{L}}=\left(D_{\mathrm{L},1},...,D_{\mathrm{L},N_{\mathrm{D}}}\right)^{\top}$ denotes the column vector of the luminosity distance measurements and $\boldsymbol{D}_{\mathrm{L}}(a,\boldsymbol{c})$ is the column vector containing the luminosity distances at the same $a$ as determined by Equation (12)333This ansatz assumes that the errors in $a$ are negligible, which is realised by incorporating the uncertainties and biases in the redshift measurements in the covariance matrix $\Sigma_{\mu}$, see Betoule et al. (2014)..

Inserting Equations (13), (14), and (16) into Equation (18), we obtain

[TABLE]

Hence, the scaled and scale-free optimisation problems, and thus, their solutions, are of the same form and yield the unbiased, consistent, efficient, and asymptotically normal generalised least-squares estimator444under the assumption that $\mathbb{E}\left[\boldsymbol{D}_{\mathrm{L}}-\boldsymbol{D}_{\mathrm{L}}(a,\boldsymbol{c})\,|\,\boldsymbol{D}_{\mathrm{L}}\right]=0$

[TABLE]

If the deviations of $\boldsymbol{D}_{\mathrm{L}}$ to $\boldsymbol{D}_{\mathrm{L}}(a,\boldsymbol{c})$ are normally distributed, it is also the maximum-likelihood estimator. For the sake of convenience, unless specified otherwise, we will drop the tilde and refer to the scale-free solution as $\hat{\boldsymbol{c}}$.

Compared to other methods that expand $E(a)$ in a set of basis functions, the ansatz pursued in Equation (12) has the advantage that it is easy to show that Equation (18) has a single global optimum which can be efficiently determined by Equation (21).

2.4 Reconstruction of the Hubble function

In the next step, we use the basis function expansion of the scale-free supernovae distances to determine $E(a)$. By construction, $E(a)$ is determined up to an overall scale. Therefore, we employ the normalisation condition $E(a=1)=1$, such that the overall scale of $H(a)$, and the consequently the scale of the distance measures in Equations (8) and (9), is given by $H_{0}$. Hence, $H_{0}$ remains the only free parameter left in Equations (8) and (9).

Luminosity distances of supernovae are measured with respect to $a=1$ today. Hence, when inserting $D_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ on the left-hand side of Equation (8), we also have to insert $a_{1}=1$ and $a_{2}=a$ on the right-hand side. If not explicitly specified, we abbreviate $D_{\mathrm{L}}(a)\equiv D_{\mathrm{L}}(1,a)$ for the luminosity distance defined by Equation (8) and analogously for the angular diameter distance measure.

Following Starobinsky (1998), we can solve Equation (8) for $E(a)$ by first isolating the integral of the right-hand side

[TABLE]

subsequently deriving both sides by $a$

[TABLE]

and then obtain $E(a)$ as

[TABLE]

We restrict the discussion to flat universes with $K=0$, so that $f_{K}(r)=r$, which is in agreement with the measurements of Planck Collaboration et al. (2016). The cases for non-vanishing curvature are analogous and treated in Appendix C for the sake of completeness.

Inserting $D_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ for $D_{\mathrm{L}}(a)$ into Equation (24), we obtain the expansion function up to an overall scale, $\tilde{E}(a)$. If the supernova sample were calibrated such that $c/H_{0}=10^{-M/5}$, the reconstructed expansion function would be normalised, $\tilde{E}(a)=E(a)$. Yet, the Pantheon sample is scale-free, i.e. not calibrated by the Hubble scale $H_{0}$, nor by an absolute magnitude $M$, so that we have to insert the data-based scale-free series expansion $\tilde{D}_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ into Equation (24). In this way, we obtain a normalised expansion function by

[TABLE]

Inserting $\tilde{D}_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ and its derivative into Equation (24), dropping all scale factors, and normalising the result, $E(a)$ is given by

[TABLE]

Consequently, when introducing the overall scale to determine $H(a)$, we use

[TABLE]

Since, for the Pantheon sample, $a_{\mathrm{max}}=0.99\approx 1$ and $H(a_{\mathrm{max}})=68.06$ km/s/Mpc for the Planck cosmological parameters (Planck Collaboration et al. (2016)) is contained in the confidence interval of $H_{0}=67.74\pm 0.46$ km/s/Mpc, we may also approximate $H(a_{\mathrm{max}})$ by $H_{0}$ in Equation (27).

Since the basis functions can be determined at all $a\in\left[a_{\mathrm{min}},1\right]$ (either analytically or numerically), Equations (26) and (27) yield a data-based expansion and Hubble function at any $a\in\left[a_{\mathrm{min}},1\right]$.

2.5 Reconstruction of distances

Inserting Equation (27) into Equations (8) or (9), we can now determine distances between arbitrary scale factors without specifying a parametrisation like in Equation (7).

If an unbiased, global measurement of an absolute distance for a standardised supernova $M$ were available, we could simply insert $M$ into Equation (14) and solve Equation (18) for $\hat{\boldsymbol{c}}$ (i.e. the scaled quantity as defined on the left-hand side of Equation (21)) to determine data-based distance measures. As we will briefly show in Section 4.4, this would lead to a higher degree of accuracy and tighter confidence bounds. Until such a measurement might become feasible, Sections 4.3 and 4.4 show that employing Equations (8) or (9) with (27) currently is the most robust way to determine data-based distances. Recent progress on the determination of $M$ can be found in Richardson et al. (2014).

2.6 Reconstruction of the lensing distance ratio

With the establishment of data-based distance measures as detailed in Section 2.5, we can determine the lensing distance ratio for all lenses and sources located within the redshift555While we use $a$ as variable for the reconstructions of the Hubble function and the distance measures, we convert $a$ to the redshift $z$ for the reconstruction of $D$ to be consistent with standard lensing notation. range of the reconstruction, as defined in Equation (1). In Sections 4.3 and 4.4, we will specify a lens redshift $z_{\mathrm{l}}$ and calculate the lensing distance ratio for all sources at $z_{\mathrm{s}}>z_{\mathrm{l}}$.

3 Implementational details

3.1 Analytic bases

3.1.1 Einstein-de-Sitter basis

As a set of orthonormal basis functions, we use the Einstein-de-Sitter basis as introduced in Mignone & Bartelmann (2008). It consists of functions

[TABLE]

that are recursively orthonormalised over the scale factor interval $\left[a_{\mathrm{min}},1\right]$ by Gram-Schmidt orthonormalisation to obtain the $\phi_{\alpha}(a)$ introduced in Equation (12). This basis is designed to recover the luminosity distances in an Einstein-de-Sitter universe ($\Omega_{m0}=1$, $\Omega_{\Lambda}=0$) with the first two basis functions (see Mignone & Bartelmann (2008) for the derivation). It is thus a sparse basis decomposition in the matter-dominated era of the cosmic evolution, for scale factors approximately ranging from $a\approx 10^{-3}$ to $a\approx 0.7$.

3.1.2 Variations of the Einstein-de-Sitter basis

In the late universe ($a\gtrsim 0.7$), the term of the cosmological constant, $\Omega_{\Lambda}$, dominates in Equation (7) and the universe expands exponentially. Determining the luminosity distances in such a cosmology

[TABLE]

yields $D_{\mathrm{L}}(a)\propto 1/a^{2}-1/a$, such that a basis obtained from

[TABLE]

is physically motivated analogously to the one obtained from Equation (28).

Since Equation (30) does not include the square-root functions contained in Equation (28), a combination of both bases can be obtained from functions of the form

[TABLE]

For the reconstruction of $H(a)$, it turns out that the term in the last bracket of Equation (24) belonging to $c_{0}$ cancels out when using the Einstein-de-Sitter basis. This implies that the expansion function is effectively reconstructed by $N_{\mathrm{B}}-1$ basis functions and coefficients, such that we also set up a modified Einstein-de-Sitter basis without the first basis function from functions of the form

[TABLE]

Table 1 summarises the functions that form the basis sets of Sections 3.1.1 and 3.1.2 after Gram-Schmidt orthonormalisation and also lists the first four functions that are employed to set up the bases. For a flat $\Lambda$CDM model in the late universe,

[TABLE]

The integral on the right-hand side cannot be simplified to a few polynomial terms as in the previous cases. Consequently, expanding the luminosity distance in such a cosmology requires an increasing amount of polynomial basis functions to achieve a highly accurate approximation. The quality of approximation for different $N_{\mathrm{B}}$ is investigated in Section 4.4.

3.2 Further numerical bases

Even sparser bases that require fewer coefficients for the reconstruction exist, e.g. the principal component basis as introduced in Ishida & de Souza (2011) or Maturi & Mignone (2009). Yet, these approaches are less suitable for our purpose than an analytic basis. The latter assumes a physically motivated decomposition into basis functions that can be efficiently determined at any arbitrary point between $\left[a_{\mathrm{min}},a_{\mathrm{max}}\right]$, while the numerical bases require to be represented by a lot of sampling points. The dense sampling increases the run-time. It also requires numerical imprecisions to be carefully taken into account and an interpolation between the sampling points to be defined (see Section 4.4 for a comparison of a numerical and the analytic implementation of the Einstein-de-Sitter basis). Apart from requiring a high accuracy and precision in the lensing distance ratio, the increasing amount of data in a data set with an increasing interval of scale factors also requires a fast reconstruction of $H(a)$ and the distance measures. This favours closed form expressions over numerical ones, especially when the confidence bounds are determined by a Monte-Carlo simulation as detailed in Section 3.3.

3.3 Confidence bounds

The uncertainties on $\hat{\boldsymbol{c}}$ due to the covariances of the $D_{\mathrm{L},i}$ are determined by a Monte-Carlo simulation from the $\chi^{2}$ in Equation (20): We simulate 1000 data sets with $N_{\mathrm{D}}$ supernovae, each at the same scale factors as are listed in the original data set by drawing $D_{\mathrm{L},i}$, $i=1,...,N_{\mathrm{D}}$ from a Gaussian distribution around the measured $D_{\mathrm{L},i}$ with a width $\sigma$ corresponding to the measured uncertainty as listed in the original data set. In this choice of error model, we assume an ideal standardisation process, such that the $D_{\mathrm{L},i}$ determined by the light-curve fits are only subject to statistical uncertainties and are not subject to any bias anymore and that the true $D_{\mathrm{L},i}$ are equal to the reconstructed $D_{\mathrm{L},i}$ up to a linear transformation. Other choices of error models are based on drawing Gaussian random samples around a true $D_{\mathrm{L},i}$ of a simulated data set. In any case, the choice of an error model heavily relies on assumptions that are very difficult to corroborate in practice.

Correlations between different data points are neglected because we focus on the imprecision that is caused by the measurement uncertainties in the $D_{\mathrm{L},i}$. The correlations that arise due to the compilation of several inhomogeneous data sets to one are usually much smaller than the statistical uncertainties. For the Pantheon sample, we find that the correlations, i.e. the off-diagonal entries, in $\Sigma_{\mu}$ are all smaller than 1% of the statistical uncertainties.

For each of the 1000 simulated data sets, $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ are reconstructed. From the ensemble of all 1000 reconstructions of each quantity, the 68%, 95% and 99% confidence intervals, corresponding to 1-$\sigma$, 2-$\sigma$, and 3-$\sigma$ confidence intervals for a Gauss distribution are calculated. In addition, we determine the standard deviation of each quantity from the 1000 simulated data sets.

While the confidence bounds on $\hat{\boldsymbol{c}}$ can also be derived from the Fisher-matrix method, the Monte-Carlo simulation is required to propagate the confidence bounds on $\hat{\boldsymbol{c}}$ accurately into $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},\mathrm{s})$ as these reconstructions contain non-linear transformations of $\hat{\boldsymbol{c}}$. By construction, the Monte-Carlo simulation yields the same confidence bounds on $\hat{\boldsymbol{c}}$ as can be derived from the Fisher-matrix, i.e. it reaches the Cramér-Rao lower bound as detailed in Appendix D.

3.4 Quality assessment of a basis

A priori, a Friedmann-parameter-free reconstruction of $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ can be performed with any set of orthonormal basis functions. To rank and compare different bases with different number of basis functions, $N_{\mathrm{B}}$, we assess their reconstruction quality by calculating the reduced $\chi^{2}$ defined as

[TABLE]

and we determine the relative imprecision of the reconstruction given by the confidence bounds determined according to Section 3.3.

The optimal reconstruction has $\chi_{\nu}^{2}=1$. If $\chi_{\nu}^{2}>1$, the basis does not fully capture the information contained in the data or the covariances have been underestimated. For $\chi_{\nu}^{2}<1$, the basis overfits the data or the covariances have been overestimated. Hence, monitoring $\chi^{2}_{\nu}(N_{\mathrm{B}})$ for an increasing number of basis functions, we can determine the basis with maximum information about $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$, i.e. which has the maximum number of $N_{\mathrm{B}}$ retrievable for a given measurement precision and which fulfils $\chi^{2}_{\nu}(N_{\mathrm{B}})\approx 1$. Bases with small confidence bounds and small $N_{\mathrm{B}}$ are preferred.

In addition, we determine the reconstruction accuracy of $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ for a basis in a simulation of a flat $\Lambda$CDM model. Since the Pantheon sample can be explained by such a cosmological model, the result is a good estimate for the reconstruction accuracy of the true, unknown cosmology of our universe. The reconstruction inaccuracies of $E(a)$ in the $\Lambda$CDM model simulation should lie within the confidence bounds of $E(a)$ caused by the measurement precision. Assuming that the inaccuracy of the data-driven reconstruction of $E(a)$ in the true cosmology is of the same order of magnitude as the inaccuracy in the $\Lambda$CDM model simulation, we do not expect the data-driven reconstruction to be biased.

To summarise, we search for the basis which comes closest to $\chi_{\nu}^{2}=1$ for the maximum number of $N_{\mathrm{B}}$ retrievable for a given measurement precision, shows the smallest bias in a $\Lambda$CDM model simulation and has the smallest confidence bounds (i.e. relative imprecisions) that encompass the relative inaccuracies of the $\Lambda$CDM model simulation at the 68% confidence level.

3.5 Run-time enhancements

Our implementation is based on MATLAB, employing the full covariance matrix including the correlations between the data points. As detailed in Section 4.2, we choose the Einstein-de-Sitter basis as orthonormal set of basis functions. Analytically performing the Gram-Schmidt orthonormalisation, we obtain closed-form expressions for $\Phi$ and subsequently closed-form expressions for $D_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$, $H(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ up to $N_{\mathrm{B}}=4$. Higher order coefficients and basis functions can be retrieved numerically. Yet, for the Pantheon sample, maximally four coefficients are significantly determined given the covariances (see Section 4.2).

A major speed-up in run-time is obtained by replacing the standard inv-function by mldivide() to calculate $\hat{\boldsymbol{c}}$ and by replacing all for-loops by matrix operations. The overall run-time to reconstruct $H(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ (the latter for one fixed $z_{\mathrm{l}}$, as detailed in Section 2.6) including confidence bounds according to Section 3.3 for the Pantheon sample amounts to roughly 110 seconds on a standard notebook (MacBook Pro, 2.2 GHz Intel Core i7, 8 GB 1333 MHz DDR3 RAM).

4 Application to data

4.1 Synopsis of the data

As observational data, we use the Pantheon data set, Scolnic et al. (2018), for the selection of the optimum basis function (Section 4.2) and to determine the reconstruction precision for $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ (Section 4.3). To investigate the reconstruction accuracies for different implementations and for varying $N_{\mathrm{B}}$ (Section 4.4), we generate Pantheon-like simulations. We simulate luminosity distances at the scale factors of the Pantheon sample from a $\Lambda$CDM model as parametrised by Planck Collaboration et al. (2016), based on measurements from the cosmic microwave background, (see first column of Table 2) and divide them by an arbitrary scale to obtain scale-free $\tilde{D}_{\mathrm{L},i}$. The parametrisation in the second row of Table 2 by Scolnic et al. (2018) is based on a fit to the Pantheon sample and measurements of $H_{0}$ in the local neighbourhood (employing $H_{0}$ from Riess et al. (2018)). We use these two parametrised $\Lambda$CDM models in Section 4.3 and in Section 4.5 for the comparison with our reconstructions.

4.2 Selection of the optimal basis

We first perform a comparison of the bases listed in Table 1 to find the optimal basis and number of basis functions, $N_{\mathrm{B}}$, given the quality assessment criteria of Section 3.4. Table 3 shows the quality measures for different configurations of the basis sets. For all bases 1 to 4, we test 2 to 6 numerically implemented basis functions. For the numerical implementation, we equidistantly sub-sample $a$ at $N_{\mathrm{S}}=10000$ points between $a_{\mathrm{min}}$ and 1. We employ the full covariance matrix in Equation (20) and determine confidence bounds as detailed in Section 3.3 from the data of the Pantheon sample.

Comparing $\hat{\boldsymbol{c}}$ for basis 1 for all $N_{\mathrm{B}}$ with the $\hat{\boldsymbol{c}}$ obtained when using the diagonal covariance matrix with the statistical uncertainties only, we find that the coefficients deviate only on the order of $10^{-4}$. Hence, correlations between the different data in the Pantheon sample play a minor role for the reconstructions of $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$, as already mentioned in Section 2.3.

For all bases, $N_{\mathrm{B}}>3$ yields $\chi^{2}_{\nu}$ slightly smaller than one, as listed in the third column of Table 3, so that we conclude that these basis configurations capture the information contained in the data well. Contrary to that $N_{\mathrm{B}}=2$ has $\chi^{2}_{\nu}>1$ for all bases, indicating that these configurations are not suitable to represent $D_{\mathrm{L}}(a)$.

Next, we determine the expected inaccuracies from the simulated Pantheon-like data set in a flat $\Lambda$CDM model simulation, as detailed in Section 4.1. We observe that the relative inaccuracies are contained within the 68% confidence bounds for configurations of $N_{\mathrm{B}}>3$ for bases 1,3, and 4 and for $N_{\mathrm{B}}>5$ for basis 2. For these configurations, we compare the standard deviations of $\hat{\boldsymbol{c}}$ with $\hat{\boldsymbol{c}}$ for all $N_{\mathrm{B}}$ to find the best basis configuration for each of the four basis sets as follows:

• •

basis 1 with $N_{\mathrm{B}}=4$,

• •

basis 2 with $N_{\mathrm{B}}=5$,

• •

basis 3 with $N_{\mathrm{B}}=4$,

• •

basis 4 with $N_{\mathrm{B}}=4$.

For basis 2, more than 6 significant coefficients could be determined. Yet, monitoring the continuously declining $\chi^{2}_{\nu}(N_{\mathrm{B}})$ for this basis and the increasing width of the confidence bounds, this basis is practically of much less use than the others, so that we consider the version of $N_{\mathrm{B}}=5$ for the sake of completeness in the following.

[FIGURE:]

Basis 1 with $N_{\mathrm{B}}=4$ turns out to have the smallest confidence bounds and thus the least relative imprecision for the reconstructions of $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ from the Pantheon sample. Figure 1 shows the relative inaccuracies of the Pantheon-like $\Lambda$CDM model simulation and relative imprecisions for the reconstructions of the Pantheon sample in form of the standard deviations of the Monte-Carlo simulation as obtained according to Section 3.3. For the plots showing the relative inaccuracies, the model-based reconstructions from $\Lambda$CDM (Planck Collaboration et al. (2016)) are subtracted from our reconstructions of the simulated Pantheon-like data set.

Thus, unless mentioned otherwise, all reconstructions in the following are determined with the standard settings for our MATLAB code as detailed in Section 3, i.e. we employ

• •

the full scale-free covariance matrix, $\tilde{\Sigma}$ from the Pantheon sample (see Equation (16)),

• •

the scale-free Einstein-de-Sitter basis of Section 3.1.1 with $N_{\mathrm{B}}=4$ basis functions (see Equation (14)),

• •

closed-form expressions for the basis functions, their derivatives to reconstruct $E(a)$, $D_{\mathrm{L}}(a)$, and $D(z_{\mathrm{l}},z_{\mathrm{s}})$ (Equations (26), (8), (1)) as detailed in Section 2,

• •

the scaling according to Equation (27) with $H(a_{\mathrm{max}})$ as determined by the $\Lambda$CDM model based on Planck Collaboration et al. (2016),

• •

1000 samples in a Monte-Carlo simulation drawn from the Pantheon data set to calculate confidence bounds as detailed in Section 3.3. We also use 1000 samples for the Pantheon-like simulations of flat $\Lambda$CDM models.

4.3 Reconstruction precision from the Pantheon sample

Running the MATLAB code on the Pantheon sample with all specifications as described at the end of Section 4.2, we obtain $E(a)$, $D_{\mathrm{L}}(a)$, and $D(0.5,z_{\mathrm{s}})\equiv D(z_{\mathrm{s}})$ as shown in Figure 2 (left column) and the relative reconstruction imprecisions (right column).

We compare our reconstruction to the respective quantities of a $\Lambda$CDM model as parametrised by Planck Collaboration et al. (2016) and by Scolnic et al. (2018) as summarised in Table 2. In order to compare $E(a)$ on equal footage, we normalise it to the value at $a_{\mathrm{max}}$ for the $\Lambda$CDM models as well.

For $E(a)$, both model-based reconstructions lie well within the confidence bounds of our reconstruction. The same applies to the model-based reconstructions for $D_{\mathrm{L}}(a)$ and $D(0.5,z_{\mathrm{s}})$ as parametrised by Planck Collaboration et al. (2016). As expected, the tension in $H_{0}$ between Planck Collaboration et al. (2016) and Riess et al. (2018) causes the $D_{\mathrm{L}}(a)$ and $D(0.5,z_{\mathrm{s}})$ as determined by the parametrisation of Scolnic et al. (2018) to lie below the 99% confidence bounds of our reconstruction because we employ $H_{0}$ as derived by Planck Collaboration et al. (2016).

4.3.1 Fitting a flat $\Lambda$CDM model as consistency check

Not finding any significant tensions with Planck Collaboration et al. (2016) and Scolnic et al. (2018), we perform an additional consistency check and employ our $\chi^{2}$-parameter-estimation function defined by Equation (18) to fit a flat $\Lambda$CDM to the luminosity distance data of the Pantheon sample. We use the average standard absolute magnitude $M=-19.25$ from local supernovae as determined in Richardson et al. (2014) to scale the $\tilde{D}_{\mathrm{L},i}$ in Equation (13). Subsequently, we insert the luminosity distance data of the Pantheon and the luminosity distance measure as given by Equation (33) and parametrised by $H_{0}$ and $\Omega_{\mathrm{m}0}$ ($\Omega_{\Lambda}=1-\Omega_{\mathrm{m}0}$) into Equation (18). The resulting non-linear least-squares optimisation problem can be solved by the standard lsqcurvefit routine in MATLAB. As $M$ and $H_{0}$ are arbitrary but dependent scales, we vary the value of $M$ to find that it only changes the fitted value of $H_{0}$ and leaves the resulting value for $\Omega_{\mathrm{m}0}$ invariant, as we expected. We obtain $\Omega_{\mathrm{m}0}=0.2870$ at $\chi^{2}_{\nu}=0.9084$. To determine confidence bounds on this value, we use the same 1000 Monte-Carlo-simulated Pantheon-like samples (see Section 3.3) that we used to determine the confidence bounds for our approach. Fitting these samples to the $\Lambda$CDM model set up above, we obtain for the median and the subsequent 68%, 95%- and 99% confidence bounds

[TABLE]

Thus, our fit agrees with Scolnic et al. (2018) within the 68% confidence bound and with Planck Collaboration et al. (2016) within the 95% confidence bound.

4.4 Reconstruction accuracy from simulated data

Having determined the reconstruction precision for the Pantheon sample, we investigate the reconstruction accuracy for a Pantheon-like simulated data set in a flat $\Lambda$CDM model cosmology, as detailed in Section 4.1. All remaining input to the reconstruction is taken from the specifications listed at the end of Section 4.2. Figure 3 shows the results. The reconstruction of $\Lambda$CDM (Planck Collaboration et al. (2016)) (or of Scolnic et al. (2018)) is subtracted from our reconstruction.

The plot on the left of Figure 3 shows the increase in accuracy for $E(a)$ with an increasing amount of basis functions. For $N_{\mathrm{B}}=3,4$, the closed-form basis functions are used, while for $N_{\mathrm{B}}=5$, we employ the numerical implementation. Using the numerical implementation of $N_{\mathrm{B}}=5$ for $N_{\mathrm{B}}=3,4$ as well, differences between the analytic and the numerical reconstructions of $E(a)$ only differ on the order of $0.001$, i.e. invisibly in this plot. While reconstructions with $N_{\mathrm{B}}=3,4$ favour the parametrisation of Scolnic et al. (2018), for $N_{\mathrm{B}}=5$ the parametrisation of the underlying simulation of Planck Collaboration et al. (2016) is finally preferred.

Reconstructing $D_{\mathrm{L}}(a)$ with $N_{\mathrm{B}}=4$, we compare the implementation with the numerical and the analytic basis functions in the central plot of Figure 3. To test the impact of the implementation on the reconstruction accuracy, we employ two different numerical implementations: one which only assumes $N_{\mathrm{S}}=1048$ sampling points for each basis function at the same scale factors as the Pantheon sample. The second one, also used in Section 4.2, assumes $N_{\mathrm{S}}=10000$ equidistantly distributed sampling points of $a\in\left[a_{\mathrm{min}},1\right[$ for each basis function.

For small scale factors, the 1% gain in accuracy between the numerical implementation sampled at the scale factors of the Pantheon sample (i.e. $\Phi$ is evaluated at $N_{\mathrm{S}}=1048$ sampling points) and the analytic basis functions can also be achieved for a numerical basis function with $N_{\mathrm{S}}=10000$ sampling points. Yet, this increases the run-time of the $D_{\mathrm{L}}$-routine by more than a factor of 4 and still causes numerical instabilities when the scale factor approaches 1.

For the lensing distance ratio $D$, determined for a typical lens at redshift $z_{\mathrm{l}}=0.5$ as a representative example, we obtain the relative inaccuracies as shown in the plot on the right-hand side of Figure 3. It shows that the reconstruction with numerical basis functions is slightly worse than the one employing analytic basis functions. In addition, we plot the relative reconstruction inaccuracies that arise when we first reconstruct $\tilde{D}_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ with Equation (14), insert the $M$ as scaling factor that the simulated data has been previously divided by, and subsequently determine $D$ from the $D_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ (using Equation (9) to convert $D_{\mathrm{L}}$ to $D_{\mathrm{A}}$ to be inserted into Equation (1)). As stated in Section 2.2, using a global measurement of $M$ as scaling would thus yield a higher accuracy for sources close to the lens (for $z_{\mathrm{l}}=0.5$ the sources should be located between $z_{\mathrm{s}}=0.5$ and 1).

4.5 Comparison to model-based reconstructions

Next, we compare the confidence intervals of our Friedmann-parameter-free reconstructions with the confidence intervals of the model-based reconstructions parametrised by Planck Collaboration et al. (2016), which is the most precise one, and by Scolnic et al. (2018), based on the Pantheon sample. The parametrisations are shown in Table 2, the confidence intervals for the parameters are summarised in Table 4.

Since the absolute scale, e.g. $H_{0}$ (see Section 2.2), cannot be determined from the supernova sample and is thus subject to the same confidence intervals for all reconstructions, we only focus on the confidence intervals of the $\Omega_{i}$. Evaluating relative uncertainties, i.e. $\sigma E(a)/E(a)$, $\sigma D_{\mathrm{L}}(a)/D_{\mathrm{L}}(a)$ and $\sigma D(z_{\mathrm{s}})/D(z_{\mathrm{s}})$, $H_{0}$ cancels out. To determine confidence bounds on $E(a)$ for the parametrisations according to Planck Collaboration et al. (2016) and Scolnic et al. (2018), we employ the definition of $E(a)$ in Equation (6) and draw 1000 Pantheon-like data sets at the scale factors of the Pantheon sample from a Monte-Carlo simulation of a flat $\Lambda$CDM model. Each Pantheon-like data set is generated with a different Friedmann parametrisation $(\Omega_{m0},\Omega_{\Lambda}=1-\Omega_{m0})$ drawn from a Gaussian distribution around the values listed in Table 2 with a standard deviation given by the confidence intervals shown in Table 4. Confidence intervals for $E(a)$ parametrised by Planck Collaboration et al. (2016) and Scolnic et al. (2018) are then derived from this Monte-Carlo simulation in the same way as for the Friedmann-parameter-free reconstruction (see Section 3.3). Subsequently, the confidence intervals on the $\Lambda$CDM-parametrised $D_{\mathrm{L}}(a)$ and $D(0.5,z_{\mathrm{s}})$ are determined. Figure 4 (left) shows the relative imprecision in the reconstructions based on the parametrisation according to Planck Collaboration et al. (2016), Figure 4 (right) shows the same plots for the reconstructions based on the parametrisation according to Scolnic et al. (2018).

We find that the imprecisions for both $\Lambda$CDM parametrisations are of the same order with tighter confidence bounds for Planck Collaboration et al. (2016) for $E(a)$ and $D_{\mathrm{L}}(a)$, which is expected from the smaller confidence bounds on $\Omega_{m0}$ (see Table 4).

As we can observe from a comparison of the plots in Figure 4 with the right-hand side of Figure 2, the precision of $E(a)$ of both model-based reconstructions is about one order of magnitude higher than the precision of $E(a)$ of the Friedmann-parameter-free reconstruction. For $D_{\mathrm{L}}(a)$, the precision of the model-based reconstructions is about a factor of three higher for the parametrisation by Planck Collaboration et al. (2016) and on equal footing for the parametrisation by Scolnic et al. (2018). For $D(0.5,z_{\mathrm{s}})$, the model-based reconstruction is about a factor of 20-50 more precise than our reconstruction.

5 Synopsis of results compared to other sources of imprecision in lensing

Summarising the results from Section 4, we find that the lensing distance ratio $D$ for a typical lens redshift $z_{\mathrm{l}}$ can be reconstructed without specialising a parametrisation for the underlying Friedmann model with a relative imprecision of the order of 1-2% (68% confidence level), 2.5-5% (95% confidence level), and 3.5-8% (99% confidence level). Parametrising the Friedmann model with the parameter values of Planck Collaboration et al. (2016) or Scolnic et al. (2018) (see Table 2), the relative imprecisions are below 0.41% for all three confidence levels.

In the model-independent approach to characterise gravitational lenses as developed in Wagner (2017), Wagner & Tessore (2018), and Wagner (2018a) so far, the lensing distance ratio only enters in the time-delay equation, Equation (2), while $D$ does not enter the equations to determine lens properties (locally constrained reduced shear and ratios of potential derivatives) by positions and shapes of multiple images in the lens plane. Therefore, we only have to consider the impact of data-based distances in Equation (2).

Having measured a time delay between two multiple images, the difference in the lensing potential between those images can be determined with Equation (2). In Wagner (2018a), we showed that time-delay measurements between multiple images fix the enclosed mass density for a given cosmological model, i.e. for a known lensing distance ratio $D$. Hence, inserting $E(a)$ as reconstructed by the supernovae and a given $H_{0}$ into Equation (2), we can uniquely determine $\Delta\phi$, given the mathematical prerequisites on $\phi$ detailed in Wagner (2018a) are fulfilled.

Assuming that the distance and time-delay measurements are performed independently, the relative uncertainty of $\Delta\phi$ is given by

[TABLE]

We conservatively estimate the imprecision of the redshift, $\sigma_{z}$, to 1%, considering spectroscopic analysis, as e.g. performed in Scodeggio et al. (2018). If the redshift of the supernova is not acquired from the supernova itself but from its host, the relative uncertainty can be of the order of 0.1%, Scolnic et al. (2018). Assuming that the time delay is determined between multiple images of quasars, the relative uncertainty, $\sigma_{\tau}$, amounts to 1-5%, see e.g. Wagner et al. (2019) for an overview on galaxy-cluster scale and Liao et al. (2015) for a systematic analysis on galaxy scale. Compared to these estimates, the relative imprecision of $D$, $\sigma_{D}$, as obtained by our reconstruction is of the same order of magnitude, while the relative imprecision of $D$ obtained by a model-based reconstruction is one order of magnitude smaller, so that we obtain

[TABLE]

where PF and PB stand for our Friedmann-parameter-free (PF) reconstruction and the parametrised Friedmann model (PB), respectively. Hence, generalising the reconstruction of the gravitational lensing potential by not specifying a parametrisation for the Friedmann model, the precision only deteriorates by 0.3%. This is the loss in precision that is caused by dropping the model assumption of individual parameters $\Omega_{i}$ in Equation (7).

Using a lens model, it is often not necessary to measure a time delay to constrain the mass density profile. Fitting the observables to the lens model already fixes the mass density profile in frequently occurring cases or additional observables as the velocity dispersions along the line of sight are employed, such that $\Delta\phi$ is determined from the lens model with its most-likely model parameters. Further details on this issue, including an analysis where $D$ enters in lens reconstructions with lens models, can be found in the follow-up paper, when we discuss different possibilities to determine $H_{0}$ from Equation (2) and their degeneracies. State-of-the-art estimates for the currently achievable accuracy and precision for $\Delta\phi$ by lens models can be found in the paper series by Suyu et al. (2017).

6 Conclusion

We investigated to which precision and accuracy it is possible to determine the lensing distance ratio $D=D_{\mathrm{l}}D_{\mathrm{s}}/D_{\mathrm{ls}}$ in a generalised Friedmann universe being agnostic about its constituents and their individual contributions. Using the latest compilation of supernovae, the Pantheon sample (Scolnic et al. (2018)), $D$ can be reconstructed for any combination of lens and source redshifts in the range of 0 to 2.3 with a relative uncertainty on the order of percent. To arrive at this result, we expanded the luminosity distance into a set of orthonormal basis functions. Due to the limited amount of supernovae and their uncertainties, relative inaccuracies in the reconstructed $D$ arise, which most probably lie within the 68% confidence bounds of the measurement uncertainties when we use an Einstein-de-Sitter basis, as introduced in Mignone & Bartelmann (2008), with 4 basis functions. Compared to three other orthonormal basis sets obtained from analytic functions, the Einstein-de-Sitter basis was shown to be the optimal basis set: It has the tightest confidence bounds that encompasses the relative inaccuracies of a flat $\Lambda$CDM model simulation and shows only a slight overfitting in the reconstruction.

As any other method to reconstruct distances based on standardised, precompiled ensembles of supernovae, our approach also relies on the methods used to fit and calibrate the light curves, which might still be dependent on a cosmological background model and which are based on additional assumptions like the one that uncertainties in the redshift estimates can be incorporated in the uncertainty of the distance modulus (see Section 2 for further details). Consequently, improvements in the standardisation process imply changes in our reconstruction method, as we do not work on the raw data directly. For instance, taking into account the redshift uncertainties separately requires our optimisation problem to reconstruct the luminosity distances to change from a least-squares formulation to a total-least-squares formulation. The determination of confidence bounds in our approach can be greatly improved, if the publicly available SNe compilations additionally provide the luminosity distances with their uncertainties.

Compared to model-based reconstructions of $D$ based on the parametrisations of the $\Lambda$CDM model by Planck Collaboration et al. (2016) and Scolnic et al. (2018), our Friedmann-parameter-free reconstruction is a factor of 20-50 less precise than the model-based reconstructions for typical lens and source redshifts between redshifts of 0.5 and 1.0. This loss can be considered as an estimate for the impact of the cosmological model on $D$, as similarly investigated in Williams & Liesenborgs (2019) for the influence of lens model assumptions on the determination of $H_{0}$ from supernova Refsdal. Propagating the uncertainties through the time-delay equation to determine the difference in the lensing potential between the two multiple images, $\Delta\phi$, we find that $\Delta\phi$ is less than 0.5% less precise for the Friedmann-parameter-free reconstruction than for a $\Lambda$CDM model by Planck Collaboration et al. (2016) or Scolnic et al. (2018). This holds for the typical redshift and time-delay uncertainties on the order of percent that are assumed to be the same for the Friedmann-parameter-free and the parametrised reconstructions.

Thus, for the model-independent characterisation of gravitational lensing configurations, we conclude that dropping the parametrisation of the Friedmann model in the lensing distance ratio $D$ in favour of a data-based reconstruction allows us to greatly generalise the method at the cost of a small and tolerable additional imprecision. As a drawback, configurations with redshifts larger than 2.3 cannot benefit from Friedmann-parameter-free lensing distance ratios yet and may require the combination of several cosmic probes.

The usage of the reconstructed $E(a)$, $D_{\mathrm{L}}(a)$, $D_{\mathrm{A}}(a)$ is not limited to applications in strong and weak gravitational lensing, distances that are determined from a data-based expansion function can also be employed in any kind of astrophysical context. Furthermore, the Friedmann-parameter-free reconstruction of $E(a)$ can be used to determine the linear growth factor, as done in Haude et al. (2019), and both, $E(a)$ and the linear growth factor are extensively used in the recently developed approach by Bartelmann et al. (2017) to set up propagators for particles moving in the phase space of an expanding universe.

A parametrised Friedmann model yields a higher reconstruction precision, yet, there is a plethora of possible parametrisations that cannot be excluded by current observations, e.g. different types of dark energy or modified gravity models, as summarised e.g. in Amendola et al. (2013) and investigated, e.g. in Benitez-Herrera et al. (2012) and Moews et al. (2019). We have already analysed the analogous situation in the previous papers of this series for the case of specifying a gravitational lens model. Hence, when considering specific parametrisations of Friedmann models, we should – analogously to specifying gravitational lens models – marginalise over all possible parametrisations to obtain the confidence bounds that represent our knowledge most realistically. Those confidence bounds should be compared with the ones obtained by the Friedmann-parameter-free method.

Acknowledgements

We would like to thank Matthias Bartelmann, Sophia Haude, Bettina Heinlein, Bruno Leibundgut, Henrik Nersisyan, Dan Scolnic, and the Galaxy Cluster Group at the Institute for Theoretical Astrophysics for helpful discussions and comments. JW gratefully acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG) WA3547/1-3. SM gratefully acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG) BA 1359/20-1.

Appendix A Validity restrictions for the series expansion

Searching for a $D_{\mathrm{L}}(a,\boldsymbol{c})$ that best fits a given set of measured data at known scale factor positions $a_{i}$, $i=1,...,N_{\mathrm{D}}$, can be treated as a sampling problem. Since the scale factor positions belong to a limited interval, the Fourier transform of the data is limited in bandwidth. Furthermore, we assume a random, not necessarily equidistant, spacing of the $a_{i}$. Then, a generalised version of the Nyquist-Shannon sampling theorem, Landau (1967), states that the function underlying the $D_{\mathrm{L},i}$ can be exactly reconstructed, if the average sampling rate is (at least) twice the occupied bandwidth of the signal.

In our case, the bandwith amounts to $2\pi(1/a_{\mathrm{min}}-1/a_{\mathrm{max}})\approx 15$, which is much smaller than the average sampling rate of $N_{\mathrm{D}}/(a_{\mathrm{max}}-a_{\mathrm{min}})\approx 1500$. Hence, for infinitely many coefficients $c_{\alpha}$ in Equation (12), the luminosity distance can be uniquely determined, while, in practice, the precision of the measured $D_{\mathrm{L},i}$ limits the number of coefficients, $N_{\mathrm{B}}$, that can be determined.

Appendix B Derivation of $\Sigma_{ij}$

We determine the uncertainty on $D_{\mathrm{L},i}$ as

[TABLE]

such that

[TABLE]

In order to identify $\Delta D_{\mathrm{L},i}\Delta D_{\mathrm{L},j}$ with the symmetric matrix $\Sigma_{ij}$, we symmetrise the right-hand side by identifying

[TABLE]

Inserting Equation (42) into Equation (41), Equation (15) is obtained.

Appendix C $E(a)$ for $K\neq 0$

For a given curvature $K\neq 0$, we perform the following transformations:

[TABLE]

Starting from Equation (8), the analogous equation to Equation (22) for $K\neq 0$ is set up in the transformed coordinates using

[TABLE]

and

[TABLE]

It reads

[TABLE]

in which we introduced the new scaling $R_{0}=c/(H_{0}|K|)$. Deriving this equation with respect to $\tilde{a}$,

[TABLE]

and solving for $E(\tilde{a})$, we arrive at

[TABLE]

Hence, employing the transformations of Equation (43), we obtain $E(\tilde{a})$ as determined by Equation (49), in which $D_{\mathrm{L}}(\tilde{a})$ is scaled by $R_{0}$, so that we can proceed with the overall normalisation as described in Section 2.4. The subsequent steps of the method remain the same as for $K=0$, taking into account the transformation of Equation (43) and $f_{\tilde{K}}$.

Appendix D Comparison of the Monte-Carlo simulation and the Fisher-matrix method

Given the expression for $\chi^{2}$ in Equation (20), the Fisher information matrix is determined by

[TABLE]

Given the unbiased estimator $\hat{\boldsymbol{c}}$ for the vector of coefficients (see Equation (21)), the Cramér-Rao bound states that the covariance of $\hat{\boldsymbol{c}}$, $\text{cov}(\hat{\boldsymbol{c}})$, is bounded from below as

[TABLE]

which means that $(\text{cov}(\hat{\boldsymbol{c}})-\mathcal{I}^{-1}(\boldsymbol{c}))$ is positive semi-definite.

Assuming we sample from a multivariate normal distribution of deviations of $\tilde{D}_{\mathrm{L},i}$ from $\tilde{D}_{\mathrm{L}}(a,\hat{\boldsymbol{c}})$ in our Monte-Carlo simulation, such that

[TABLE]

with the diagonal matrix $\tilde{\Sigma}_{\mathrm{stat}}$. On its diagonal are the variances caused by the measurement uncertainties of the $\tilde{D}_{\mathrm{L},i}$. The $\text{cov}(\hat{\boldsymbol{c}})$ based on this multivariate Gaussian $\boldsymbol{X}$ is given by

[TABLE]

with

[TABLE]

Reducing the full covariance matrix $\tilde{\Sigma}$ of the data sample to the statistical uncertainties, $\tilde{\Sigma}_{\mathrm{stat}}$, we arrive at

[TABLE]