Model Independent Expansion History from Supernovae: Cosmology versus Systematics

TL;DR

This study uses a model-independent approach to analyze supernova data, revealing apparent deviations from standard cosmology at high redshifts that can be mitigated with simple corrections, though statistical significance remains uncertain.

Contribution

It demonstrates that apparent high-redshift deviations from $ ext{Lambda}$CDM in supernova data can be explained by simple corrections, emphasizing the importance of model-independent analysis.

Findings

Deviations at $z extgreater1$ can be reduced with a linear correction.

Model-independent methods do not find statistically significant deviations.

Corrected data remains consistent with $ ext{Lambda}$CDM.

Abstract

We examine the Pantheon supernovae distance data compilation in a model independent analysis to test the validity of cosmic history reconstructions beyond the concordance $\Lambda$CDM cosmology. Strong deviations are allowed by the data at $z\gtrsim1$ in the reconstructed Hubble parameter, $Om$ diagnostic, and dark energy equation of state. We explore three interpretations: 1) possibility of the true cosmology being far from $\Lambda$CDM, 2) supernovae property evolution, and 3) survey selection effects. The strong (and theoretically problematic) deviations at $z\gtrsim1$ vanish and good consistency with $\Lambda$CDM is found with a simple Malmquist-like linear correction. The adjusted data is robust against the model independent iterative smoothing reconstruction. However, we caution that while by eye the original deviation from $\Lambda$CDM is striking, $\chi^2$ tests do not show the extra linear correction parameter is statistically significant, and a model-independent Gaussian Process regression does not find significant evidence for the need for correction at high-redshifts.

Tables

Table 1. Table 1 : Δ ​ χ 2 Δ superscript 𝜒 2 \Delta\chi^{2} and Ω m subscript Ω m \Omega_{\text{m}}  after correction for a potential magnitude bias due to the host galaxy metallicity evolution of Eq. ( 9 ).

	$\Delta {\chi }^2$	${\Omega }_{\text{m}}$
$\Lambda$CDM	0	0.298
$\Delta M_B\left(z\right)$ correction	-0.040	0.293
$\Delta M_B\left(z\right)$ correction+0.5	-0.109	0.282
$\Delta M_B\left(z\right)$ correction+1	-0.176	0.251
Correction (13)	-0.902	0.293

Table 2. Table 2 : Crossing statistic mode 0 ( T 0 subscript 𝑇 0 T_{0} ) and mode N ( χ 2 superscript 𝜒 2 \chi^{2} ), for the 23 data points with z > z c 𝑧 subscript 𝑧 c z>z_{\text{c}} , for different corrected (c) and uncorrected (uc) parametric ( Λ Λ \Lambda CDM) and non-parametric (smoothing) “theory" values, and data. Italics indicates that the fitting is run with the opposite set (corrected vs uncorrected) as the data, and should be used only as a crosscheck.

	Data-uc		Data-c
	$T_0$	${\chi }^2$	$T_0$	${\chi }^2$
$\Lambda$CDM-uc	6.94	16.35	3.47	15.40
$\Lambda$CDM-c	11.33	16.48	1.27	15.35
smoothing-uc	1.06	15.60	12.03	16.15
smoothing-c	25.12	16.95	0.26	15.45

Equations

$$\mathcal{D}(z)$$

$$h^2(z)$$

$$+(1-{\Omega }_{\text{m}}-{\Omega }_k)\exp \left(3{\int }_0^z\frac{1+w(x)}{1+x}\mathrm{d}x\right)$$

$${\hat{\mu }}_{n+1}(z)$$

$${\mathbf{W}}_i(z)$$

$${\delta \mu }_n{|}_i$$

$$Om(z)=\frac{h^2(z)-1}{(1+z{)}^3-1}.$$

$$w(z)=-1+\frac{1}{3}\frac{d\ln [h^2-{\Omega }_{\text{m}}(1+z{)}^3]}{d\ln (1+z)}.$$

$$\Delta M_{\mathrm{B}}$$

$$Z_b(z)$$

$${\rho }_{\ast }(z)$$

$${\rho }_{0\text{c}}$$

$$\psi (z)$$

$$\mu \to \mu +\delta \mu ,$$

$$\delta \mu =\{\begin{array}{ll}A(z-z_c),& \mathrm{if}z>z_c\\ 0,& \mathrm{otherwise},\end{array}$$

$$T_0$$

$$\left[\begin{array}{c}\mathbf{y}\\ {\mathbf{f}}_{\ast }\end{array}\right]$$

$$k(\mathbf{u},\mathbf{v};{\sigma }_f,\ell )$$

$$\log p(\mathbf{y}|X)$$

Full text

††pubyear: 2018††pagerange: Model Independent Expansion History from Supernovae: Cosmology vs Systematics–A

1 Introduction

Type Ia supernovae (SNIa) distance indicators have proved to be one of the most successful probes of the cosmic expansion history, leading to the discovery of its acceleration (Perlmutter et al., 1999; Riess et al., 1998). The currently most complete SNIa distance compilation, Pantheon (Scolnic et al., 2018), has more than a thousand spectroscopically confirmed SNIa. Upcoming surveys such as Euclid (Laureijs et al., 2011), Large Synoptic Survey Telescope (LSST, Ivezic et al., 2008), or Wide-Field Infrared Survey Telescope (WFIRST, Spergel et al., 2015) will increase the number of potential SNIa candidates and extend the data to higher redshifts. While the statistical error will improve greatly with the next generation data, it is important to improve our understanding of systematic uncertainties accordingly.

This is especially true when seeking to verify or falsify the current concordance cosmology of $\Lambda$CDM, dominated by a cosmological constant and cold dark matter. In order to detect robustly small variations in cosmological model, we must control the astrophysical and observational aspects of the data. Use and interpretation of SNIa survey data must account for systematics such as the Malmquist bias, population evolution, gravitational lensing, and dust extinction. These have been reasonably well characterized at $z\lesssim 1$ (Conley et al., 2011; Scolnic et al., 2014; Rubin et al., 2015; Brout et al., 2018; Hinton et al., 2018), but are more uncertain at higher redshift. Some are complicated functions of the survey characteristics, for example the Malmquist bias is a well-known selection effect where at redshifts near the survey magnitude threshold, brighter supernovae are more likely to be detected, biasing the effective luminosity towards brighter values, therefore inferring shorter distances.

Cosmologically, SNIa constrain the shape of the distance-redshift relation and hence the Hubble expansion rate that depends on the energy density contents of the universe. For a Friedmann-Lemaître-Robertson-Walker (FLRW) metric, the dimensionless comoving distance $\mathcal{D}$ is defined for any curvature parameter $\Omega_{k}$ by

[TABLE]

is the Hubble expansion history (squared), $\Omega_{\text{m}}$ and $\Omega_{k}$ are the matter and curvature density parameters at $z=0$, and $w(z)$ is the equation of state of dark energy.

Given a specific cosmological model, i.e. $\Omega_{\text{m}}$, $\Omega_{k}$, $w(z)$, one can assess the goodness of fit to the data. Certainly one can test the concordance flat $\Lambda$CDM model with $\Omega_{k}=0$, $w=-1$. However, it is also of interest to explore the cosmological data in a model independent manner, without claiming to know the energy density constituents and their behaviors (especially given our ignorance of dark energy properties). That is, we use Eq. (1) in terms of the Hubble parameter $h(z)$, without using its specific construction in terms of components in Eq. (2). For clarity, and motivated by theoretical considerations from inflation and observational considerations from the cosmic microwave background (CMB), we do adopt spatial flatness $\Omega_{k}=0$, so $\mathcal{D}=\int dz/h(z)$. We will call this the model independent approach, meaning it does not assume a dark energy model, though it does assume FLRW and flatness.

Shafieloo et al. (2018a) reconstructed the expansion history from the Pantheon supernova compilation in such a model independent manner. The reconstructions show validity of rapidly growing expansion history at high redshifts, corresponding to a flattening of the luminosity distances (shorter distances than expected from $\Lambda$CDM) or negative distance moduli relative to $\Lambda$CDM. While in that paper the focus was on consistency with baryon acoustic oscillation (BAO) distances and the CMB sound horizon, exploring $\mathcal{D}$ vs $h$ to test the FLRW and flatness assumptions, and comparing growth vs expansion to test general relativity, here we focus on the expansion history and comparison with concordance $\Lambda$CDM cosmology, and exploring implications of any deviations for cosmology, supernova properties, or survey systematics (see also Ferramacho et al., 2009; Tutusaus et al., 2018).

In Sec. 2 we briefly review the implications of the data on model independent reconstruction of cosmological quantities. We then examine three avenues of explanation for deviations from $\Lambda$CDM: cosmology (Sec. 3), source properties (Sec. 4), and survey characteristics (Sec. 5). We conclude in Sec. 6.

2 High Redshift Deviation

The behaviour of the smooth reconstructions in Shafieloo et al. (2018a) shows that a substantial deviation from $\Lambda$CDM is allowed by the data at high redshift.

We applied the iterative smoothing algorithm (Shafieloo et al., 2006; Shafieloo, 2007; L’Huillier & Shafieloo, 2017), taking into account the covariance matrix of the Pantheon sample (Shafieloo et al., 2018a). Starting from the best-fit $\Lambda$CDM model, we stop the procedure after 200 iterations and kept all reconstructions yielding a better $\chi^{2}$ than the best-fit $\Lambda$CDM model, which correspond to a subsample of the reconstructions in Shafieloo et al. (2018a). For the sake of completeness, we briefly explain the method here. Starting from some initial guess for $\hat{\mu}_{0}(z)$, we iteratively obtain the reconstructed $\hat{\mu}_{n+1}(z)$ at iteration $n+1$ as follow:

[TABLE]

and $\textbf{{C}}_{\text{SN}}$ is the covariance matrix of the Pantheon data. In case of uncorrelated data ($C_{ij}=\delta_{ij}\sigma_{i}^{2})$), we recover the formula used in Shafieloo (2007) and L’Huillier & Shafieloo (2017). We thus end up with a collections of expansion histories all yielding a better $\chi^{2}$ to the data than the best-fit $\Lambda$CDM model, and as such, are a non-exhaustive sample of plausible expansion histories indistinguishable from $\Lambda$CDM. The top-left panel of Fig. 1 shows the residuals in blue with respect to the best-fit $\Lambda$CDM model. The different smoothing iterations are shown in colour-codes, from dark blue (first iteration) to yellow (last iteration). Each curve is a viable fit to the data, with $\chi^{2}$ equal to or better than the $\Lambda$CDM fit. We emphasize that these should be viewed as sample viable cosmologies and not a confidence region.

The data at $z\gtrsim 1$ yield negative residuals on average, pulling the reconstructed $\hat{\mu}(z)$ toward lower values with respect to the $\Lambda$CDM best-fit. These low distance moduli in turn drive the corresponding sharp rise in the Hubble parameter $h$ seen in the upper right panel. The bottom-left panel shows the $Om$ diagnostic, a model independent statistic defined as (Sahni et al., 2008)

[TABLE]

For a flat $\Lambda$CDM universe, $Om$ is constant and equal to the matter density parameter $Om(z)\equiv\Omega_{\text{m}}$. In the panel $Om$ can clearly deviate from constancy, increasingly so at high redshift for successive iterations. That is, the data tend to pull the smooth reconstructions away from $\Lambda$CDM, reflecting higher $h$, lower distances, and, in a less model independent sense, more dark energy at higher redshifts, corresponding to less negative equation of state $w$.

We explore three origins for the possible deviations in the model independent reconstructions in the following sections: 1) the data accurately reflect the underlying (non-$\Lambda$CDM) cosmology and theoretical possibilities, 2) the data reflect differences in the high and low redshift supernova samples, i.e. some shift in source properties, and 3) the data are influenced by survey properties at high redshift such as declining detection efficiency (e.g. incompletely corrected Malmquist bias). In other words, if the data are taken at face value, then it has the following implications for cosmology (Sec. 3) or source properties (Sec. 4), while Sec. 5 explore the possibility of the survival of residuals after the correction.

3 Cosmological Interpretation

The reconstructed model-independent Hubble parameter closely traces the $\Lambda$CDM behavior out to $z\approx 1$ and then it can rapidly increase above it at higher redshift. This implies an extra energy density relative to $\Lambda$CDM at high redshift, which can be interpreted as a dark energy that fades away less strongly than a cosmological constant. In order to preserve the $\Lambda$CDM behavior at $z\lesssim 1$, the dark energy equation of state has to be near $w\approx-1$ at low redshift before becoming less negative at early times.

The apparent severity of the cosmology shift is perhaps more clearly appreciated if we relax model independence to the extent of defining a separate matter component with $\Omega_{\text{m}}$ given by the best fit $\Lambda$CDM and plotting the effective dark energy equation of state $w(z)$ in the lower right panel of Fig. 1. We can compute the equation of state behavior $w(z)$ from the Hubble parameter by

[TABLE]

Clearly the results show that this is not a simple evolution as for a freezing quintessence, or even a transition from an early dark energy scaling as matter.

Such $w(z)$ behavior does not correspond to any standard dark energy model. A sharp transition is required at $z\approx 1$ from nearly a cosmological constant behavior ($w=-1$) to large positive values of $w$ (this is independent of the exact value of $\Omega_{\text{m}}$ taken). Even phenomenological models with a rapid (much faster than Hubble time) transition to early dark energy matter behavior ($w=0$, Shafieloo et al., 2018b), or more extremely stiff matter behavior ($w=1$, Zeldovich, 1972), do not provide close fits. Moreover, any such model that approaches $w\approx 0$ by $z\approx 1.5$ would unviably alter the growth of structure and the integrated Sachs-Wolfe effect in the cosmic microwave background.

A purely cosmological interpretation of the high-redshift behaviour carries grave challenges. Such extreme expansion histories require dark energy models that would violate matter domination, and thus do not seem viable overall. Therefore, rather than cosmology, we look for systematics in the data and consider in the next section source and survey characteristics as origins of the model-independent deviation from $\Lambda$CDM.

There is one other cosmological source for the deviation worth considering first. Gravitational lensing changes the magnitudes of sources, increasingly so for higher redshift sources. While for a large number of sources at a given redshift the apparent distance is unaffected at linear order, for a small number of sources they will be preferentially brightened as a selection effect. The average dispersion effect is corrected for in the Pantheon compilation. Let us briefly consider how much the correction would need to be off to undo the deviation from $\Lambda$CDM in Fig. 1.

The Pantheon analysis assumes the form $\sigma_{\rm lens}=0.055z$ from Jönsson et al. (2010). If instead we took the form $0.093z$ from Holz & Linder (2005), or the form more complete at higher redshift, $0.1z/(1+0.07z)$ from Aldering et al. (2007), then we would expect an uncorrected deviation of 0.053 at $z=1.5$. Interestingly, we see in Sec. 5 that the deviation from $\Lambda$CDM corresponds to a deviation of roughly 0.06 at $z=1.5$. We do not claim this is the cause of the deviation – every source would need to be magnified at the $\sim 1\sigma$ level – but we urge caution in treating $z>1$ distance indicators without a robust lensing probability distribution function, and look forward to surveys that will put sufficient $z=1$–1.7 SNIa on the Hubble diagram to allow for safety in numbers from lensing.

4 Supernova Property Interpretation

In this section, we focus on supernovae or environmental effects that might bias the absolute magnitude through drift of the observed supernova population. One possibility is the metallicity of the host galaxy. Moreno-Raya et al. (2016) found this could affect the peak magnitude of a supernova as

[TABLE]

where $Z_{\odot}$ is solar metallicity. The magnitude bias as a function of the host metallicity is shown in the left-hand panel of Fig. 2.

A constant offset would not affect supernova cosmology, but an unaccounted for redshift evolution could. The cosmic mean metallicity $Z_{b}$ as a function of redshift is given by (Madau & Dickinson, 2014)

[TABLE]

The yield $y$ and the return rate $R$ depend on the initial mass function, and for a Salpeter IMF, are about $y=0.02$ and $R=0.27$.

Combining equations (8) and (9), we obtain the absolute magnitude offset as a function of redshift, shown in the right-hand panel of Fig. 2. Following Belczynski et al. (2016), we shifted the relation by {0,+0.5,+1} dex in blue (lower), orange (middle), and green (upper), accounting for the uncertainty in the redshift zero metallicity expression in comparison with data.

None of these potential residual corrections match the shape of the empirical correction Eq. (13) of the next section. However, we report in Table 1 the impact of such metallicity corrections on $\Omega_{\text{m}}$ and their associated $\Delta\chi^{2}$ with respect to the $\Lambda$CDM uncorrected best-fit. It is worth noting that none of these corrections yield a statistically significant improvement to $\chi^{2}$. In addition, note the case with the shift of +1 dex, which is closest in magnitude to the empirical correction of the next section, yields a best-fit $\Omega_{\text{m}}=0.251$, which is significantly lower. Therefore, these considerations do not seem to favour unaccounted-for metallicity evolution as the cause for the blowup of the model-independent reconstructions of $h$ at $z\geq 1$.

5 Survey Property Interpretation

Survey selection effects also can affect the Hubble diagram of the distance-redshift relation, and the cosmological quantities derived from it. The most well known example is Malmquist bias, where supernovae near the upper redshift limit of the survey have a preferential selection effect for the intrinsically brighter members of the population due to the survey flux and signal to noise limits. Other possible effects include filter zeropoints and color corrections, though we note the Pantheon compilation has put great effort into calibration (Jones et al., 2018; Scolnic et al., 2018). Recent techniques for dealing with sample selection generally are discussed in Hinton et al. (2017, 2018).

To mock up the effect of a high redshift selection bias (which we do not suggest is a Malmquist bias, simply that high redshift is the most likely place for bias), we consider a simple linear trend above a certain redshift, where the selection starts to bias the magnitude. In other words, we are considering the possibility of residuals after the current best efforts at correction. To choose the transition redshift we examine the data. Figure 3 shows the data density as a function of redshift; one might expect that selection effects are increased by a sparse sampling with just a few supernovae per redshift interval. The data density drops to low values for $z\gtrsim 1$, with just a few points per interval of $\Delta z=0.1$ for $z=1$–2. In addition, in the lower-left panel of Fig. 1, $Om$ is consistent with a constant level at low-$z$, and starts to depart from constant after $z\gtrsim 1$, suggesting that any new element enters only there. Therefore, we apply the following empirical correction at $z>z_{c}=1$ (corresponding to 23 points):

[TABLE]

where

[TABLE]

and where $A$ is a free parameter to fit.

Figure 4 shows $\Delta\chi^{2}=\chi^{2}-\chi^{2}_{0}$, where $\chi^{2}_{0}$ is the $\chi^{2}$ of the best-fit flat-$\Lambda$CDM model to the uncorrected Pantheon data. The best $\chi^{2}$ is located at $A=0.12$ with $\Delta\chi^{2}=-0.9$, which is not statistically significant for one added degree of freedom (dof). Since the correction only affects the last 23 data points, it has little impact on the best-fit cosmology ($\Omega_{\text{m}}$= 0.293 versus 0.298). In the following, we use this value of $A=0.12$ to “correct" the data.

After finding the best-fit flat $\Lambda$CDM cosmology for the $\mu$ corrected according to Eq. (13), we implement the smoothing method to the corrected data using the best-fit value for $A$ to obtain the cosmological quantities. Figure 5 shows the same plots as Fig. 1, with the corrected data (for $A=0.12)$. The corrected data points are shown in orange in the top-left panel. In the other three panels, the reconstructions then show tighter agreement in behavior out to higher redshift, and do not lead to a blow-up in $h$ or $Om$ for almost all the redshift range. We caution that the dark energy equation of state $w$, being a second derivative of the data, is always more sensitive to statistical fluctuations.

To quantify the significance of the deviation at high-redshift, we use the crossing statistics (Shafieloo et al., 2011; Shafieloo, 2012). The 0th mode of the statistics is

[TABLE]

which is related to the pull, while the $N$th mode is the $\chi^{2}$. Table 2 shows $T_{0}$ and $\chi^{2}$ for the 23 points with $z>z_{\text{c}}$.

The smoothing applied to the corrected data yields a similar $\chi^{2}$, but its associated $T_{0}$ is much lower: the data are more evenly distributed around 0. (A useful crosscheck is that the uncorrected data prefer the uncorrected smoothing, and the corrected data prefer the corrected smoothing.) To assess the significance of the obtained $T_{0}$, we simulated 1000 realizations of a Gaussian random variable $\sim\mathcal{N}(0,1)$, and measured the distribution of $T_{0}$. The expected $1\sigma$ of the distribution of $T_{0}$ is about 32. Therefore, all values reported here are within $1\sigma$. We note that the $\chi^{2}$ for the smoothing of the corrected data (15.45) is slightly larger than for the $\Lambda$CDM fit to the same data (15.35). The reason is that the smoothing algorithm considers the total $\chi^{2}$ to the whole data (to improve the fits iteratively), while here we are just comparing the $\chi^{2}$ for the last 23 data points. The correction (13) has in fact made $\Lambda$CDM such a good fit to the data at $z>1$ that smoothing is not really improving much in this range.

As a further, independent test to assess whether the data prefer any correction, we also applied Gaussian Processes to the residuals. The marginal likelihood does not have a clear peak for the correlation length and does not prefer a nonzero amplitude, meaning that the data are consistent with $\Lambda$CDM and do not suggest need for a Malmquist-like correction despite the visual impact of Fig. 1 (see Appendix A).

6 Summary and Conclusion

Model independent reconstructions of the expansion history from the Pantheon supernovae compilation allow substantial deviation from the standard $\Lambda$CDM model at redshifts greater than 1. We studied three possibilities for these deviations.

For cosmology, the most radical theoretical possibilities, that dark energy can behave as matter or stiff matter at early stages, do not behave in as extreme a manner as cases of the reconstructions. Similarly, while gravitational lensing can cause apparent cosmology deviations, the magnitude of the apparent effect seems too large for this, though better knowledge of the lensing distribution function is needed. This makes it highly unlikely that the cause of the extreme deviations from the $\Lambda$CDM model should be interpreted as cosmological.

As for source properties, we studied the possible effect of the host galaxy metallicity redshift evolution on the peak magnitude of the type Ia supernovae and their effect on the reconstructions. We derived the magnitude bias - redshift relation, and investigated its dependence on the still incompletely determined zero redshift offset. This mechanism also could not reproduce the reconstructions at high redshifts while preserving the low redshift cosmology, though it could have more modest effects.

Regarding survey characteristics, we explored an empirical high-redshift correction of the distance moduli in the Pantheon compilation. Such selection effects, such as Malmquist bias, are often difficult to fully correct the data for. Using a simple linear correction for $z>z_{c}=1$, we showed the blow-up of the expansion history at high redshifts could be strongly reduced. While visually the effect of the corrections on the reconstructions are substantial, the improvement on the $\chi^{2}$ fit (around $0.9$) is marginal. Furthermore, using Gaussian Process regression as well as crossing statistics, we examined the significance of the improvement.

Our results show that the correction is not required statistically: current Pantheon data is consistent with $\Lambda$CDM despite significant deviations also being allowed in a model independent expansion history. We eagerly look forward to more data at $z\gtrsim 1$ that would have greater leverage on confirming $\Lambda$CDM – or pointing toward significant effects in cosmology, gravitational lensing, source property evolution, or survey selection effects.

Acknowledgements

This work benefited from the Supercomputing Center/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2016-C2-0035 and KSC-2017-C2-0021) and the high performance computing clusters Polaris and Seondeok at the Korea Astronomy and Space Science Institute. A.S. would like to acknowledge the support of the National Research Foundation of Korea (NRF- 2016R1C1B2016478). A.S. would like to acknowledge the support of the Korea Institute for Advanced Study (KIAS) grant funded by the Korea government. EL is supported in part by the Energetic Cosmos Laboratory and by the U.S.  Department of Energy, Office of Science, Office of High Energy Physics, under Award DE-SC-0007867 and contract no. DE-AC02-05CH11231.

Appendix A GP regression

We applied Gaussian Processes (GP) regression (Rasmussen & Williams, 2006) to the residuals (with respect to the best-fit, uncorrected $\Lambda$CDM model) at $z>z_{c}$. GP have been widely used in the community to reconstruct the expansion history, the equation of state of dark energy, the cosmic growth rate, and more (see, e.g., Holsclaw et al., 2010a, b, 2011; Shafieloo et al., 2012, 2013).

Starting from a training set of points $(\mathbf{x},\mathbf{y}=f(\mathbf{x})+\mathbf{\varepsilon})$ where $\varepsilon$ is a Gaussian noise with mean zero and covariance C, and $\mathbf{y}$ are the measured values, we can model $f$ as a stochastic process with covariance K, and reconstruct $\mathbf{f}_{*}(\mathbf{x}_{*})$ at the test points $x_{*}$.

The joint-distribution of the training outputs $\mathbf{y}$ and the test output $\mathbf{f}_{*}$ is a Gaussian distribution given by

[TABLE]

where C is the covariance of the data.

For a given kernel, the covariance between pairs of random variables $\mathbf{u}$ and $\mathbf{v}$ is thus given by$\textbf{{K}}(f(\mathbf{u}),f(\mathbf{v}))=k(\mathbf{u},\mathbf{v})$, where $k(\mathbf{u},\mathbf{v})$ is the covariance kernel. We use the squared exponential kernel defined as

[TABLE]

where $(\sigma_{f},\ell)$ are two hyperparameters controlling the amplitude and the correlation scale.

The log marginal likelihood (LML) is given by

[TABLE]

A clear preference for $\sigma_{f}>0$ would mean that the data suggest that a correction to the mean function (i.e. $\Lambda$CDM) is needed. Fig. 6 shows the log marginal likelihood (LML) of the hyperparameters $(\sigma_{f}^{2},\ell)$. In all cases of the uncorrected and corrected, and unsmoothed and smoothed data, the shapes of the LML are similar, with the corrected case having a slightly higher likelihood. It is maximal at low values of $\sigma_{f}$, i.e.  consistent with no correction, suggesting no deviation from the mean function in both corrected and uncorrected cases.