Changing the prior: absolute neutrino mass constraints in nonlocal gravity

TL;DR

This study investigates how extending the neutrino sector in a nonlocal gravity model can resolve tensions with cosmological data, leading to non-zero neutrino masses and improved parameter fits.

Contribution

It introduces a neutrino extension in a nonlocal gravity model, demonstrating how it alleviates tensions and aligns with observational data, including neutrino mass constraints.

Findings

Neutrino extension resolves the $H_0$ tension.

Nonlocal gravity model with neutrinos fits data well.

Neutrino mass inferred to be around 0.21 eV at 2 sigma.

Abstract

Prior change is discussed in observational constraints studies of nonlocally modified gravity. In the latter, a model characterized by a modification of the form $\sim m^2 R\Box^{-2}R$ to the Einstein-Hilbert action was compared against the base $\Lambda$CDM one in a Bayesian way. It was found that the competing modified gravity model is significantly disfavored (at $22 \,$:$\, 1$ in terms of betting-odds) against $\Lambda$CDM given CMB+SNIa+BAO data, because of a dominant tension appearing in the $H_0 \,$-$\, \Omega_M$ plan. We identify the underlying mechanism generating such a tension and show that it is mostly caused by the late-time, quite smooth, phantom nature of the effective dark energy described by the nonlocal model. We find possible solutions for it to be resolved and explore a given one that consists in extending the initial baseline from one massive neutrino eigenstate to three degenerate ones, whose absolute mass $\sum m_\nu \, / \, 3$ is allowed to take values within a reasonable prior interval. As a net effect, the absolute neutrino mass is inferred to be non-vanishing at $2 \sigma$ level, best-fitting at $\sum m_\nu \approx 0.21 {\, \rm eV}$, and the Bayesian tension disappears rendering the nonlocal gravity model statistically equivalent to $\Lambda$CDM, given recent CMB+SNIa+BAO data. We also discuss constraints from growth rate measurements $f \sigma_8$ whose fit is found to be improved by a larger massive neutrino fraction as well. The $\nu$-extended nonlocal model also prefers a higher value of $H_0$ than $\Lambda$CDM, therefore in better agreement with local measurements.

Tables

Table 1. Table 1: Summary of the means, standard deviations and (effective) χ 2 superscript 𝜒 2 \chi^{2} goodness-of-fit values for the one-dimensional marginalized likelihood distributions of the Λ ​ CDM Λ CDM \Lambda{\rm CDM} , ν ​ Λ ​ CDM 𝜈 Λ CDM \nu\Lambda{\rm CDM} , R ​ R 𝑅 𝑅 RR and ν ​ R ​ R 𝜈 𝑅 𝑅 \nu RR models obtained with the Planck dataset. The Δ ​ χ 2 Δ superscript 𝜒 2 \Delta\chi^{2} values are taken with respect to the Λ ​ CDM Λ CDM \Lambda{\rm CDM} χ 2 superscript 𝜒 2 \chi^{2} values for each dataset, where χ 2 ≡ − 2 ​ ln ​ ℒ superscript 𝜒 2 2 ln ℒ \chi^{2}\equiv-2\,{\rm ln}\mathcal{L} , with ℒ ℒ \mathcal{L} being the likelihood function. All bounds shown correspond to 1 ​ σ 1 𝜎 1\sigma unless explicitly stated otherwise.

	$Planck$	$Planck$	$Planck$	$Planck$
	$\Lambda \mathrm{CDM}$	$\nu \Lambda \mathrm{CDM}$	$RR$	$\nu RR$
$100{\omega }_b$	${2.225}_{-0.016}^{+0.016}$	${2.220}_{-0.017}^{+0.017}$	${2.227}_{-0.016}^{+0.016}$	${2.222}_{-0.017}^{+0.017}$
${\omega }_{cdm}$	${0.1194}_{-0.0015}^{+0.0014}$	${0.1198}_{-0.0016}^{+0.0015}$	${0.1191}_{-0.0015}^{+0.0014}$	${0.1196}_{-0.0016}^{+0.0015}$
$H_0$	${67.50}_{-0.66}^{+0.65}$	${66.12}_{-1.2}^{+2.1}$	${71.51}_{-0.84}^{+0.81}$	${69.57}_{-1.6}^{+2.5}$
$\ln \left({10}^{10}{A}_s\right)$	${3.064}_{-0.025}^{+0.025}$	${3.080}_{-0.034}^{+0.030}$	${3.047}_{-0.025}^{+0.026}$	${3.071}_{-0.035}^{+0.032}$
$n_s$	${0.9647}_{-0.0049}^{+0.0048}$	${0.9637}_{-0.0050}^{+0.0050}$	${0.9649}_{-0.0049}^{+0.0049}$	${0.9639}_{-0.0052}^{+0.0051}$
$\tau$	${0.06530}_{-0.014}^{+0.014}$	${0.07312}_{-0.018}^{+0.016}$	${0.05733}_{-0.014}^{+0.014}$	${0.06905}_{-0.018}^{+0.017}$
$\sum m_{\nu }\left[\mathrm{eV}\right]$	$0.06\left(\mathrm{fixed}\right)$		$0.06\left(\mathrm{fixed}\right)$
${\sigma }_8$	${0.8171}_{-0.0089}^{+0.0089}$	${0.7949}_{-0.016}^{+0.033}$	${0.8487}_{-0.0096}^{+0.0097}$	${0.8212}_{-0.020}^{+0.038}$
$\Delta {\chi }_{Planck}^2$	$0\left({\chi }^2=12943.30\right)$	$-0.04$	$-1.6$	$-1.6$

Table 2. Table 2: As for Table I but obtained with the BAPJ dataset.

	$BAPJ$	$BAPJ$	$BAPJ$	$BAPJ$
	$\Lambda \mathrm{CDM}$	$\nu \Lambda \mathrm{CDM}$	$RR$	$\nu RR$
$100{\omega }_b$	${2.228}_{-0.015}^{+0.014}$	${2.229}_{-0.015}^{+0.014}$	${2.213}_{-0.015}^{+0.014}$	${2.221}_{-0.015}^{+0.014}$
${\omega }_{cdm}$	${0.1190}_{-0.0011}^{+0.0011}$	${0.1189}_{-0.0011}^{+0.0011}$	${0.1210}_{-0.0010}^{+0.0010}$	${0.1197}_{-0.0012}^{+0.0012}$
$H_0$	${67.67}_{-0.50}^{+0.47}$	${67.60}_{-0.55}^{+0.66}$	${70.44}_{-0.56}^{+0.56}$	${69.49}_{-0.80}^{+0.79}$
$\ln \left({10}^{10}{A}_s\right)$	${3.066}_{-0.026}^{+0.019}$	${3.071}_{-0.029}^{+0.026}$	${3.027}_{-0.023}^{+0.027}$	${3.071}_{-0.032}^{+0.032}$
$n_s$	${0.9656}_{-0.0043}^{+0.0041}$	${0.9661}_{-0.0043}^{+0.0043}$	${0.9601}_{-0.0039}^{+0.0040}$	${0.9635}_{-0.0045}^{+0.0043}$
$\tau$	${0.06678}_{-0.013}^{+0.011}$	${0.06965}_{-0.015}^{+0.014}$	${0.04516}_{-0.012}^{+0.014}$	${0.06880}_{-0.017}^{+0.017}$
$\sum m_{\nu }\left[\mathrm{eV}\right]$	$0.06\left(\mathrm{fixed}\right)$		$0.06\left(\mathrm{fixed}\right)$	${0.219}_{-0.084}^{+0.083}$
${\sigma }_8$	${0.8170}_{-0.0095}^{+0.0076}$	${0.8157}_{-0.011}^{+0.013}$	${0.8443}_{-0.0099}^{+0.010}$	${0.8215}_{-0.017}^{+0.017}$
$\Delta {\chi }_{Planck}^2$	$0\left({\chi }^2=12943.42\right)$	$-0.14$	$-0.14$	$-1.52$
$\Delta {\chi }_{BAO}^2$	$0\left({\chi }^2=4.42\right)$	$0$	$2.48$	$2.38$
$\Delta {\chi }_{JLA}^2$	$0\left({\chi }^2=683.2\right)$	$-0.12$	$3.56$	$2.5$
$\Delta {\chi }_{total}^2$	$0\left({\chi }^2=13631.04\right)$	$-0.26$	$5.9$	$3.36$

Table 3. Table 3: Scale used for comparing model ℳ 1 subscript ℳ 1 \mathcal{M}_{1} against model ℳ 0 subscript ℳ 0 \mathcal{M}_{0} in this article, i.e. for interpetating their BIC difference Δ ​ χ 01 2 ≡ χ 1 2 − χ 0 2 Δ subscript superscript 𝜒 2 01 subscript superscript 𝜒 2 1 subscript superscript 𝜒 2 0 \Delta\chi^{2}_{01}\equiv\chi^{2}_{1}-\chi^{2}_{0} and their log-Bayes factors ln ⁡ B 01 subscript 𝐵 01 \ln B_{01} . Positivity of the latter tends to favor ℳ 0 subscript ℳ 0 \mathcal{M}_{0} . These scales are taken as a rule of thumb inspired by Secs. 2.6-2.10 of Ref.  Burnham and Anderson ( 2002 ) in accordance with the (more conservative) Jeffreys’ scale of Ref.  Trotta ( 2008 ) (see also Ref.  Efstathiou ( 2008 ) for a comparison of the latter with the original scale proposed by Jeffrey).

	Interpretation	$\Delta {\chi }^2$	$\ln B_{01}$
	“inconclusive”	$0$–$\mathrm{\hspace{0.17em}2}$	$0$–$\mathrm{\hspace{0.17em}1}$
	“weak”	$2$–$\mathrm{\hspace{0.17em}6}$	$1$–$\mathrm{\hspace{0.17em}2.5}$
	“moderate-to-strong”	$6$–$\mathrm{\hspace{0.17em}10}$	$2.5$–$\mathrm{\hspace{0.17em}5}$
	“strong”	$>10$	$>5$

Equations

$$S_{RR}=\frac{1}{16\pi G}\int {\mathrm{d}}^{4}x\sqrt{-g}\left[R-\frac{m^2}{6}R{\square }^{-2}R-{\mathcal{L}}_m\right],$$

$$S_{RR,\mathrm{loc}}$$

$$\square U$$

$$\square V$$

$$U$$

$$V$$

$$G_{\mu \nu }-\frac{m^2}{6}{K}_{\mu \nu }=8\pi G{T}_{\mu \nu },$$

$$\square U=-R,$$

$$\square V=-U,$$

$$K_{\mu \nu }\equiv 2V{G}_{\mu \nu }-2{\nabla }_{\mu }{\nabla }_{\nu }V-2{\nabla }_{(\mu }U{\nabla }_{\nu )}V$$

$$+\left(2\square V+{\nabla }_{\alpha }U{\nabla }^{\alpha }V-\frac{U^2}{2}\right)g_{\mu \nu },$$

$$\displaystyle\theta_{\rm base}=\big{(}H_{0},100\,\omega_{b},\omega_{cdm},{\rm ln}(10^{10}A_{s}),n_{s},\tau\big{)},$$

$$\displaystyle P\big{(}M_{\nu}|\Lambda{\rm CDM})=\delta\big{(}M_{\nu}/{\rm eV}-0.06\big{)}\,,$$

$$\displaystyle\theta_{\rm\nu}=\big{(}H_{0},100\,\omega_{b},\omega_{cdm},{\rm ln}(10^{10}A_{s}),n_{s},\tau,\Sigma\,m_{\nu}\big{)}\,,~{}~{}~{}$$

$$\displaystyle H(z)=H_{0}\big{[}\Omega(z)+\Omega_{\rm DE}(z)\big{]}^{1/2}\,,$$

$${\Omega }_{\mathrm{DE}}(z)$$

$${\Omega }_{\mathrm{DE}}(z\approx 0)$$

$$r_{\ast }$$

$$D_A(z_{\ast })$$

$$c_s=1/\sqrt{3\left[1+3{\Omega }_b/(4{\Omega }_{\gamma })\right]},$$

$$D_L(z)\equiv (1+z){\int }_0^z\frac{\mathrm{d}{z}^{\prime }}{H(z^{\prime })},$$

$$\mathrm{BIC}\equiv {\chi }^2+k\ln N,$$

$$B_{\nu \Lambda ,\nu RR}$$

$$=$$

$$\displaystyle B_{i,\nu i}=\frac{P\big{(}\sum m_{\nu}=0.06\,\big{|}\,d,\mathcal{M}_{\nu i}\big{)}}{P\big{(}\sum m_{\nu}=0.06\,\big{|}\,\mathcal{M}_{\nu i}\big{)}}\,,$$

$$B_{\nu \Lambda ,\nu RR}$$

$$\mathrm{BIC}:$$

$$\mathrm{Bayes}:$$

$$\displaystyle B_{i,\nu i}=\frac{P\big{(}\sum m_{\nu}=0.06\,\big{|}\,d,\mathcal{M}_{\nu i}\big{)}}{P\big{(}\sum m_{\nu}=0.06\,\big{|}\,\mathcal{M}_{\nu i}\big{)}}=\frac{L_{i}}{V_{i}}P_{i}\,,$$

$$L_{RR}$$

$$L_{\Lambda }$$

$$B_{\Lambda ,\nu \Lambda }=12.5\times 5=62.5,B_{RR,\nu RR}=1\times 5,$$

$${\chi }_{\Lambda \mathrm{CDM}}^2=3.9,{\chi }_{RR}^2=6.5,{\chi }_{\nu RR}^2=5.2,$$

Full text

I Introduction

Modern cosmology has undergone fast developments over the past two decades. The evidence for accelerated expansion from observations of distant Type Ia supernovae Riess et al. (1998); Perlmutter et al. (1999) led to the abandon of the inflationary cold dark matter (CDM) paradigm and to the birth of the $\Lambda{\rm CDM}$ one, introducing a cosmological constant $\Lambda$ into the model. Such an extension adds an overwhelming, thereby constant, dark energy component to the corresponding homogeneous and isotropic Friedmann-Lemaître-Robertson-Walker (FLRW) universe and at the same time raises fundamental theoretical questions about its origin, late-time domination and naturalness Carroll (2003); Bianchi and Rovelli (2010). The idea of modelling differently such a component for addressing these issues keeps to be actively explored nowadays through different perspectives, in particular in attempting to modify General Relativity (GR) (see e.g. Nojiri and Odintsov (2011); Clifton et al. (2012); Joyce et al. (2015); Berti et al. (2015); Koyama (2016) for recent reviews). Among these, nonlocal modifications to GR have attracted much interest over the years. They are motivated by radiative corrections induced by gravity or light-matter components on curved spacetime Barvinsky and Vilkovisky (1987, 1990); Gorbar and Shapiro (2003), also responsible for the trace anomaly Antoniadis and Mottola (1992); Barvinsky et al. (1995) or non-perturbative effects giving raise to a scale dependent Newton’s constant Arkani-Hamed et al. (2002); Barvinsky (2003); Hamber and Williams (2005); Lopez Nacir and Mazzitelli (2007); Barvinsky (2012a), but also by renormalizability and singularity issues in higher-derivative theories of gravity Biswas et al. (2012); Modesto (2012) or string inspired scenarios Biswas et al. (2006); Calcagni and Modesto (2015). Models have also been proposed in a bottom-up approach for appreciating the effects of various types of nonlocal operators within the cosmological context Wetterich (1998); Espriu et al. (2005); Nojiri and Odintsov (2008); Deser and Woodard (2007); Koivisto (2008); Deser and Woodard (2013); Maggiore (2014); Ferreira and Maroto (2013); Maggiore and Mancarella (2014); Mitsou (2015); Cusin et al. (2016a, b); Nersisyan et al. (2017); Vardanyan et al. (2017).

Two nonlocal models of this sort have been recently proposed, one of which is characterized by the addition of a term \sim m^{2}\big{(}g_{\mu\nu}\Box^{-1}R\big{)}^{T} to Einstein’s equations Maggiore (2014), whereas a second one modifies the Einstein-Hilbert action by a term $\sim m^{2}R\Box^{-2}R$ Maggiore and Mancarella (2014). Their effective dark energy phenomenology has been studied in Ref. Dirian et al. (2014) where it was found that both models describe a quite smooth (i.e. whose perturbations are small), phantom dark energy component emerging at late-times, with an equation-of-state today of $w_{\rm DE}=-1.04$ and $w_{\rm DE}=-1.15$ respectively. Furthermore, the models feature a fifth force that enhances the clustering of linear structures compared to that in $\Lambda{\rm CDM}$: at $\lesssim 6\%$ level in the linear matter power spectrum around the BAO scale (see also Ref. Barreira et al. (2014a) where nonlinear structure formation through N-body simulation has been studied for the model of Ref. Maggiore and Mancarella (2014)). Observational constraints using Bayesian techniques were then carried out in Refs. Dirian et al. (2015, 2016) in a complementary perspective Tegmark et al. (1998); Eisenstein et al. (1999), i.e. using cosmological data from Cosmic Microwave Background (CMB), Baryon Acoustic Oscillations (BAO) and Type Ia supernova (SNIa) observations. Constraints on the linear growth rate of structures were also obtained a posteriori, that is, constraining the quantity $f\sigma_{8}$ derived from both nonlocal models on their respective CMB+SNIa+BAO bestfit with Redshift-Space Distortions (RSD) measurements. The joined CMB+SNIa+BAO constraints showed that, provided a prior parametrization fixed on the so-called Planck 2015 baseline Ade et al. (2015a), the model of Ref. Maggiore (2014) is indistinguishable from standard $\Lambda{\rm CDM}$ with a Bayes factor of $1.0$, whereas the one presented in Ref. Maggiore and Mancarella (2014) is significantly disfavored with a Bayes factor of $22.7$. The latter discrepancy was shown to result from a severe CMB-SNIa tension appearing in the $H_{0}\,$–$\,\Omega_{M}$ plane, breaking the nonlocal models’ concordance given the data.

In this article, we analyze in more details the origin of this tension and find a solution for resolving it. In particular, we will see that changing the neutrino sector of the aforementioned baseline from one dominant active mass-eigenstate to three degenerated ones, whose absolute mass is taken as a free parameter, restores the concordance of the nonlocal model to a non-negligible extent. Effectively, such a resolution exploits degeneracies between modified gravity effects of the nonlocal model and those caused by a more massive neutrino component. Similar degeneracies have already been noticed in local modified gravity theories, for instance at linear level in TeVeS Skordis et al. (2006), covariant galileons Barreira et al. (2014b), $K$-mouflage Barreira et al. (2015) and recently in Horndeski models Bellomo et al. (2017), but also at the nonlinear one through N-body simulations of $f(R)$ scenarios in Ref. Baldi et al. (2014).

The rest of this paper is organized as follows. The specific nonlocal gravity model we consider will briefly be reviewed in section II . In section III , we expose how the latter is embedded into a statistical model, given specific data that we also present. Furthermore, we specify our prior parametrization and justify such a choice from empirical evidences. Section IV presents the results drawn from observational constraints, attempts for a comprehensive analysis of the resolution of the aforementioned tension and quantifies it using Bayesian model comparison methods. We present a summary of our work in section V before we conclude.

II The $RR$ nonlocal gravity model

In this section, we briefly introduce the nonlocal gravity model originally proposed by Ref. Maggiore and Mancarella (2014) which, following Ref. Dirian et al. (2016) for conveniency, will be referred to as the $RR$ model. Our conventions, notations and strategy to solve the equations are similar to that of Refs. Dirian et al. (2016, 2014), that the interested reader is invited to consult for more details.

The $RR$ model is defined through an extension of the Einstein-Hilbert action reading,

[TABLE]

where $\mathcal{L}_{m}$ is the Lagrange density of minimally coupled matter fields and $\Box^{-1}$ is a formal notation for a Green’s function of the curved-space d’Alembert operator $\Box\equiv\nabla^{\mu}\nabla_{\mu}$. Since the kernel of $\Box$ is non-trivial such a Green’s function is not unique and its precise structure in Eq. (1) is not known, as well as the origin of the mass scale $m$. This originates from the fact that this model was built through a bottom-up approach and its embedding into a more fundamental framework still remains under investigation Maggiore (2015, 2016a), see also Ref. Maggiore (2016b) for a recent review. From a more technical point of view, several subtleties arise from the presence of $\Box^{-1}$ into an action. One of them concerns the causal character of the classical evolution governed by the corresponding Euler-Lagrange equations and has been discussed in details in the literature (see e.g. Refs. Soussa and Woodard (2003); Barvinsky (2003)). Another one ties to the “localized” version of the action (1), most conveniently used for computing the corresponding cosmology without having to deal with less manipulable integro-differential systems. The present work makes use of such a local version, taken to be111See also Ref. Dirian and Mitsou (2014) where another localization has been proposed.

[TABLE]

where $U$ and $V$ are two “auxiliary” scalar fields and $\xi_{1,2}$ are Lagrange multipliers enforcing the constraints,

[TABLE]

sourcing the latter. Invoking a given (left) inverse, one can solve the latter formally in simply writing,

[TABLE]

which allows one to integrate $U$, $V$ out from the action (and $\xi_{1,2}$ in a similar way), leading back to the nonlocal action (1). The same procedure can be applied at the level of the equations of motion, that are

[TABLE]

with

[TABLE]

where $T^{\mu\nu}\equiv\left(2/\sqrt{-g}\right)\delta\left(\mathcal{L}_{m}\sqrt{-g}\right)/\delta g_{\mu\nu}$ is the matter source. In writing Eqs. (5) and (6) at the level of the equations of motion, causality of the classical dynamics requires that the used source-convolving kernels are all of the retarded kind and, in accordance with their artificial nature, that the auxiliary fields possess no homogeneous (i.e. free) solutions Koshelev (2009); Koivisto (2010); Barvinsky (2012b); Deser and Woodard (2013); Maggiore (2014); Foffa et al. (2014a, b); Maggiore and Mancarella (2014); Dirian and Mitsou (2014), at least into the in state. This uniquely fixes the inverse of $\Box$ used in integrating out $U$ and $V$ and catches the classical solutions that we are interested in in that work. At the “localized” level, trivial homogeneous modes for the auxiliary fields is realized in fixing vanishing values for $U$ and $V$ and their time derivatives on the initial hypersurface. Here we interpret these requirements as theory-level data that one should supplement to the local action (2) from the beginning, and changing these prescriptions changes the underlying nonlocal model.

According to the procedure used until now Dirian et al. (2014, 2015, 2016), we set the initial hypersurface deep into the radiation dominated era (RD). During this period, at the background level, the curvature scalar is sub-dominating compared to the overall energy scale of the process, so one can naively set the space-averaged value $\left.\bar{R}\right|_{\mathrm{RD}}\simeq 0$ and see that the resulting $\left.\bar{U},\bar{V}\right|_{\mathrm{RD}}$ are therefore not sourced, and remain small. Linear perturbations $\left.\delta U,\delta V\right|_{\mathrm{RD}}$ induced by gravitational ones remain also small Dirian and Mitsou (2014). However, once matter starts to dominate over radiation, the latter quantities acquire a non-trivial dynamics leading to the emergence of a late time, quite smooth dynamical dark energy component Dirian et al. (2014) driving the accelerated expansion of the Universe which ends up into a Big Rip Dirian and Mitsou (2014); Caldwell et al. (2003). For being ultimately legitimate, this choice of initial conditions needs to assume that the secular growth of $U$ and $V$ during earlier stages of the Universe is mild enough for not strongly affecting the vanishing values chosen while starting in RD. For the case of inflationary scenarios, it can be argued that the additional term in Eq. (1) is suppressed by the corresponding inflationary energy scale, since it only modifies the theory into the far-infrared. In any case this question deserves a special attention and a more quantitative analysis is needed (see Refs. Maggiore (2015); Cusin et al. (2016b); Codello and Jain (2016) where such a scenario was considered). Of course, the configuration reached by the auxiliary fields in RD will depend on the particular assumed model of inflation but also on other processes affecting the scalar curvature within the primordial Universe, such as the electroweak or QCD phase transition, conformal anomalies, or during RD itself through the presence of thermalized Standard Model massive particles Caldwell and Gubser (2013). This issue was recently anticipated in Ref. Nersisyan et al. (2016), where the authors studied the effect of varying the auxiliary fields’ initial conditions deep into RD over a broad range of values, unveiling in particular the existence of another phenomenologically viable cosmology of the nonlocal model (1).

We close this description in emphasizing that, from a statistical perspective which will be taken in the following, fixing vanishing initial conditions for the auxiliary fields is a modelling assumption or theoretical prior (much as pretending that $\Lambda$-dark energy is replaced by the $RR$ one in our context). Therefore, the only remaining free parameter added to the Einstein-Hilbert action is the mass scale $m$, whose value is to be fixed from spatial flatness condition – similarly of fixing the value of $\Lambda$ in flat $\Lambda{\rm CDM}$ models. The model has therefore the same number of free parameters as the $\Lambda{\rm CDM}$ model.

III Cosmological models

The cosmological models incorporating $\Lambda$- and $RR$-dark energies considered in this article are one-parameter extensions of those studied in Ref. Dirian et al. (2016) and are similarly tested within a Bayesian framework.

In this section, we therefore only provide a brief, self-contained review of the datasets and parametrization used and refer the reader to Ref. Dirian et al. (2016), and references therein, for further details about their structure and construction. Technical details relative to the equations and numerical implementation used can be found in Sec. 2 and App. A. of Ref. Dirian et al. (2016).

III.1 Datasets

For performing our global fit we utilize Cosmic Microwave Background (CMB) data complemented with distant Type Ia supernovae (SNIa) observations and distance measurements from several Baryon Acoustic Oscillations (BAO) surveys. Together with CMB lensing data, SNIa and BAO observations allow to apply constraints on the late time expansion of the Universe since they break further primary CMB degeneracies. We will then discuss a posteriori constraints on the growth rate from Redshift-Space Distortions (RSD) measurements.

The CMB dataset we consider comes from the Planck satellite mission and is made of the lowTEB power spectra data for multipoles $\ell\leq 29$, the high-$\ell$ TT, TE and EE ones for $\ell>29$ Aghanim et al. (2016a), as well as the data relative to the power spectrum of the reconstructed CMB lensing potential Ade et al. (2015b). In the following, this dataset will be referred to as the Planck dataset. For the BAO data, we include the isotropic measurements reported in Ref. Beutler et al. (2011) at $z_{\rm eff}=0.106$ from 6dfGS and Ref. Ross et al. (2015) at $z_{\rm eff}=0.15$ from SDSS-MGS DR7, as well as the anisotropic ones of Ref. Anderson et al. (2014) from the LOWZ ($z_{\rm eff}=0.32$) and CMASS ($z_{\rm eff}=0.57$) samples from the BOSS DR11 release. The SNIa data we use are those of the SDSS-II/SNLS combined analysis (also known as JLA), which comprises $740$ objects at $z\lesssim 1$ Betoule et al. (2014).

Below we study constraints with two dataset combinations: (i) Planck and (ii) the combined Planck, BAO and SNIa data, which we shall refer to as the BAPJ dataset to preserve the nomenclature of Ref. Dirian et al. (2016).

III.2 Parametrization and prior specification

In accordance with the use of the above measurements, the complement of the dark energy sector of the models considered in this paper are defined following the so-called Planck 2015 baseline Ade et al. (2015a). Regarding the radiation/matter content of the Universe, ionization history and primordial initial conditions, the baseline assumes a particular modelling that has now become the standard cosmological setting on which is built the base $\Lambda{\rm CDM}$ model. Predictions from such a model are computed in using a cosmological linear Einstein-Boltzmann code such as CAMB Lewis et al. (2000) or CLASS Blas et al. (2011). In our study, we make use of the latter and provide a modified version of it on GitHub (see git for the link) which includes the $RR$ gravity model introduced in the previous subsection. Bayesian parameter extraction and model selection are carried out with the Markov Chain Monte Carlo (MCMC) code MONTEPYTHON Audren et al. (2013) originally interfaced with CLASS.

The baseline specifies a continuous 6-dimensional parametrization which can be provided by the vector,

[TABLE]

with $H_{0}\equiv 100h~{}{\rm km/s/Mpc}$, the Hubble expansion rate today; $\omega_{i}\equiv\Omega_{i}h^{2}$, where $\Omega_{i}$ is the present energy density fraction of baryons ($i=b$) and cold dark matter ($i=cdm$); $A_{s}$ and $n_{s}$ are the amplitude and tilt of the power spectrum of primordial fluctuations respectively and $\tau$ is the optical depth to reionization. In the following, we will also use derived quantities as $\Omega_{m}$ the total matter density fraction today, $\sigma_{8}$ the root mean square linear matter fluctuations in a sphere of radius $8{\rm\,Mpc/h}$ at $z=0$ and the dark energy density fractions $\Omega_{\Lambda}$, $\Omega_{RR}$ corresponding to the $\Lambda{\rm CDM}$ and $RR$ models respectively. The latter are determined in requiring a vanishing spatial curvature of the Universe $\Omega_{K}=0$ and this is done adequately tuning the parameter controlling the dark energy density within the model specified: $\Lambda$ in $\Lambda{\rm CDM}$ or $m$ in $RR$. In this work we use improper flat prior with edges everywhere unbounded expect for the lower bound of the optical depth $\tau$ taken to be $0.01$, in accordance with Gunn-Peterson trough observations (see e.g. Ref. Becker et al. (2001)).

In the neutrino sector, the baseline assumes one dominant active neutrino species of fixed mass $M_{\nu}=0.06{\,\rm eV}$, i.e. with prior

[TABLE]

whose temperature is tuned for reproducing predictions from neutrino decoupling computations Mangano et al. (2005). The massive neutrino species is supplemented with a massless neutrino one whose effective number of degrees-of-freedom is adjusted so as to reproduce a total effective number of relativistic components of $N_{\rm eff}=3.046$ in the early Universe (see e.g. Refs. Archidiacono et al. (2013); Lesgourgues and Pastor (2014)).

As will be justified more clearly below, the extension we consider in this work assumes three active massive species with degenerated mass-eigenstates, i.e. the sum of their masses respects $\sum m_{\nu}\equiv 3m_{0}$, where the absolute neutrino mass $m_{0}$ is taken to be a free parameter. Their temperature and the additional massless neutrino component are fixed respecting the same conditions as those of the initial baseline. We therefore adopt the following parametrization,

[TABLE]

and the prior interval on the additional parameter is taken to be uniform and compact $\sum m_{\nu}\in[0,5]{\,\rm eV}$, consistent with existing data (see subsection IV.6 for more details and references on the latter choice).

The cosmological models denoted by $\Lambda{\rm CDM}$ and $RR$ are parametrized by $\theta_{\rm base}$ while their extended versions build out of $\theta_{\rm\nu}$ will be called $\nu\Lambda{\rm CDM}$ and $\nu RR$ respectively.

IV Observational constraints and model comparison

In this section, we present and study our observational constraints and model comparison results.

We start by analyzing some structural features of the $RR$ model and, comparing them to the ones of standard $\Lambda{\rm CDM}$, we explain the origin of the tensions first found in Ref. Dirian et al. (2015), given the nonlocal model. Then, we present the parameter constraints on the $\nu RR$ model introduced above and discuss Bayesian model selection through the computation of various Bayes factors, comparing pairs of the four above-mentioned models. We finally also discuss the capabilities of growth rate measurements to constrain further the cosmologies at hand.

IV.1 The origin of the CMB-SNIa tension given $RR$

IV.1.1 $Planck$ dataset

Table 1 summarizes the constraints on the $\Lambda{\rm CDM}$ and $RR$ models obtained with the Planck dataset (second and fourth columns). The parameter shifts between the two models can be understood in comparing the relevant features in both models with each other.

A first noteworthy point is that the differences between the results in the two cosmologies are statistically non-significant ($\lesssim 1\sigma$) for all parameters in $\theta_{\rm base}$, with the exception of the background-related parameter $H_{0}$ which undergoes the most significant shift ($\sim 5\sigma$). Notice that $\sigma_{8}$, derived at the linear level, also undergoes a significant change ($\sim 3\sigma$). As already noticed in Ref. Barreira et al. (2014a), the latter results from the enhanced clustering in the nonlocal model, mostly induced by a lower expansion rate (see top panel of Fig. 1) that reduces the Hubble friction to matter perturbations. This effect is supplemented by, although to a milder level, a higher late-time gravitational strength modelled by a time-dependent effective Newton constant $G_{\rm eff}(z,k)$, exhibiting scale-dependence in the far-infrared Dirian et al. (2016). More clustering also increases the lensing power (see bottom panel of Fig. 1) and in turn smooths out temperature fluctuations more efficiently. Moreover, this requires a smaller primordial amplitude $A_{s}$ which comes together with a delayed reionization epoch, given that the CMB damping tail constrains well the combination $A_{s}\,e^{-2\tau}$ at high-$\ell$ 222Observe that the preference for a lower optical depth to reionization $\tau$ within the nonlocal model compared to $\Lambda{\rm CDM}$ is consistent with the results found from the new analysis of the Planck HFI data Aghanim et al. (2016b)..

As already mentioned above, the Planck dataset constrains $H_{0}$ to be larger in $RR$ ($H_{0}\approx 71.51\pm 0.84\ {\rm km/s/Mpc}$) than in $\Lambda{\rm CDM}$ ($H_{0}\approx 67.50\pm 0.66\ {\rm km/s/Mpc}$). This preference for higher $H_{0}$ a fortiori originates from the late-time emerging, quite smooth and phantom nature of the $RR$ effective dark energy compared to that modelled by a cosmological constant $\Lambda$. In the following we attempt to provide a comprehensive explanation of this fact. The Friedmann equation including a dynamical dark energy component with equation of state $w_{\rm DE}(z)$ reads

[TABLE]

where $\Omega_{\rm DE}(z)$ denotes the dark energy density fraction present in the Universe at redshift $z$ and $\Omega(z)$ includes all the other components, that is, the density fraction of cold dark matter $\Omega_{cdm}(z)$ and of baryons $\Omega_{b}(z)$ in the case of the baseline, but also other ingredients available in extensions of it such as the density fraction of massive neutrinos $\Omega_{\nu}$. From the dark energy conservation equation, $\Omega_{\rm DE}(z)$ can be written in terms of the dark energy equation of state,

[TABLE]

which, at low redshift when the dark energy is dominant, can be approximated by

[TABLE]

where we wrote $w_{\rm DE}(z\approx 0)\simeq-1+\delta w_{0}$, with $|\delta w_{0}|\ll 1$ and constant, which is sensible enough for the present discussion. One can see that for the case of a phantom dark energy, such as the one featured by the $RR$ model for which $\delta w_{0}<0$, the predicted dark energy density $\Omega_{RR}(z)$ is generically smaller than $\Omega_{\Lambda}$ at $z\gtrsim 0$. This explains why $H_{RR}(z)$ is lower compared to $H_{\Lambda}(z)$, fixing the other parameter values. This is illustrated by the green line in the upper panel of Fig. 1, which shows the time evolution of the Hubble rate in the $RR$ model when assuming the cosmological parameters of the best-fitting $\Lambda{\rm CDM}$ model to the Planck dataset. Cosmological constraints on the $RR$ model that are sensitive to $H(z)$ will then generically infer primarily higher values of $H_{0}$, $\Omega_{cdm}$, $\Omega_{b}$, etc, for trying to compensate the change induced by the different dark energy modelling. Such a change clearly depends on the particular data set used to derive the constraints. We explore the case of the Planck data into this subsection whereas we examine the SNIa JLA ones into the next one.

From the point of view of Planck data, modifications to $H(z)$ at low redshift alter the predicted angular acoustic scale $\theta_{*}$, which determines the position of the acoustic peaks of the CMB temperature power spectrum. $\theta_{*}$ is measured with a very good precision by Planck ($\lesssim 0.1\%$ at $1\sigma$ in the case of base $\Lambda{\rm CDM}$) and is robust under cosmology change. It is expressed as $\theta_{*}\equiv r_{*}/D_{A}(z_{*})$, where $r_{*}$ is the sound horizon at the redshift of recombination $z_{*}$ and $D_{A}(z_{*})$ the comoving angular diameter distance to recombination defined as,

[TABLE]

where $c_{s}$ is the sound speed of the primordial plasma

[TABLE]

with $\Omega_{\gamma}$ the photon density fraction today. For fixed cosmological parameter values, $r_{*}$ does not change significantly from $\Lambda{\rm CDM}$ to $RR$ because it is a function of early time background configurations and does not depend on the particular dark energy modelling. At late time however, such a modelling becomes important and the lower expansion rate in the $RR$ model leads to a larger $D_{A}(z_{*})$, which in turn lowers $\theta_{*}$. The lower acoustic scale shifts the CMB temperature power spectrum towards higher multipoles $\ell$, yielding the poor fit to the data seen in the middle panel of Fig. 1.

In the case of the Planck baseline, this discrepancy can be resolved in shifting either the background quantities $H_{0}$, $\omega_{b}$ or $\omega_{cdm}$ (or equivalently $\Omega_{b}$ or $\Omega_{cdm}$) present into the $H(z)$ expression Eq. (14). However, the shape information of the first CMB peaks such as their relative position and their relative height provide strong, model-independent constraints on both $\omega_{b}$ and $\omega_{cdm}$ (see e.g. Refs. Ade et al. (2014); Aghanim et al. (2016c)) and there is therefore only significant room for $H_{0}$ to vary. Consequently, since the dark energy featured by the $RR$ model is phantom, $H_{0}$ is doomed to increase. The blue curves in Fig. 1 show the same as the green ones, but with $H_{0}$ adjusted to $H_{0}=71.31\ {\rm km/s/Mpc}$ so as to yield the same $\theta_{*}$ as in the best-fitting $\Lambda{\rm CDM}$ model to Planck. This yields a cosmological scenario that is very similar to the best-fitting $RR$ model to Planck (red curves in Fig. 1), whose goodness-of-fit to Planck is better than the base $\Lambda{\rm CDM}$ with $\Delta\chi^{2}=1.6$. This is mostly because of the lower power in the low-$\ell$ part of the CMB temperature power spectrum, induced by a smaller ISW effect dominating at large-scales. Such a preference being “inconclusive” according to the classification reported in our Table 3, both models are therefore statistically equivalent given Planck 2015 CMB data.

IV.1.2 BAPJ dataset

The agreement with observations of the $RR$ model however degrades when it is confronted against the BAPJ dataset. The observational tensions that arise when one includes the BAO and SNIa data in the analyzes are better illustrated in Fig. 2. The figure shows the 2d marginalized constraints in the $H_{0}\,$–$\,\Omega_{m}$ plane for $\Lambda{\rm CDM}$ (upper panel) and $RR$ (lower panel), obtained individually using the Planck dataset (red), SNIa data (grey) and BAO data (green). Contrary to $\Lambda{\rm CDM}$, for the $RR$ model the marginalized posterior suggests a $\sim 3\,$–$\,4\sigma$ level tension between Planck and SNIa data. According to the above discussion, this can be understood in looking at the luminosity distance relevant for SNIa lightcurves,

[TABLE]

where the expression for $H(z)$ is found in Eq. (14). SNIa measurements only constrain the total matter density $\Omega_{m}$, whereas $H_{0}$ has been integrated out via marginalization on the absolute magnitude. Here the fact that the dark energy in the $RR$ model is on the phantom side has the net effect of raising $\Omega_{m}$ towards higher values in $RR$ than in $\Lambda{\rm CDM}$, for fixed luminosity distance. This shift has already been reported in Ref. Dirian et al. (2014) and is illustrated by the grey contours in Fig. 2. We find $\left.\Omega_{m}\right|_{\Lambda{\rm CDM}}=0.298\pm 0.035$ and $\left.\Omega_{m}\right|_{RR}=0.343\pm 0.033$, exhibiting a $\sim 1\sigma$ shift between the two models, given SNIa JLA data.

This trend is inconsistent with the one inferred from Planck: this explains the origin of the tension quantified at $\Delta\chi^{2}=5.9$, which appears when constraining the $RR$ model using joined Planck and SNIa JLA data.

IV.2 Changing the prior: From $RR$ to $\nu RR$

From the above discussion we have deduced that the late-time phantom nature of the $RR$ effective dark energy induces an increase in $H_{0}$ given Planck data, as this helps to resolve the mismatch with the CMB peaks position constraining $\theta_{*}$. However, the same fact also induces an increase in $\Omega_{m}$ given SNIa data, which is inconsistent with Planck’s preference since the latter provides tight constraints on $\sim\Omega_{m}h^{2}$ which, together with an increase in $H_{0}$, forces $\Omega_{m}$ to go down333This behavior is a generic one for reasonably smooth phantom dark energy, as indicated by the degenerate directions in the $\Omega_{m}\!$ – $\!w_{0}$ plane while constraining ${w_{0}\rm CDM}$ models given equivalent CMB and SNIa data as those used in this work (see e.g. Fig. 16 of Ref. Betoule et al. (2014)).. This results in an overall dominant CMB-SNIa tension that renders the $RR$ model non-concordant and disfavors it with respect to $\Lambda{\rm CDM}$ by $\Delta\chi^{2}=5.9$, given BAPJ data. Furthermore, a posteriori constraints from RSD data have also been studied in Ref. Dirian et al. (2016) and they increase the overall tension even more, up to $\left.\Delta\chi^{2}\right|_{\rm bapj}+\left.\Delta\chi^{2}\right|_{\rm rsd}^{\rm post}=8.5$ ($\sim 3\sigma$) compared to $\Lambda{\rm CDM}$, ruling out the nonlocal model given BAPJ+(post)RSD data.

This problem can be solved by considering extensions of the initial model, that is, not only changing the dark energy parametrization as in Ref. Dirian et al. (2016), but also allowing other physically relevant parameters, otherwise fixed, to vary. Adding such new components will open new possibilities in the global parameter space and possibly provide an access to a new global maximum of the posterior probability distribution. The consequence of such a procedure is however the introduction of new degeneracies in the extended cosmological model coming together with a loss of constraining power for fixed data combination.

As already mentioned here-above, in this work we allow the sum of the neutrino masses to vary. This is justified for three reasons. First, constraints on $\sum m_{\nu}$ coming from terrestrial experiments are very weak, therefore there no obvious reason to fix $\sum m_{\nu}=0.06\,{\rm eV}$ on empirical grounds. Indeed, this value only corresponds to the smallest mass-splitting measured by oscillations experiments (see subsection IV.6 for more details and references). Second, an increase in neutrino masses would not alter the expansion rate at early times if the neutrinos are still relativistic at photon decoupling, so that CMB anisotropies remain unaffected, and raise the energy density of pressureless matter (after they turn non-relativistic), which therefore increases the expansion rate during the matter dominated era. Third, the free-streaming behavior exhibited by a more massive neutrino component helps to tame the growth rate of structures and therefore potentially lowers the additional discrepancy caused by the inclusion of RSD data (see e.g. Refs. Archidiacono et al. (2017); Vagnozzi et al. (2017) for analyses of the effects induced by massive neutrinos on cosmological observables). Thus, we can expect that the $RR$ model would prefer a higher value of $\sum m_{\nu}$ than $0.06\,{\rm eV}$, with a corresponding decrease of the CMB-SNIa tension. This is what is discussed in more details in the following.

The cyan curves in Fig. 1 illustrate such facts in showing the same as the blue curves, but instead of adjusting $H_{0}$ to give the same $\theta_{*}$ as in $\Lambda{\rm CDM}$, one adjusts the neutrino masses to $\sum m_{\nu}=0.42\ {\rm eV}$. The middle panel of Fig. 1 also confirms that this would help to drastically improve the goodness-of-fit to Planck data compared to the $RR$ model when using the $\Lambda{\rm CDM}$ best-fitting parameters. This is a clear sign that the degenerate effects of $H_{0}$, $\sum m_{\nu}$ and a late-time phantom dark energy on $\theta_{*}$ can therefore be exploited to try to reconcile the Planck and SNIa JLA constraints given the $RR$ nonlocal gravity.

Figure 2 illustrates better the beneficial impact of varying $\sum m_{\nu}$ in the constraints given the $RR$ model compared to the $\Lambda{\rm CDM}$ one. The blue contours show the Planck constraints on the $H_{0}\!$ – $\!\Omega_{m}$ plane when $\sum m_{\nu}$ is a free parameter (quoted Planck (free $\sum m_{\nu}$) for definiteness). For the nonlocal gravity model, the Planck contour is now overlapping the SNIa one. This illustrates the fact that allowing $\sum m_{\nu}$ to vary weakens the CMB-SNIa tension, as quantified by the corresponding individual $\Delta\chi^{2}$ values reported in Table 2.

A remarkable aspect of the combination of the Planck, SNIa and BAO data in the constraints given $RR$ nonlocal gravity is the evidence for non-vanishing neutrino masses. Figure 3 shows that $\sum m_{\nu}>0$ at $\sim 2\sigma$ level, with the best-fit value $\sum m_{\nu}\approx 0.21\ {\rm eV}$. As depicted above, such a shift is primarily caused by the relatively smooth, late-time and phantom nature featured by the effective dark energy described by the $RR$ nonlocal model. These constraints are very different than in $\nu\Lambda{\rm CDM}$ for which the data only sets an upper bound on $\sum m_{\nu}$. In the following, we will see that such a preference of the $\nu\Lambda{\rm CDM}$ model for lower values of $\sum m_{\nu}$ reflects one of its weakness in a Bayesian model comparison context, and therefore opens room for alternative dark energy models including similarly a varying $\sum m_{\nu}$ to compete with it.

IV.3 Bayesian model comparison

In the following, we compare the $\nu\Lambda{\rm CDM}$ and $\nu RR$ models given BAPJ data computing the associated Bayes factor $B_{\nu\Lambda,\nu RR}$ and set it side by side with the Bayesian Information Criterion (BIC)444The Bayesian Information Criterion is given by (see e.g. Ref. Trotta (2008) for more details),

$\displaystyle{\rm BIC}\equiv\chi^{2}+k\ln N\,,$

(21) where $k$ is the number of parameters and $N$ the number of data points. The lower the BIC the better the model. Since $k$ and $N$ are equal within the models that we compare throughout this work, we can therefore only use the difference between the models $\chi^{2}$ as a BIC diagnostic. The latter criterion originates from an approximation of the Bayesian evidence assuming gaussianity of the posterior, a likelihood dominated regime and weak correlations between parameters. The $\chi^{2}$ goodness-of-fit values computed throughout this work are obtained constructing low-temperature MCMC as described in more details in Sec. 2.2 of Ref. Dirian et al. (2016). differences that are reported in Table 2. Degrees of significance used in this work are reported in Table  3 for definiteness of the discussion. A Bayes factor $B_{01}$ comparing model $\mathcal{M}_{0}$ against model $\mathcal{M}_{1}$ can be though of as telling betting odds of $B_{01}\,$:$\,1$ in favor of the former given the data.

For computing $B_{\nu\Lambda,\nu RR}$ we use a combination of statistical coherence and the Savage-Dickey density ratio (SDDR) (see Ref. Dirian et al. (2016) and references therein for details) that exploits the nested structure of the overall models discussed in that work. This allows one to get $B_{\nu\Lambda,\nu RR}$ in a rather economic way in writing,

[TABLE]

where $P(d|\mathcal{M}_{i})$ is the marginal likelihood (evidence) of the data $d$ given the model $\mathcal{M}_{i}$. The factor $B_{\Lambda,RR}$ appearing above is one of the main results of Ref. Dirian et al. (2016) and has been computed to be $B_{\Lambda,RR}=22.7$. The remaining factors are computed for the model $i=\Lambda,RR$ through the SDDR,

[TABLE]

which is the ratio of the marginalized one-dimensional posterior distribution to the marginalized one-dimensional prior of $\sum m_{\nu}$ obtained from the extended model, evaluated at the point where the simpler model is nested inside the extended model (i.e. at $\sum m_{\nu}=0.06$). Since we chose the same prior for both $\nu\Lambda{\rm CDM}$ and $\nu RR$ parameter spaces (in particular on $\sum m_{\nu}$) their contribution simplifies. Eq. (24) then yields,

[TABLE]

and we find, given BAPJ data,

[TABLE]

which is one of the main results of this article. It tells that $\nu RR$ is statistically equivalent to $\nu\Lambda{\rm CDM}$ with odds of $1.8\,$:$\,1$ for the latter, instead of being “moderately-to-strongly” disfavored with odds $22.7\,$:$\,1$ when the neutrino mass is fixed. The result is invariant under prior changes on $\sum m_{\nu}$ (as long as they are assumed to be equal) and leads to several implications. This result shows that allowing $\sum m_{\nu}$ to vary within $[0,5]{\rm eV}$ helps to reconcile $RR$-gravity with the data as already noticed above, but it also has the effect of penalizing the $\Lambda{\rm CDM}$ cosmology. This can be seen through the fact that applying the BIC method to compare $\nu\Lambda{\rm CDM}$ against $\nu RR$ given BAPJ leads to biased results. Indeed, comparing the shifts endured by the results of both methods in varying $\sum m_{\nu}$ we obtain,

[TABLE]

all in favor of $\Lambda{\rm CDM}$ given BAPJ data. Referring to Table 3, one can see a significant discrepancy between the results from the BIC differences and Bayes factors. While the former announces a reduction from “weak”/“moderate-to-strong” only to “weak” evidence in favor of $\nu\Lambda{\rm CDM}$, Bayesian model comparison tells that the latter “moderate-to-strong” evidence is in fact comfortably reduced to an “inconclusive” one. These two results are discrepant because of the loss of validity of the assumptions made in computing the BIC, which is only an approximation of the Bayes evidence (see footnote 4). In particular as, beside of the net maximum likelihood shift encapsulated in $\Delta\chi^{2}$ favoring the $RR$ model, Occam’s razor further penalizes the $\nu\Lambda{\rm CDM}$ one. Obviously one should therefore trust the result of Eq. (29).

In what follows, we address a rough analysis for trying to understand to which extent allowing the absolute neutrino mass to vary is beneficial for the $RR$ model, given BAPJ data. Coming back to Eq. (25), in our present context we can write,

[TABLE]

where $L_{i}$ is the value of the marginalized 1d posterior for $\sum m_{\nu}$ (normalized to its maximum) at the nesting point, $V_{i}$ is the volume of it and $P_{i}$ is the upper bound of the prior on the sum of the neutrino masses, $\sum m_{\nu}\in[0,P_{i}]$. We find,

[TABLE]

where we have set $P_{i}=5$ in either cases. One can then compute,

[TABLE]

which in particular shows that the $\Lambda$-based model provides “moderate-to-strong” evidence with odds of $62.5\,$:$\,1$ for fixing $\sum m_{\nu}=0.06{\rm~{}eV}$. Although the latter has a non-negligible contribution coming from the prior, it also has a non-negligible one from the likelihood (essentially originating from boundary effects), which is a handicap when compared against models preferring higher neutrino masses such as the nonlocal one studied in this work. As can be seen from Fig. 3, this is because the $\nu\Lambda{\rm CDM}$ marginalized posterior on $\sum m_{\nu}$ hits the lower bound of the prior at $0.06{\rm~{}eV}$, which involves a loss of posterior volume and a waste of prior one. In the $RR$ case the situation is different, since non-vanishing neutrino masses are preferred at $2\sigma$ level, exploiting therefore better the BAPJ data. This contributes to Occam’s razor effect intrinsically taken into account in Bayesian model comparison and partially explains why the $RR$ nonlocal model undergoes a favorable and significant change when compared against $\Lambda{\rm CDM}$ after allowing $\sum m_{\nu}$ to vary (Eq. (29)). Moreover, this also explains why the BIC difference effectively fails when comparing $\nu\Lambda{\rm CDM}$ against $\nu RR$ given BAPJ data, because it is only sensitive to the maximum of the posteriors, not to their entire volume.

The $RR$ nonlocal model described by the action (1) is therefore statistically equivalent (given BAPJ data) to Einstein gravity supplemented by a cosmological constant when reconsidering the prior on the neutrino sector, that is when one changes the cosmological parametrization from (11) to (13). This has been made possible exploiting an apparent degeneracy at the background level between $H_{0}$, $\sum m_{\nu}$ and the nature of the effective dark energy described by the nonlocal model, which was illustrated in Fig. 2. In what follows we provide an outlook motivating the use of additional data, in particular coming from galaxy surveys, for being able to make a distinction between the $\nu\Lambda{\rm CDM}$ and the $\nu RR$ cosmological models.

IV.4 Constraints a posteriori from Redshift-Space Distortions data

Apart from secondary CMB anisotropies such as ISW or lensing effects, the constraints on $\sum m_{\nu}$ mostly come from background-geometrical features when considering the BAPJ dataset. However, massive neutrinos give rise to characteristic inhomogeneous and anisotropic signatures induced by their thermal velocity flow. In particular, they do not cluster inside regions delimited by their free-streaming scale. Below that scale, the neutrino perturbations are smoothed out and this causes a suppression of the late-time matter power spectrum at mid-to-small cosmological scales, a decrease of the lensing power and of the growth of structures in a scale-dependent manner within the linear regime Lesgourgues and Pastor (2014); Hu et al. (1998), as well as non-linear effects Adamek et al. (2016); Inman and Pen (2017). Additional data putting stronger constraints on these features are therefore relevant to include into the global fit. Nevertheless, the presence of an appreciable fraction of massive neutrinos can have partial degenerate effects with a positive fifth force present in modified gravity scenarios. A fifth force should be present at late-times in the $RR$ model and it was found in Refs. Barreira et al. (2014a); Dirian et al. (2016) that it enhances the growth of linear and non-linear structures compared to the one described by $\Lambda{\rm CDM}$. Into the latter reference, constraints on the linear growth rate of structure modelled by $f\sigma_{8}$ were derived given RSD data. It was found that the $\Lambda{\rm CDM}$ model was favored over the nonlocal one and this was quantified by a BIC of $\Delta\chi^{2}=2.6$. This value was computed a posteriori, that is, given that $f\sigma_{8}$ was derived from each model on its BAPJ bestfit. In this part, we study the impact of a massive neutrino component on the linear growth rate of structure using the same method.

The degenerate effects present between a massive neutrino fraction and linear growth rate are well-illustrated from the degeneracy direction observed in the $\sigma_{8}\,$–$\,\sum m_{\nu}$ plane in Fig. 3, where one can see that they are anti-correlated: the higher the massive neutrino fraction $\Omega_{\nu}\sim\sum m_{\nu}$, the lower $\sigma_{8}$. Given Planck data, the mean value inferred on $\sigma_{8}$ for $\Lambda{\rm CDM}$ is smaller than the one provided by $RR$, in agreement with the higher growth within the nonlocal model, and their mean values are generically smaller in the $\nu$-extended case. For Planck only, we find that the departure of the best-fit value of $\sum m_{\nu}$ from the lower bound of the prior in $RR$ cosmology is caused by the addition of the Planck CMB lensing power spectrum which is sensitive to a weighted projection of density fluctuations along the line-of-sight. Joining BAO+SNIa data pulls the total matter density fraction $\Omega_{m}$ to higher values, involving a stronger increase in the absolute neutrino mass that preserves the value of $\sigma_{8}$ close to the one inferred in $\Lambda{\rm CDM}$. Focusing on the growth, Fig. 4 shows the time evolution of $f\sigma_{8}$ for the best-fitting $\Lambda{\rm CDM}$, $RR$ and $\nu RR$ models to the BAPJ dataset. As anticipated, the growth rate is lower in $\nu RR$ compared to the $RR$ model. The figure also displays the most recent observational determinations of $f\sigma_{8}$ from the DR12 BOSS analysis Alam et al. (2016) (grey symbols with errorbars). Using these data, the reduced $\chi_{\rm red}^{2}$ values for the $\Lambda{\rm CDM}$, $RR$ and $\nu RR$ models are, respectively555These values do not consider the mid-redshift data point. This is because the associated galaxy sample completely overlaps with those of the other two points which are independent. The number of degrees of freedom is therefore two., $\chi^{2}_{\rm red}=0.58$, $\chi^{2}_{\rm red}=1.38$ and $\chi^{2}_{\rm red}=0.97$. More pragmatically for comparison with previous results, we compute the corresponding $\chi^{2}$ values using the same data points as in Ref. Dirian et al. (2016), that is, the ones collected from 6dF GRS Beutler et al. (2012) at $f\sigma_{8}(0.067)=0.423\pm 0.055$, SDSS LRG Oka et al. (2014) at $f\sigma_{8}(0.3)=0.49\pm 0.08$ , SDSS MGS Howlett et al. (2014) at $f\sigma_{8}(0.15)=0.63^{+0.027}_{-0.24}$, BOSS LOWZ Chuang et al. (2016) at $f\sigma_{8}(0.32)=0.371\pm 0.091$, BOSS CMASS Samushia et al. (2014) at $f\sigma_{8}(0.57)=0.441\pm 0.0434$ 666Replacing the latter BOSS data by those of Ref. Alam et al. (2016) does not significantly affect our statistical conclusions., WiggleZ Blake et al. (2012) at $f\sigma_{8}(0.44)=0.413\pm 0.08$, $f\sigma_{8}(0.6)=0.39\pm 0.063$, $f\sigma_{8}(0.73)=0.437\pm 0.072$) and VIPERS de la Torre et al. (2013) at $f\sigma_{8}(0.8)=0.47\pm 0.08$. The corresponding goodness-of-fit read

[TABLE]

which shows that the fit is indeed improved in going from $RR$ to $\nu RR$ with BIC values changing from $\Delta\chi^{2}=2.6$ to $\Delta\chi^{2}=1.3$ in favor of $\Lambda{\rm CDM}$. Therefore, we can conclude that allowing $\sum m_{\nu}$ to be a free parameter helps to decrease the discrepancy of the $RR$ nonlocal gravity model with growth rate measurements and brings down the total discrepancy from $\left.\Delta\chi^{2}\right|_{\rm bapj}+\left.\Delta\chi^{2}\right|_{\rm rsd}^{\rm post}=8.5$ ($\sim 3\sigma$) to $\left.\Delta\chi^{2}\right|_{\rm bapj}+\left.\Delta\chi^{2}\right|_{\rm rsd}^{\rm post}=4.6$ ($\sim 2\sigma$) given BAPJ+(post)RSD, which induces a significant change in the (although approximated) statistical conclusion.

In turn, this shows that the data considered in this work do not possess enough constraining power to clearly distinguish between $\Lambda{\rm CDM}$ and $RR$ cosmologies. Nevertheless, the situation is expected to be different for a survey like Euclid Amendola et al. (2013). This is illustrated by the black symbols in Fig. 4, which show an estimate of the forecast errorbars for this future mission (taken from Fig. 3 of Ref. Majerotto et al. (2012)), centred around the $\Lambda{\rm CDM}$ prediction. One notes that the difference between $\Lambda{\rm CDM}$ and $\nu RR$ is larger than the forecast precision of Euclid for $z<1$, from which we can conclude that, despite partial degeneracies between the effects of massive neutrinos and the $RR$-modifications to gravity, there is still room for future RSD data to be used to help distinguishing between $\Lambda{\rm CDM}$ and $RR$ cosmologies. In addition, constraints using weak gravitational lensing data such as those of CFHTLenS Heymans et al. (2012) or KiDS Kuijken et al. (2015) also prove to be of particular interest for constraining the $\nu$-extended models. Indeed, measuring the cosmic shear induced by the large-scale structure, these data allow to put constraints on the nature of the dark energy as well as on the absolute neutrino mass Joudaki et al. (2016), although to a smaller extent than the BAPJ dataset considered in this work. However, weak lensing measurements reported in Refs. Heymans et al. (2012); Kuijken et al. (2015) are in tension with Planck CMB observations, given the $\Lambda{\rm CDM}$ model, and systematic issues first need to be addressed before these data can be used in combination with Planck for constraining modified gravity models.

As a final remark, a larger fraction of massive neutrinos in cosmological models contributes to an increased scale-dependence in the linear growth of structure. This may raise some concerns when confronting models like the best-fitting $\nu RR$ model to BAPJ against $f\sigma_{8}$ values, because the latter are usually extracted from galaxy survey data using RSD models assuming the growth to be scale-independent (see e.g. Ref. Johnson et al. (2014) for an exception to this fact and Ref. Barreira et al. (2016) for a validation study of RSD modelling in DGP gravity which exhibits scale-independent linear growth). If the scale-dependence in the $\nu RR$ is non-negligible compared to the precision targeted, extra care is required in the analysis of the data before observational constraints can be performed. To test such a fact, we plot in Fig. 4 the evolution of $f\sigma_{8}$ in the $\nu RR$ model for $k=0.01\ h/{\rm Mpc}$ (red solid) and $k=0.5\ h/{\rm Mpc}$ (red dashed). One notes that the $k$-dependence is small compared to the expected precision of Euclid, which suggests that standard methods can be used to constrain the $\nu RR$ model.

IV.5 A word on $H_{0}$

Another interesting outcome of the constraints on the $RR$ model relates to the preferred values of $H_{0}$. For the best-fitting $\Lambda{\rm CDM}$ model to the BAPJ dataset, one finds $H_{0}=67.67^{+0.47}_{-0.50}\ {\rm km/s/Mpc}$ (cf. Table 1), which lies $\sim 1\sigma$ below the determination from local measurements discussed in Ref. Efstathiou (2014), which sets $H_{0}=70.6\pm 3.3\ {\rm km/s/Mpc}$ (note that this value becomes $H_{0}=72.5\pm 2.5\ {\rm km/s/Mpc}$ if other assumptions are made into the analysis). More recently, the work of Ref. Riess et al. (2016) sets a higher value $H_{0}=73.24\pm 1.74\ {\rm km/s/Mpc}$ (see also Refs. Riess et al. (2011); Humphreys et al. (2013)). Furthermore, recent determinations of $H_{0}$ using hyperparameters, $H_{0}=73.75\pm 2.11\ {\rm km/s/Mpc}$ Cardona et al. (2017), or from gravitational lensing time delay methods, $H_{0}=71.9^{+2.4}_{-3.0}\ {\rm km/s/Mpc}$ Bonvin et al. (2017), are also significantly away from the $\Lambda{\rm CDM}$ bestfit.

The seriousness of the above-mentioned $H_{0}$ tensions is still subject to current debates and one still needs to understand better the role of systematics before claiming the need of new physics (see e.g. Refs. Di Valentino et al. (2016); Luković et al. (2016); Bernal et al. (2016)). Nevertheless, taking the current measurements at face value, one notes that for the best-fitting $\nu RR$ model to the BAPJ dataset one has $H_{0}=69.49^{+0.79}_{-0.80}\ {\rm km/s/Mpc}$, which significantly ameliorates the agreement with the local determinations and would therefore improve further the global fit.

IV.6 The importance of terrestrial determinations of $\sum m_{\nu}$

The current constraints on neutrino masses that are independent of cosmology arise from terrestrial experiments. The lower bounds on $\sum m_{\nu}$ come from neutrino oscillations experiments which, assuming a massless eigenstate, set $\sum m_{\nu}\gtrsim 0.05\ {\rm eV}$ and $\sum m_{\nu}\gtrsim 0.1\ {\rm eV}$ for normal and inverted mass hierarchies respectively. The current best upper bounds are obtained by analysing the high-energy part of the spectrum of Tritium $\beta$-decay in experiments such as MAINZ and TROITSK and set the electron neutrino mass to $m_{\nu_{e}}\lesssim 2.2\ {\rm eV}$ ($2\sigma$) which corresponds to $\sum m_{\nu}\lesssim 6.6\ {\rm eV}$ in our context. Future Tritium $\beta$-decay experiments such as KATRIN will be sensitive to mass scales $\sum m_{\nu}\lesssim 0.6\ {\rm eV}$ at $90\%$ confidence level. The sensitivity can be even better if neutrinos turn out to be Majorana particles, in which case neutrinoless double $\beta$ decay experiments should be able to probe the region corresponding to $\sum m_{\nu}\gtrsim 0.3\ {\rm eV}$ with high precision777Note that the quoted upper bounds assume specific models of nuclear matrix elements. (see e.g. Refs. Drexlin et al. (2013); Dell’Oro et al. (2016); Vergados et al. (2016); Engel and Menéndez (2017) for reviews). These forecast sensitivities can therefore be proven useful for confirming cosmological observations. As it has been shown throughout our study, the determination of the absolute neutrino mass scale from cosmological probes depends on the assumed cosmological model. As such, if terrestrial neutrino experiments will detect non-minimal neutrino masses, we will need to modify the standard $\Lambda{\rm CDM}$ cosmological model.

V Summary & Conclusion

We have revisited the observational constraints of Ref. Dirian et al. (2016) where the nonlocal model of modified gravity described in Sec. II was found to be disfavored against $\Lambda{\rm CDM}$ with odds of $22:1$, given Planck+SNIa+BAO data. Such a discrepancy was noticed to be mostly caused by a CMB-SNIa tension in the $\Omega_{m}\,$–$\,H_{0}$ plane, that we analyzed in more details into our Sec. IV. We have found that it results from the quite smooth, late-time and phantom nature of the effective dark energy described by the nonlocal model which, for fixed parameter values, induces a decrease on the late-time Hubble expansion rate $H(z\approx 0)$ compared to that described by $\Lambda{\rm CDM}$. Such a fact generically implies a smaller acoustic scale $\theta_{*}$ for the CMB which is corrected by the inference of a higher value of $H_{0}$ given Planck data, as well as a larger luminosity distance that is compensated by a larger $\Omega_{m}$ given SNIa data. Since the shape information from CMB temperature power spectrum constrains well $\omega_{m}$, which is a multiplicative combination of $H_{0}$ and $\Omega_{m}$, the trends inferred from the nonlocal model are contradictory and a tension appears, as illustrated in Fig. 2.

We have then shown that allowing the absolute neutrino mass to be a free parameter in the nonlocal gravity model resolved the tension. The $\nu$-extended nonlocal model, denoted $\nu RR$, ends up to be statistically equivalent to $\nu\Lambda{\rm CDM}$ given Planck+SNIa+BAO data with odds of $1.8:1$ in favor of the latter. We have shown that the compatibility between $\nu RR$ and $\nu\Lambda{\rm CDM}$ was caused by a better fit of the nonlocal model to the data, but also by the Occam’s razor effect penalizing the $\nu\Lambda{\rm CDM}$ model because of its preference for small absolute neutrino masses. As a result, the absolute neutrino mass is inferred to be non-zero $\sum m_{\nu}>0$ at $\sim 2\sigma$ level given the nonlocal model, with the best-fitting value $\sum m_{\nu}\approx 0.21\ {\rm eV}$. We have then placed constraints from RSD data a posteriori on both models, i.e. considering the matter power spectra corresponding to their respective Planck+SNIa+BAO bestfit. These constraints have been shown to be improved as well by the presence of a higher neutrino fraction $\Omega_{\nu}$ into the nonlocal cosmology. Further determinations from local measurements of $H_{0}$ were also discussed, as these are in better agreement with the nonlocal gravity model inferring a value of $H_{0}=69.49^{+0.79}_{-0.80}\ {\rm km/s/Mpc}$, which is $\sim 2\sigma$ above $H_{0}=67.67^{+0.47}_{-0.50}\ {\rm km/s/Mpc}$, inferred from $\Lambda{\rm CDM}$.

In conclusion, letting the absolute neutrino mass to be a free parameter allowed the resulting $\nu RR$ nonlocal gravity model to fit current Planck+SNIa+BAO data as well as $\nu\Lambda{\rm CDM}$. Given these data, the nonlocal model provides $\sum m_{\nu}>0$ at $\sim 2\sigma$ with bestfit $\sum m_{\nu}\approx 0.21\ {\rm eV}$. This is in disagreement with the value inferred from $\nu\Lambda{\rm CDM}$ which prefers $\sum m_{\nu}=0\ {\rm eV}$, corresponding to the lower extreme value of the prior chosen in this work. Furthermore, allowing $\sum m_{\nu}$ to take higher values within the nonlocal cosmology also improved the fit to RSD data done a posteriori. This provides one more example showing that the cosmological constraints on the absolute neutrino mass depends on the assumed cosmological model, because of degenerate effects between modifications to gravity and massive neutrinos. Still, our study also suggests that the use of additional data coming from future galaxy redshift surveys could reduce such a degeneracy in the studied case and potentially discriminate between the $\nu\Lambda{\rm CDM}$ and $\nu RR$ models. We have provided an illustration to this fact in considering forecast constraints from Euclid RSD data. A more quantitative analysis is left for future work.

Acknowledgements.

The author sincerely thanks Alexandre Barreira for his useful contributions to early stages of this work and for fruitful discussions. Many thanks to Martin Kunz and Michele Maggiore for their useful comments about the manuscript, to Leïla Haegel for sharing her expertise in experimental neutrino research and to the anonymous referee for many valuable comments contributing to the clarity of the manuscript. Numerical work presented in this publication used the Baobab cluster of the University of Geneva. YD is supported by the Fonds National Suisse.