Null test for interactions in the dark sector

TL;DR

This paper introduces a null test for interactions in the dark sector based on the evolution of the dark energy and dark matter densities, using model-independent machine learning methods on recent cosmological data.

Contribution

It proposes a new null test for dark sector interactions and applies machine learning techniques to observational data to assess deviations from standard cosmology.

Findings

Standard $ m extLambda$CDM is consistent within 3σ

Some deviations at low redshift suggest possible dark sector interactions

Potential link to H0 tension in cosmology

Abstract

Since there is no known symmetry in Nature that prevents a non-minimal coupling between the dark energy (DE) and cold dark matter (CDM) components, such a possibility constitutes an alternative to standard cosmology, with its theoretical and observational consequences being of great interest. In this paper we propose a new null test on the standard evolution of the dark sector based on the time dependence of the ratio between the CDM and DE energy densities which, in the standard $\Lambda$CDM scenario, scales necessarily as $a^{-3}$. We use the latest measurements of type Ia supernovae, cosmic chronometers and angular baryonic acoustic oscillations to reconstruct the expansion history using model-independent Machine Learning techniques, namely, the Linear Model formalism and Gaussian Processes. We find that while the standard evolution is consistent with the data at $3\sigma$ level, some deviations from the $\Lambda$CDM model are found at low redshifts, which may be associated with the current tension between local and global determinations of $H_0$.

Figures

Tables

Table 1. Table 1: Angular BAO data.

Catalog	$z$	$\theta \left(z\right)$	${\sigma }_{\theta \left(z\right)}$	Ref.
SDSS-DR7	0.235	9.06	0.23	Alcaniz:2016ryy
SDSS-DR7	0.365	6.33	0.22	Alcaniz:2016ryy
SDSS-DR10	0.450	4.77	0.17	Carvalho:2015ica
SDSS-DR10	0.470	5.02	0.25	Carvalho:2015ica
SDSS-DR10	0.490	4.99	0.21	Carvalho:2015ica
SDSS-DR10	0.510	4.81	0.17	Carvalho:2015ica
SDSS-DR10	0.530	4.29	0.30	Carvalho:2015ica
SDSS-DR10	0.550	4.25	0.25	Carvalho:2015ica
SDSS-DR11	0.570	4.59	0.36	Carvalho:2017tuu
SDSS-DR11	0.590	4.39	0.33	Carvalho:2017tuu
SDSS-DR11	0.610	3.85	0.31	Carvalho:2017tuu
SDSS-DR11	0.630	3.90	0.43	Carvalho:2017tuu
SDSS-DR11	0.650	3.55	0.16	Carvalho:2017tuu
SDSS-DR12Q	2.225	1.77	0.31	deCarvalho:2017xye

Table 2. Table 2: Result of the statistical analysis with the type Ia SN data from the Pantheon sample.

	$H_0$	${\omega }_b$	$M$	${\Omega }_{c0}$
R18/BBN	$73.52$	$0.0223$	$-{19.25}_{-0.01}^{+0.01}$	${0.256}_{-0.023}^{+0.022}$
Planck	$67.66$	$0.02242$	$-{19.42}_{-0.01}^{+0.01}$	${0.249}_{-0.023}^{+0.022}$

Table 3. Table 3: Analysis based on the χ ν 2 subscript superscript 𝜒 2 𝜈 \chi^{2}_{\nu} for the full CC dataset.

${\alpha }_{\max }$	${\chi }_{\nu }^2$	$3\sigma$ interval	$\nu$
0	4.1	$[0.37,1.90]$	30
1	0.57	$[0.36,1.92]$	29
2	0.53	$[0.35,1.94]$	28
3	0.48	$[0.34,1.96]$	27

Table 4. Table 4: Analysis based on the χ ν 2 subscript superscript 𝜒 2 𝜈 \chi^{2}_{\nu} for the full type Ia SNe dataset.

${\alpha }_{\max }$	${\chi }_{\nu }^2$	$3\sigma$ interval	$\nu$
1	5.3	$[0.42,1.79]$	38
2	1.2	$[0.42,1.80]$	37
3	1.1	$[0.41,1.81]$	36
4	1.2	$[0.40,1.83]$	35

Table 5. Table 5: Analysis based on the χ ν 2 subscript superscript 𝜒 2 𝜈 \chi^{2}_{\nu} for the full angular BAO dataset.

${\alpha }_{\max }$	${\chi }_{\nu }^2$	$3\sigma$ interval	$\nu$
1	1.0	$[0.14,2.53]$	12
2	1.1	$[0.12,2.62]$	11

Equations

$${\dot{\rho }}_c+3H{\rho }_c$$

$${\dot{\rho }}_x$$

$$Q=3H\gamma R({\rho }_c,{\rho }_x),$$

$$r\equiv \frac{{\rho }_c}{{\rho }_x}\Rightarrow \dot{r}=r\left(\frac{{\dot{\rho }}_c}{{\rho }_c}-\frac{{\dot{\rho }}_x}{{\rho }_x}\right).$$

$$\dot{r}+3Hr\left(\gamma R\frac{{\rho }_c+{\rho }_x}{{\rho }_c{\rho }_x}+1\right)=0.$$

$$\dot{r}+3Hr\big{[}\gamma f\left(r\right)+1\big{]}=0\,.$$

$${\rho }_d$$

$$p_d$$

$${\rho }_d={\rho }_c\left(1+\frac{1}{r}\right)\text{mbox}or{\rho }_d={\rho }_x\left(1+r\right),$$

$$p_d=w_d{\rho }_d\text{with}{w}_d\equiv -\frac{1}{1+r},$$

$${\dot{\rho }}_d+3H\left(1+w_d\right)=0.$$

$$r(z)=r_0(1+z{)}^3\text{with}{r}_0=\frac{{\Omega }_{c0}}{{\Omega }_{x0}},$$

$$\frac{H^2}{H_0^2}={\Omega }_{d0}\frac{1+r_0{\left(1+z\right)}^3}{1+r_0}+{\Omega }_{b0}{\left(1+z\right)}^3,$$

$$r_0(z)=\frac{1-{\Omega }_{b0}+{\Omega }_{b0}{\left(1+z\right)}^3-H^2/H_0^2}{H^2/H_0^2-{\left(1+z\right)}^3}.$$

$$\mu \left(z\right)=5\log \left[\frac{d_L\left(z\right)}{1\text{ Mpc}}\right]+25,$$

$$d_L\left(z\right)=\frac{c}{H_0}\left(1+z\right)\mathcal{D}\left(z\right),$$

$$\mathcal{D}\left(z\right)=\frac{H_0}{c}{\left(1+z\right)}^{-1}10^{\frac{\mu \left(z\right)}{5}-5}.$$

$$\mathcal{D}\left(z\right)={\int }_0^z\frac{d\tilde{z}}{E\left(\tilde{z}\right)}\Rightarrow E\left(z\right)=\frac{1}{{\mathcal{D}}^{\prime }\left(z\right)}.$$

$$\theta \left(z\right)=\frac{r_s}{d_A\left(z\right)\left(1+z\right)}\left(\frac{180}{\pi }\right),$$

$$d_A\left(z\right)=\frac{c}{H_0\left(1+z\right)}\mathcal{D}\left(z\right).$$

$$\mathcal{D}\left(z\right)=\frac{H_0}{c}\frac{r_s}{\theta \left(z\right)}\left(\frac{180}{\pi }\right).$$

$$H_0$$

$${\omega }_b$$

$$H_0$$

$${\omega }_b$$

$$r_s=146.6\pm 4.1\text{ Mpc}.$$

$$r_s=147.21\pm 0.23\text{ Mpc}.$$

$$t(z,\{c_{\alpha }\})=\sum _{\alpha =0}^{{\alpha }_{\max }}{c}_{\alpha }{g}_{\alpha }(z),$$

$$d_i=t_i+e_i,$$

$$F_{\alpha \beta }=g_{\beta i}{\Sigma }_{ij}^{-1}{g}_{\alpha j}$$

$$c_{\alpha ,\mathrm{bf}}=F_{\alpha \beta }^{-1}{B}_{\beta }\equiv {\Sigma }_{\alpha \beta }{B}_{\beta },$$

$${\sigma }_{r_0}^2=J_{\alpha }{S}_{\alpha \beta }{J}_{\beta },$$

Full text

Dark Energy, Dark Matter, Cosmological Parameters, Distance Scale

I Introduction

According to the standard model of cosmology about 5% of the energy content of the universe is made of particles belonging to the standard model of particle physics. The remaining 95% is attributed to the so-called dark sector. Roughly 25% are thought to consist of a yet-undetected cold dark matter component (CDM), while dark energy (DE), the fuel that drives the current cosmic acceleration, would be responsible for the missing 70%. The fact that DE and CDM have comparable energy densities today – the so-called coincidence problem – has motivated the study, at great depth, of dynamical models of DE that feature interactions between dark energy and dark matter (see e.g. Copeland:2006wr, ; 2010deto.book…..A, ; Li:2011sd, ; Chen:2010dk, , and references therein), hoping to shed light on the nature of the dark sector. In this context, an important topic of research that has been extensively explored relies on considering some specific models for this interaction between the dark components in order to assess its cosmological consequences Amendola:1999er ; CalderaCabral:2008bx ; Alcaniz:2012mh ; Marttens:2016cba ; Marra:2015iwa ; Zimdahl:2014jsa ; Gonzalez:2018rop .

Here, we adopt a different approach, i.e., instead of constraining the interaction parameter of a specific model we present a model-independent way to investigate whether or not such a interaction in the dark sector really exists. For this, we introduce a new null test that is sensitive to the existence of a possible interaction between the dark components. Equivalently, if this null test is failed, then one may suspect that there may be new physics beyond the standard model and, in particular, that dark matter and dark energy are not independent entities. In other words, this null test has the ability to extract information that one may miss when the analysis performed is restricted to parameter estimation within a specific class of interacting dark energy models.

The proposed null test is based on the time dependence of the ratio between CDM and DE energy densities, i.e., $r(z)=\rho_{\rm{CDM}}/\rho_{\rm{DE}}$, which in the $\Lambda$CDM model is given by $r\left(z\right)=r_{0}\left(1+z\right)^{3}$, where $r_{0}$ is the current value of this ratio. Since an interaction in the dark sector affects the dynamics of the components involved, this quantity is directly sensitive to the existence of such interaction. In order to carry out this new null test we will reconstruct the expansion history of the universe, in a model-independent way, using Machine Learning (ML) techniques applied to cosmic chronometers (CC) measurements, type Ia supernovae (SNe Ia) data and also angular Baryon Acoustic Oscillation (BAO) determinations. In particular we will use the Linear Model formalism (LM) and Gaussian Processes (GP). Regarding LM, we improve the so-called “learning curve” methodology by generalizing the “Mean Square Error” (MSE) and “Mean Square Prediction Error” (MSPE) to the case of data that have an arbitrary covariance matrix and by taking into account the covariance matrix on the model parameters obtained from the training set. We call this generalization the “calibrated learning curves”.

The outline of this paper is as follows. In Section II we present the theoretical description of a rather general class of unified/interacting models, which is used to derive the $r\left(z\right)$ null test in Section III. The Section IV is devoted to present the datasets used to perform the null test and to discuss about the priors on the further “external” parameters. In the Sections V and VI the two methods used to perform the null test, i.e., LM and GP are discussed. The results are presented in Section VII, while Section VIII is dedicated to conclusions. The paper has also two appendices: Appendix A where we present the theoretical basis of the calibrated learning curves, and Appendix B where the python script learning_curve is released as an automatic tool to compute and plot the (calibrated) learning curves for any given dataset with a covariance matrix.111The python scripts can be downloaded from github.com/rodrigovonmarttens/learning_curve.

II Interacting/Unified description

A simple and viable alternative to the standard cosmological model is to consider an interaction between the dark components of the universe. The unknown nature of the dark sector does not allow us to provide a microphysical description of this interaction, which can only be modeled phenomenologically via a source term in the energy conservation equation,

[TABLE]

From now on the subscripts $c$ and $x$ denote the dark matter and dark energy components, respectively, the dot denotes the derivative with respect to the cosmic time, and it was assumed that the dark energy equation of state (EoS) parameter is $w_{x}=-1$. In the individual conservation equations above the interaction has opposite sign so that the total energy–momentum tensor is conserved.

Let us now assume that the interaction source can be parametrized via

[TABLE]

where the dimensionless parameter $\gamma$ gives the interaction strength and the function $R$ specifies the type of interaction (see vonMarttens:2018iav for details). It is then convenient to introduce the ratio $r$:

[TABLE]

Note that this quantity can be directly associated to the cosmic coincidence problem Velten:2014nra . Substituting equations (1) and (2) in the equation above, it is possible to obtain a differential equation for $r$,

[TABLE]

We are interested in the case in which the ratio between CDM and DE energy densities depends only on the scale factor (or, equivalently, on the redshift) and we assume, as an ansatz, that the first term in the parenthesis is a function of $r$, that is,

[TABLE]

This formalism to describe interactions in the dark sector was introduced in vonMarttens:2018iav and is at the core of the null test proposed in this work.

If equation (6) has a solution that depends only on the scale factor, then we have an interacting model that can be associated to a unified model, i.e., we can combine the CDM and DE components in order to describe a single dark fluid. For this dark fluid, we define its energy density and pressure as the sum of the energy densities and pressures of CDM and DE,

[TABLE]

In equation (7) one can express the unified dark energy density in terms of $r$ and only one of the energy densities of the dark sector’s components:

[TABLE]

and equation (8), using the second equation in (9), can be rewritten as:

[TABLE]

where we defined the dark EoS parameter $w_{d}$. Note that equation (10) is completely general, but it describes a unified model only if $r=r\left(a\right)$. Since this dark fluid must be conservative, energy conservation must be satisfied,

[TABLE]

This unified description for the dark sector has been extensivly explored in the literature Zimdahl:2010tt ; Zimdahl:2011mb ; Marttens:2017njo .

III The null test $r_{0}(z)$

For convenience, from now on, we will use the redshift $z$ instead the scale factor $a$. The $\Lambda$CDM model is recovered if

[TABLE]

and, in order to obtain a null test for interactions in the dark sector, we use (10) together with (12) so that the Friedmann equation becomes:

[TABLE]

where we have assumed spatial flatness, so that $\Omega_{d0}\equiv\Omega_{c0}+\Omega_{x0}=1-\Omega_{b0}$, and we have neglected radiation because we will consider only low-redshift data. The above equation is the Friedmann equation for $\Lambda$CDM in terms of $r_{0}$ and $\Omega_{d0}$. Note that, since we are interested in an interacting scenario in which the interaction affects only the dark components, it is necessary to describe baryons and dark matter separately Marttens:2014yja .

The null test for interacting models is obtained solving equation (13) for $r_{0}$:

[TABLE]

Within the standard flat $\Lambda$CDM model one expects a constant $r_{0}(z)={\Omega_{c0}}/{\Omega_{x0}}$. If a deviation is detected one may suspect not only that $\Lambda$CDM is falsified but also that dark matter and dark energy are not independent entities.

IV Cosmological data

In order to carry out the null test of equation (14), it is necessary to determine the three quantities $H(z)$, $H_{0}$ and $\Omega_{b0}$. The expansion history $H(z)$ will be reconstructed using CC, SNe Ia and BAO data. The parameters $H_{0}$ and $\Omega_{b0}$ will be discussed in Section IV.4.

IV.1 Cosmic Chronometers

Cosmic chronometers are passively evolving old galaxies whose redshifts are known, and the expansion history of the universe can be inferred directly from their differential ages Zhang:2012mp ; Simon:2004tf ; Moresco:2012jh ; Moresco:2016mzx ; Ratsimbazafy:2017vga ; Stern:2009ep ; Moresco:2015cya . Here we will adopt the latest data as presented in (Marra:2017pst, , Table I). Figure 1 illustrates the data points as well as the fitted Hubble function obtained using the linear model formalism presented in Section V at the first, second and third order.

IV.2 Type Ia Supernovae

The second dataset that we use to reconstruct $H(z)$ is the compressed supernova Ia Pantheon compilation (40 bins) Scolnic:2017caz .222All the data (binned and full), as well as their covariance matrices, can be downloaded from github.com/dscolnic/Pantheon. Note that, since we are performing a null test, the fact of we are using the binned catalog is not a problem in the sense of favoring the $\Lambda$CDM model.

Type Ia Supernovas provide determinations of the distance modulus $\mu$, whose theoretical prediction is related to the luminosity distance $d_{L}$ according to:

[TABLE]

where the luminosity distance is given in Mpc. In the standard statistical analysis, one adds to the distance modulus the nuisance parameter $M$, an unknown offset sum of the supernova absolute magnitude (and other possible systematics), which is degenerate with $H_{0}$. In this analysis, as will be discussed in more details in section IV.4, the value of $M$ is related to the prior on $H_{0}$. As we are assuming spatial flatness, the luminosity distance is related to the comoving distance $\mathcal{D}$ via

[TABLE]

where $c$ is the speed of light, so that, using (15), one obtains:

[TABLE]

Finally, the normalized Hubble function $E\left(z\right)\equiv H\left(z\right)/H_{0}$ can be obtained by taking the inverse of the derivative of $\mathcal{D}(z)$ with respect to the redshift (denoted with a prime):

[TABLE]

The binned Pantheon data points (subtracting $M=-19.25$) are shown in Figure 2 together with the fitted distance modulus obtained using the linear model formalism at the first, second and third order.

IV.3 BAO

The last dataset that we use to reconstruct $H(z)$ are model-independent angular BAO determinations obtained using the angular correlation function Sanchez:2010zg . In this case, we use 14 uncorrelated data points from Carvalho:2015ica ; Alcaniz:2016ryy ; Carvalho:2017tuu ; deCarvalho:2017xye , which are presented in Table 1 and illustrated in Figure 3 (with its first order fit). Model-independent determinations of the radial BAO scale were recently obtained in Marra:2018zzo .

The theoretical BAO angular scale, in degrees, is given by,

[TABLE]

where $r_{s}$ is the sound horizon of the primordial photon-baryon fluid at the drag time and $d_{A}\left(z\right)$ is the angular diameter distance, which, in a flat universe, is related to the comoving distance by:

[TABLE]

Substituting the equation above in equation (19), one obtains the following explicit relation between $\theta\left(z\right)$ and $\mathcal{D}\left(z\right)$,

[TABLE]

The normalized Hubble function $E(z)$ is then obtained using equation (18).

IV.4 Further parameters

As mentioned earlier, in order to perform the null test of equation (14), it is necessary to determine also $H_{0}$ and $\Omega_{b0}$.333When using SNe Ia data or angular BAO measurements, only $\Omega_{b0}$ is necessary as one obtains directly $H/H_{0}$. We will adopt two different sets of values for these parameters. The first set is related to model-independent measurements: the local determination of $H_{0}$ obtained from low-redshift SN Ia data calibrated with loca Cepheids Riess:2018byc , and the measurement of the baryon density parameter from big bang nucleosynthesis Pettini:2012ph ,

[TABLE]

The second set of values comes from the most recent results from Planck Aghanim:2018eyx . In this work, we use the results obtained with TT,TE,EE+lowE+lensing+BAO,

[TABLE]

When using SNe Ia, we have the additional nuisance parameter $M$. Since $H_{0}$ is fixed according to (22) or (24), we choose the respective values of $M$ from a statistical analysis of the $\Lambda$CDM model with the Pantheon data obtained by fixing $H_{0}$ to the values previously mentioned. In order to perform the analysis, we have used the statistical code MontePython Audren:2012wb ; Brinckmann:2018cvx . When $H_{0}$ is fixed to the local determination (22), we found the best-fit value $M=-19.25$, while when $H_{0}$ is fixed to the Planck result (24) we obtained $M=-19.42$. The complete result of the statistical analysis with $1\sigma$ confidence levels is shown in Table 2 and illustrated in Figure 4.

When using the angular BAO determinations, one has the sound horizon $r_{s}$ as an additional parameter that must be specified. Maintaining the idea of the two sets of values we choose for the first set the model-independent result obtained from low-redshift standard rulers Verde:2016ccp $r_{s}=\left(101.0\pm 2.3\right)h^{-1}$ Mpc. Using the local measurement of $H_{0}$ of (22), one has:

[TABLE]

For the second set of values, we use the same Planck determination from Aghanim:2018eyx :

[TABLE]

V Linear model formalism

Here, we will describe how to reconstruct the cosmological functions, and also their derivatives, using the Linear Model formalism (LM); see, also, Marra:2017pst ; Gomez-Valent:2018hwc ; Pinho:2018unz .

V.1 Linear Models

Let us choose a set of base functions $g_{\alpha}(z)$ whose linear combination will constitute the template function $t(z,\{c_{\alpha}\})$:

[TABLE]

where $\alpha$ is an integer. The assumption is that $t(z,\{c_{\alpha}\})$ can describe the actual functions that we want to reconstruct: $H(z)$, $\mu(z)$ or $\theta(z)$. Clearly, this is conditional to an appropriate choice of $g_{\alpha}(z)$ and the order $\alpha_{\rm max}$. Usually, $g_{0}$ is a constant, often unity. The template will have $\alpha_{\rm max}+1$ coefficients.

Let us then assume that the data are given by:

[TABLE]

where $t_{i}=t(z_{i},\{c_{\alpha}\})$ and $e_{i}$ are Gaussian errors with covariance matrix $\Sigma_{ij}$.

Next, we fit the template $t$ to the data and use the LM formalism to calculate the Fisher matrix relative to the parameters $c_{\alpha}$, which gives an exact description of the likelihood as the template is linear in its parameters. The Fisher matrix is:

[TABLE]

where $g_{\alpha i}=g_{\alpha}(z_{i})$, and the best-fit values of $c_{\alpha}$ are:

[TABLE]

where $B_{\alpha}=d_{i}\Sigma^{-1}_{ij}g_{\alpha j}$ and we defined the covariance matrix $\Sigma_{\alpha\beta}$ on the parameters. Summarizing, this formalism had allowed us to exactly propagate the data covariance matrix $\Sigma_{ij}$ into the parameter covariance matrix $\Sigma_{\alpha\beta}$.

V.2 Error on null test reconstruction

The null test $r_{0}(z,\{\theta_{\alpha}\})$ is a nonlinear function of the template parameters and of the additional parameters of Section IV.4. The corresponding covariance matrix $S_{\alpha\beta}$ is obtained by forming an appropriate block diagonal matrix using the covariance matrices of the corresponding parameters (e.g., $\Sigma_{\alpha\beta}$). As we have chosen independent data, correlations among different datasets are not expected to be important.

In order to compute the error on $r_{0}(z,\{\theta_{\alpha}\})$ due to the uncertainty encoded in the covariance matrix $S_{\alpha\beta}$, a straightforward approach is to apply a change of variable from $\{\theta_{\alpha}\}$ to $r_{0}$. At the first order, the error is then given by:

[TABLE]

where

[TABLE]

V.3 Base functions

In order to reconstruct $H(z)$, $\mu(z)$ and $\theta(z)$, we will adopt the following base functions, respectively:

[TABLE]

In order to choose $\alpha_{\rm max}$ we will use the so-called “learning curves”, a machine learning tool.

V.4 Calibrated learning curves

The availability of large datasets is increasingly a defining feature of modern cosmology, in which data analysis has become an important component. Computations that were not possible a few decades ago can now be performed on GPU-based laptops. Machine learning includes a set of statistical techniques that allows computer systems to learn from examples, data, and experience, rather than following pre-programmed rules.

A simple method that is commonly used to choose the template order $\alpha_{\rm max}$ is the computation of the reduced chi-square $\chi^{2}_{\nu}$:

[TABLE]

where $d_{i}$ are the $N_{\rm tot}$ data of the full dataset $D$ with covariance matrix $\Sigma_{ij}$, and $\bar{t}_{i}=t(z_{i},\{c_{\alpha,{\rm bf}}\})$ where $c_{\alpha,{\rm bf}}$ are the best-fit parameters. If one finds that $\chi^{2}_{\nu}$ is compatible with its corresponding distribution with $N_{\rm tot}-\alpha_{\rm max}-1$ degrees of freedom, then the null hypothesis that $t(z,\{c_{\alpha}\})$ is the correct model is not rejected (it is “ruled in”). While powerful in its simplicity, this method is somewhat subjective as it strongly depends on the $p$-value threshold (e.g., $p=0.01$) that one is supposed to use. For example, two or more values of $\alpha_{\rm max}$ could be acceptable.

In order to overcome this difficultly and extract more information from the data, we will study the learning curves. These usually are used in contexts in which the data covariance matrix is not available and so a performance statistics with a known distribution (like $\chi^{2}_{\nu}$) cannot be built. Therefore, we will have to first generalize the standard learning curves to what we call the “calibrated learning curves.”

Let us then consider two disjoint subsets of the data set $D$ of $N_{\rm tot}$ elements: the training set $d$ and the validation set $\tilde{d}$.444From now on $d$ will refer to the training set and not to the full data set. The basic idea behind the learning curves consists in using the training set to fit the model and then test the latter with the validation set. Within machine learning the fit is usually obtained by minimizing the “Mean Squared Error” (MSE):

[TABLE]

where $N$ is the number of data points $d_{i}$ and $\bar{t}_{i}=t(z_{i},\{c_{\alpha,{\rm bf}}\})$, where $c_{\alpha,{\rm bf}}$ are the parameters that minimize the MSE. From now on we will denote the minimized $\overline{\text{MSE}}$ with just MSE. The test on the validation set is then performed by computing the “Mean Squared Prediction Error” (MSPE):

[TABLE]

The use of the new $\tilde{N}$ data points $\tilde{d}_{i}$ justifies the alternative name “out-of-sample mean squared error”. Note that $N_{\rm tot}\geq N+\tilde{N}$.

The learning curves are then the values of the MSE and the MSPE as a function of the training-set size $N$ while keep the validation-set size $\tilde{N}$ fixed. Usually, $\tilde{N}$ is 20–30% of $N_{\rm tot}$. The expectation is that the MSE will increase as the same number of parameters will be fitted to more data, and the MSPE will decrease as the training will produce a more reliable fitted model. In particular:

• •

an under-fitting model will feature converging but high (poor) MSE and MSPE;

• •

an over-fitting model will feature low MSE but high MSPE because the model is fitting the noise in the training set which is different with respect to the validation set;

• •

an optimal model will feature converging and low MSE and MSPE. Moreover, the sooner the convergence is reached, the better. Indeed, if the MSE and MSPE converge at $N<N_{\rm tot}-\tilde{N}$ reaching a plateau, it means that there were enough data to optimally train the model.

For more details, see, for instance hastie_09_elements-of.statistical-learning ; Abu-Mostafa ; Goodfellow-et-al-2016 ; Gron:2017:HML:3153997 .

The MSE and the MSPE do not use the data covariance matrix and, therefore, it is difficult to assess statistically their values. For example, it is not clear how to define “high”, “low” and “close to each other”. Therefore, we will calibrate the learning curve method by introducing new performance estimators that can be interpreted quantitatively in a statistical way.

A natural alternative to the MSE is the reduced chi-square function $\chi^{2}_{\nu}$,

[TABLE]

where $\Sigma$ is the covariance matrix of the training set $d$ and $\alpha_{\rm max}+1$ is the number of fitted parameters. Assuming that $t(z,\{c_{\alpha}\})$ is the correct model and that the data $d_{i}$ are distributed according to a multivariate Gaussian distribution of covariance matrix $\Sigma$, it is \big{\langle}\chi^{2}_{\nu}\big{\rangle}=1.

A natural alternative to MSPE is,

[TABLE]

where $\tilde{\Sigma}$ is the covariance matrix related to the validation set $\tilde{d}$. As these data were not used to obtain the best-fit parameters $\{c_{\alpha,{\rm bf}}\}$ the denominator only contains $\tilde{N}$. However, the expected value of $\tilde{\chi}^{2}_{\nu}$ is not unity as $\bar{t}_{i}$ is not the true value. Consequently, $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\nu}$ will not converge to the same numerical value. Here, we propose a new generalization of the MSPE:

[TABLE]

whose expectation value is unity as discussed in Appendix A.

In order to obtain smooth learning curves, we compute, for a fix $N$, $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ for 2000 partitions, from which we then compute mean and standard deviation. Note that the performance estimators $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ have an expectation value of unity independently of the training set size $N$, but this is true only if the expectation value is taken using independent training sets while here the training sets all come from the same dataset. In other words, the 2000 partitions are used to extract the average behavior of a training set of size $N$ from the full dataset $D$.

For smaller $N$ it is quite likely to obtain low values of $\chi^{2}_{\nu}$ as its distribution is skewed towards lower values, while for larger $N$ one expects $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ to converge to the common value of unity. If they converge reaching a plateau it means that $t(z,\{c_{\alpha}\})$ is the correct model and that the latter has been trained optimally by the data.

In our analysis, we adopt the following criterion in order to choose the best order $\alpha_{\rm max}$: the optimal $\alpha_{\rm max}$ is the one for which $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ converge fastest to unity with a plateau. It is very important to emphasize that this learning curve procedure is completely independent of any physical assumption, depending only on data.

V.5 Learning curve results

In the following, we present the results of the learning curve analysis for the datasets of Section IV.

V.5.1 Cosmic chronometers

For the cosmic chronometers, we divide the 31 data points in a training set with $N$ up to 20 and a validation set with $\tilde{N}=11$. Figure 5 shows the learning curves obtained with the template of equation (34) with $\alpha_{\rm max}=\{0,1,2,3\}$ (top to bottom).

The case $\alpha_{\rm max}=0$ is a clear case of under-fitting and is disfavored by the data: for $N=20$ the $\chi^{2}_{\nu}$ is well outside the corresponding 3$\sigma$ interval of $[0.26,2.17]$ (relative to the $\chi^{2}$ distribution with 19 degrees of freedom). This case corresponds to a constant Hubble rate, see equation (34). For the case $\alpha_{\rm max}=1$, $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ converge with a plateau to a value close to 1 and within the corresponding 3$\sigma$ interval $[0.24,2.21]$. According to our criteria, this is an optimal value of $\alpha_{\rm max}$.

The case $\alpha_{\rm max}=2$ is similar to the case $\alpha_{\rm max}=1$: $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ converge to the expected value with a plateau. Therefore, this case is also optimal, although $\tilde{\chi}^{2}_{\delta}$ converges to a value a little higher as compared with $\chi^{2}_{\nu}$, signaling a minor over-fitting. It is worth pointing out that the training curves of Figure 5 feature error bars (relative to the mean) computed, as mentioned before, from 2000 partitions. Therefore, the fact that there is a gap between $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ is statistically significant. It is also interesting to note how the learning curves characterize the models: the case $\alpha_{\rm max}=1$ is clearly simpler than the case $\alpha_{\rm max}=2$ as less data is necessary to train it (it converges faster).

The last case $\alpha_{\rm max}=3$ shows a lack of convergence with plateau, signaling that the model is too complex to be trained by the data. Therefore, we conclude that this case is disfavored by the data. Finally, we found that $\alpha_{\rm max}=1$ and $\alpha_{\rm max}=2$ are both acceptable. If we were to use the standard analysis based on the $\chi^{2}_{\nu}$ of equation (37) for the full dataset, we would have obtained the results presented in Table 3. According to these results $\alpha_{\rm max}=3$ is also acceptable, while the learning curve analysis disfavors it.

V.5.2 Type Ia Supernovae

We divide the Pantheon sample of supernovas in a training set with $N$ up to 28 data points and a validation set with $\tilde{N}=12$. Figure 6 shows the learning curves obtained with the template of equation (35) with $\alpha_{\rm max}=\{1,2,3,4\}$ (top to bottom). The case $\alpha_{\rm max}=1$ is a clear case of under-fitting and is disfavored by the data. Note that this case coincides with the first order in the cosmographic series expansion Cattoen:2007id . For the case $\alpha_{\rm max}=2$ and $\alpha_{\rm max}=3$, $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ converge with a plateau to a value close to 1 and within the corresponding 3$\sigma$ interval $[0.3,2.0]$. According to our criteria, these are optimal values of $\alpha_{\rm max}$. Finally, $\alpha_{\rm max}=4$ shows both a lack of convergence and of plateau, signaling over-fitting.

Therefore, we found that $\alpha_{\rm max}=2$ and $\alpha_{\rm max}=3$ are both acceptable. If we were to use the standard analysis based on the $\chi^{2}_{\nu}$ of equation (37) for the full dataset, we would have obtained the results presented in Table 4. According to these results $\alpha_{\rm max}=4$ is also acceptable, while the learning curve analysis disfavors it.

V.5.3 BAO determinations

With the BAO analysis we divide the 14 data points in a training set with $N$ up to 10 data points and a validation set with $\tilde{N}=4$. Figure 7 shows the learning curves obtained with the template of equation (36) with $\alpha_{\rm max}=\{1,2\}$ (top to bottom).

For the case $\alpha_{\rm max}=1$, $\chi^{2}_{\nu}$ and $\tilde{\chi}^{2}_{\delta}$ converge with a plateau to a value close to 1. According to our criteria, this is the optimal value of $\alpha_{\rm max}$. The case $\alpha_{\rm max}=2$ shows both a lack of convergence and of plateau, signaling over-fitting. We found, therefore, that only $\alpha_{\rm max}=1$ is acceptable. If we were to use the standard analysis based on the $\chi^{2}_{\nu}$ of equation (37) for the full dataset, we would have obtained the results presented in Table 5. According to these results $\alpha_{\rm max}=2$ is also acceptable, while the learning curve analysis disfavors it.

VI Gaussian Processes

A Gaussian Process (GP) is the generalization of the gaussian distribution of a random variable to the infinite function space. This mathematical approach has been successfully used as a nonparametric reconstruction method in cosmology since the pioneering works of Holsclaw2010 ; Holsclaw2011 . For instance, it has been applied to different datasets in order to calculate the dark energy equation of state Holsclaw2010 ; Holsclaw2011 ; Seikel2012 ; Zhang:2018gjb , the Hubble constant Gomez-Valent:2018hwc ; Busti:2014dua ; Licia2015 ; Li2015 , the cosmological matter perturbations Gonzalez2016 ; Gonzalez2017b ; Yin:2018mvu , and the gas deplection factor in galaxy clusters Holanda:2017cmc , among others.

A GP as a regression method is nonparametric. This means that their predictions are not restricted to a specific functional class (e.g., polinomial), but span an infinite family of classes with properties of continuity and differentiability. As this method is based on Bayesian statistics, we need to use prior and likelihood distributions to calculate the posterior distribution. Both prior and posterior distributions are defined via a mean function and a covariance matrix. The covariance quantifies the correlation between different functional values, $f(z)$ and $f(\tilde{z})$, at arbitrary independent variable points $z$ and $\tilde{z}$.

For the prior mean function we adopt the zero function as a conservative choice (this choice is recommended to avoid biased results) and, as commonly used in the literature, we choose square exponential covariance function:

[TABLE]

The so-called hyperparameters $\sigma_{f}$ and $l$ are related with the error/variation of the reconstruction and with its smoothness, respectively. These hyperparameters can be fixed by maximizing the likelihood distribution given the observational data (for a complete description of the GP method see Seikel2012 ; Rasmussen ). To perform the GP regression, we use the python package GaPP.555Available at acgc.uct.ac.za/$\sim$seikel/GAPP/index.html.

VII Results

Now, we present the reconstructions of $r_{0}(z)$ that we obtained using the methods discussed in the previous sections. The results are divided according to the method and the data used in order to reconstruct $r_{0}(z)$. In all plots, for comparison purposes, we include a dotted line corresponding to the value of $r_{0}$ predicted for the $\Lambda$CDM model according to the last results by the Planck satellite (Aghanim:2018eyx, ).

VII.1 Linear model results

Using the formalism developed in Section V, we obtain the following results:

Cosmic chronometers:

Figure 8 shows $r_{0}(z)$ using in equation (14) the Hubble rate function reconstructed via cosmic chronometer data. For the case $\alpha_{\rm max}=1$ one detects a deviation from the standard model when the model-independent priors discussed in Section IV.4 are used. This is clearly caused by the determination of the local Hubble constant (22) which is known to be in tension with the Planck indirect determination.666See Marra:2013rba ; Camarena:2018nbr for analyses that considered the effect of cosmic variance on local $H_{0}$. Furthermore, the $r_{0}(z)$ test seems to suggest that there is tension only for $z\lesssim 0.25$ while at higher redshift the local determination of $H_{0}$ gives a reconstruction in agreement with the value expected from Planck.

However, this interesting results loses its significance when $\alpha_{\rm max}=2$ is adopted, which was also found viable. Therefore, we conclude that better CC data are needed in order to conclude on this low-redshift tension.

Type Ia Supernovae:

Figure 9 shows the $r_{0}(z)$ test using the reconstruction of the distance modulus as in equations (17) and (18). We do not detect significant deviations from the Planck reference value.

BAO:

Figure 10 shows the $r_{0}(z)$ test using the reconstruction of the BAO angular scale (19) from which $E(z)$ is obtained via equations (21) and (18). When using Planck priors we detect a tension at $z\approx 0.3$ with respect to the Planck reference value (but not with a higher constant reference value). This again suggests a low-redshift deviation from the standard model.

VII.2 Gaussian Process results

Using the nonparametric reconstruction of $H(z)$ obtained via GP regression and the $H_{0}$ and $\omega_{b}$ priors described in Section IV.4, we calculate the null test $r_{0}(z)$ and its confidence levels by Monte Carlo sampling.

Cosmic chronometers:

In this case, the Hubble rate is reconstructed directly from the CC data which reduces the error propagation and the probability of wiggles in the $H(z)$ reconstruction. Figure 11 shows the $r_{0}(z)$ calculation using GP method and Monte Carlo sampling. It is evident the similarity of these results for the two sets of $\{H_{0},\omega_{b}\}$ priors with the results obtained using the LM formalism with $\alpha_{\rm max}=2$ (see Figure 8). This emphasizes the possibility that the tension in low-$z$ for the reconstruction via LM with $\alpha_{\rm max}=1$ may be not fundamental.

Type Ia Supernovae:

Figure 12 shows $r_{0}(z)$ from type Ia SN data using GPs. In this case, to calculate the null test we need to transform the distance modulus data to comoving distance and then reconstruct the derivative of this quantity to obtain the Hubble rate (see Section IV.2). The determination of the derivative of ${\mathcal{D}}$ propagates the error and its effect can be seen at high-$z$ where the density of the data is reduced. However, the null test is compatible with the $\Lambda$CDM model in the entire redshift range.

BAO:

Because of the large gap in redshift between the first thirteen data points and the last one, it is not possible to find a suitable GP reconstruction compatible with a cosmological scenario without assuming a nontrivial prior mean function.

VIII Conclusions

Interacting models of CDM and DE constitute an alternative description of the dark sector which have been largely investigated. In this paper we proposed a new null test in which any deviation of the standard cosmological scenario indicates a non-minimal interaction between these two dark components. Using linear models and Gaussian processes the expansion rate is reconstructed from the latest CC, SNe Ia and BAO data. For each formalism, the same analysis was performed using two sets of values for the “external” parameters $\{H_{0},\Omega_{b},r_{s},M\}$: the first set is based on model-independent results whereas the second one is obtained from the latest results of the Planck collaboration.

The test performed shows compatibility with the standard $\Lambda$CDM model within $3\sigma$ confidence level, but some cases deserve a careful analysis. For the LM analysis using the CC data, except the case in which $\alpha_{\rm max}=1$ and the model-independent priors are used, all the other results are compatible with a constant value for $r_{0}\left(z\right)$. The latest Planck result is satisfied in all cases for $z\gtrsim 0.25$. For $z\lesssim 0.25$ only the result with Planck priors reaches satisfactory results. When the model-independent priors are used, there is a considerable tension with the latest Planck result at $z=0$: in the case $\alpha_{\rm max}=1$ there is a severe tension whereas for $\alpha_{\rm max}=2$ there is a $3\sigma$ tension, which is compatible with the current $H_{0}$ tension. For all cases the error range increases when $z\lesssim 0.25$, which means that a possible interaction becoming dynamically relevant at recent times may be a viable possibility.

Still in the LM approach, all the results with SNe Ia data are compatible with a constant value and are in agreement with the latest Planck result. As for the CC result, for $z<0.25$ the result are degenerate. Lastly, for the angular BAO data, the result obtained with model-independent priors is clearly inconclusive since, for all values of $z$, within a $3\sigma$-range all interval $\left[0,1\right]$ is admissible for $r_{0}\left(z\right)$. This degenerate result is related to the big error in the model-independent determination of $r_{s}$. However, using Planck priors, there is a disagreement signature when $z\approx 0.3$. For the GP analysis, all the results are consistent with a constant value and also are in agreement with the $\Lambda$CDM Planck result.

Acknowledgements.

RvM acknowledges support from CNPq and the Federal Commission for Scholarships for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS No. 2018.0443) for the academic year 2018-2019. VM thanks CNPq and FAPES for partial financial support. JEG acknowledges support from CNPq (PDJ No. 155134/2018-3). JSA acknowledges support from CNPq (Grants no. 310790/2014-0 and 400471/2014-0) and FAPERJ (Grant no. 204282).

Appendix A Calibrated MSPE

In order to obtain a performance estimator with the same expectation value of $\chi^{2}_{\nu}$, let us rewrite equation (41) as:

[TABLE]

where $t_{i}$ is the true value of $t$ computed in $z_{i}$, and in the second line we have omitted cross-product terms whose expectation value is zero because of the independence between the data used to fit the model $d_{i}$ and the $\tilde{d}_{i}$ data.777Note that, if $d$ and $\tilde{d}$ are partitions of a correlated dataset, one is neglecting the correlation between these two partitions. The expectation value of the first term is clearly $\tilde{N}$. The expectation value of the second term is not trivial:

[TABLE]

In the following we will use the notation that $\bar{t}_{i}$ is a random variable when inside expectation values while it is $t(z_{i},\{c_{\alpha,{\rm bf}}\})$ when outside, that is, we use the best-fit model in order to estimate the true model $t_{i}$.

The first term in the last equation is the variance of $\bar{t}_{i}$ which can be computed, as in Section V.2, through a change of variables using the covariance matrix on the parameters $\Sigma_{\alpha\beta}$ obtained from equation (30) using the training set:

[TABLE]

where we have used equation (28). Note that, thanks to the linearity in the parameters, the best fits $\{c_{\alpha,{\rm bf}}\}$ were not used and this computation is exact rather than only valid at the first order in a Taylor expansion.

The second term is,

[TABLE]

Using again equation (28) one finds:

[TABLE]

where we used the best-fit parameters $\{c_{\alpha,{\rm bf}}\}$ in order to estimate the true values of the parameters. Combining the equations (45), (46) and (48), the equation (44) can be rewritten as:

[TABLE]

where:

[TABLE]

It is then straightforward to obtain that:

[TABLE]

where $\tilde{\Sigma}_{\alpha\beta}$ is the covariance matrix on the parameters obtained from equation (30) using the validation set, that was not used to fit the model.

Motivated by this results, we propose a new generalization of the MSPE:

[TABLE]

whose expectation value is unity. From the previous equation it follows that for large sets the correction $\delta/\tilde{N}$ should be negligible. Indeed, as training and validation sets have usually sizes of the same order of magnitude, one has $\delta\approx\alpha_{\rm max}+1$.

Appendix B learning_curve package

All the learning curves presented in this work were obtained using the package learning_curve. The package learning_curve consists of three python scripts to compute and plot learning curves. The three python scripts are the following:

• •

learning_curve.py: general parallelized script for computing learning curves for any linear template function.

• •

learning_curve_linear.py: script for computing and plotting learning curves for some specific template functions (polynomial, log, inverse or square).

• •

plot.py: script for plotting the learning curves from the output files obtained using the script learning_curve.py.

In this work, only the learning_curve_linear.py script was used, because its template functions coincide with the template functions adopted in the present analysis. The package is available for download at github.com/rodrigovonmarttens/learning_curve.