Early dark energy constraints on growing neutrino quintessence cosmologies

TL;DR

This paper explores how early dark energy models with a scalar field coupled to neutrinos can be constrained using Planck data, analyzing background evolution with different potentials and functions.

Contribution

It provides an analytic and numerical analysis of growing neutrino quintessence models, constraining early dark energy contributions with observational data.

Findings

Relaxation of constraints on early dark energy in exponential potentials.

Constraints on the kinetic function magnitude at early times.

Analytic solutions for background evolution during radiation and matter domination.

Abstract

We investigate cosmological models in which dynamical dark energy consists of a scalar field whose present-day value is controlled by a coupling to the neutrino sector. The behaviour of the scalar field depends on three functions: a kinetic function, the scalar field potential, and the scalar field-neutrino coupling function. We present an analytic treatment of the background evolution during radiation- and matter-domination for exponential and inverse power law potentials, and find a relaxation of constraints compared to previous work on the amount of early dark energy in the exponential case. We then carry out a numerical analysis of the background cosmology for both types of potential and various illustrative choices of the kinetic and coupling functions. By applying bounds from Planck on the amount of early dark energy, we are able to constrain the magnitude of the kinetic function at early times.

Figures

Equations

$$\begin{gathered} S=\int {\text{d}}^{4}x\sqrt{-g}\left[\frac{1}{2}{M}_{\text{P}}^{2}R-\frac{1}{2}{k}^2(\phi ){\nabla }_{\mu }\phi {\nabla }^{\mu }\phi -V(\phi )\right] \\ +S_{\text{m}}[{\Psi }_{\text{m}},g_{\mu \nu }]+S_{\gamma }[{\Psi }_{\gamma },g_{\mu \nu }] \\ +S_{\nu }[{\Psi }_{\nu },C(\phi {)}^2{g}_{\mu \nu }], \end{gathered}$$

$$m_{\nu }(\phi )={\bar{m}}_{\nu }C(\phi ),$$

$$\beta (\phi )\equiv -M_{\text{P}}\frac{\text{d}\ln C(\phi )}{\text{d}\phi }.$$

$$G_{\mu\nu}=\frac{1}{M_{\text{P}}^{2}}T_{\mu\nu}+\frac{1}{M_{\text{P}}^{2}}\biggl{[}k^{2}(\varphi)\nabla_{\mu}\varphi\nabla_{\nu}\varphi\\ -\frac{1}{2}k^{2}(\varphi)\nabla_{\rho}\varphi\nabla^{\rho}\varphi g_{\mu\nu}-V(\varphi)g_{\mu\nu}\biggr{]}\,,$$

$$-k^2{\nabla }_{\mu }{\nabla }^{\mu }\phi -\frac{1}{2}\frac{\text{d}{k}^2}{\text{d}\phi }{\nabla }_{\mu }\phi {\nabla }^{\mu }\phi +\frac{\text{d}V}{\text{d}\phi }+\beta \frac{T_{\mu }^{(\nu )\mu }}{M_{\text{P}}}=0.$$

$$\text{d}{s}^2=-\text{d}{t}^2+a(t{)}^2{\delta }_{ij}\text{d}{x}^i\text{d}{x}^j,$$

$$H^2=\frac{\rho }{3M_{\text{P}}^2},$$

$$\dot{H}=-\frac{\rho +p}{2M_{\text{P}}^2},$$

$$\ddot{\phi }+3H\dot{\phi }+\frac{1}{2k^2}\frac{\text{d}{k}^2}{\text{d}\phi }{\dot{\phi }}^2+\frac{1}{k^2}\frac{\text{d}V}{\text{d}\phi }-\frac{\beta }{M_{\text{P}}}({\rho }_{\nu }-3p_{\nu })=0,$$

$${\rho }_{\phi }=\frac{k^2}{2}{\dot{\phi }}^2+V,$$

$$p_{\phi }=\frac{k^2}{2}{\dot{\phi }}^2-V.$$

$${\dot{\rho }}_{\nu }+3H({\rho }_{\nu }+p_{\nu })=-\frac{\beta }{M_{\text{P}}}({\rho }_{\nu }-3p_{\nu })\dot{\phi },$$

$${\dot{\rho }}_{\phi }+3H({\rho }_{\phi }+p_{\phi })=\frac{\beta }{M_{\text{P}}}({\rho }_{\nu }-3p_{\nu })\dot{\phi }.$$

$${\rho }_{\text{tot}}={\rho }_{\text{d}}+{\rho }_{\phi },$$

$${\rho }_{\text{d}}\propto \exp (-nN),$$

$${\rho }_{\phi }\propto \exp (-nN),$$

$${\Omega }_{\phi }=\frac{n{k}_{\text{c}}^2}{{\alpha }^2}.$$

$$\phi =M_{\text{P}}\frac{nN}{\alpha }+\hat{\phi },$$

$${\Omega }_{\phi }=\frac{n{k}^2}{{\alpha }^2}(1-\bar{u}),$$

$$\bar{u}=\frac{M_{\text{P}}}{\alpha (1-{\Omega }_{\phi })}\frac{\text{d}\log k^2}{\text{d}\phi }.$$

$$k_1^2(\phi )=\frac{M_{\text{P}}\alpha }{2\kappa (\phi -\bar{\phi })},$$

$$\bar{u}=-\frac{2\kappa {\Omega }_{\phi }}{n(1-{\Omega }_{\phi })}.$$

$$\bar{u}=-\frac{2\kappa }{n(1-{\Omega }_{\phi })},$$

$$V(\phi )=M_{\text{P}}^4\tilde{V}(M_{\text{P}}/\phi {)}^{\lambda },$$

$$\ddot{\phi }+3H\dot{\phi }+\frac{1}{2k^2}\frac{\text{d}{k}^2}{\text{d}\phi }{\dot{\phi }}^2+\frac{1}{k^2}\frac{\text{d}V}{\text{d}\phi }=0.$$

$$k^2(\phi )=\frac{M_{\text{P}}}{2\kappa \phi }.$$

$$\ddot{\phi }+3H\dot{\phi }-\frac{{\dot{\phi }}^2}{2\phi }-2\lambda \kappa M_{\text{P}}^{3+\lambda }\tilde{V}{\phi }^{-\lambda }=0.$$

$$H^2{\phi }^{\prime \prime }+H{H}^{\prime }{\phi }^{\prime }+3H^2{\phi }^{\prime }-\frac{H^2{\phi }^{\prime 2}}{2\phi }-2\lambda \kappa M_{\text{P}}^{3+\lambda }\tilde{V}{\phi }^{-\lambda }=0,$$

$$\phi =M_{\text{P}}\exp (\Phi /M_{\text{P}}),$$

$$\begin{gathered} {\Phi }^{\prime \prime }+\frac{{\Phi }^{\prime 2}}{2M_{\text{P}}}+\left(\frac{H^{\prime }}{H}+3\right){\Phi }^{\prime } \\ -2\lambda \kappa \frac{M_{\text{P}}^3\tilde{V}}{H^2}\exp (-(\lambda +1)\Phi /M_{\text{P}})=0. \end{gathered}$$

Full text

††preprint: APS/123-QED

I Introduction

Since the acceleration of the expansion of the Universe was discovered by the study of Type Ia supernovae Riess et al. (1998); Schmidt et al. (1998), uncovering the cause of this acceleration has become one of the most important problems in cosmology. The so-called concordance model of cosmology, or $\Lambda$CDM, is the simplest cosmological model that gives a good fit to the Type Ia supernovae data Riess et al. (2016) as well as the cosmic microwave background Aghanim et al. (2018) and large scale structure Abolfathi et al. (2018). In $\Lambda$CDM the late-time universal expansion is caused by a cosmological constant denoted by $\Lambda$.

Despite its experimental successes, however, $\Lambda$CDM suffers from a number of theoretical problems, notably the coincidence problem Velten et al. (2014), the fine-tuning problem Martin (2012) and the difficulty of explaining what fundamentally gives rise to the cosmological constant Polchinski (2006). This has motivated the study of a number of possible explanations of the observed late-time expansion not involving a cosmological constant, including modified gravity models Clifton et al. (2012) and dynamical dark energy Copeland et al. (2006). Quintessence is a simple form of dynamical dark energy which consists of a single scalar field $\varphi$ whose evolution is slow enough to give rise to a negative pressure, hence accelerated expansion and late time domination of the energy density Copeland et al. (2006); Tsujikawa (2013).

One novel approach to solving the coincidence problem is to give $\varphi$ a non-standard coupling to neutrinos such that the neutrino mass is $\varphi$-dependent. The scalar field evolves according to a scaling solution Wetterich (1988); Copeland et al. (1998) that ends as the neutrinos become non-relativistic and, through the scalar-neutrino coupling, provide a force on $\varphi$ that effectively stops it from evolving. At this point $\varphi$ is potential-dominated and acquires a negative equation of state, giving rise to a period of dark energy domination. Such models are known as “growing neutrino quintessence” models Wetterich (2007); Amendola et al. (2008), so called because the neutrino mass grows as the scalar field evolves. Growing neutrino quintessence models are closely related to “mass varying neutrino” (MaVaN) models, which also involve an interaction between dark energy and the neutrino mass Fardon et al. (2004); Afshordi et al. (2005); Kaplan et al. (2004); Bi et al. (2005); Brookfield et al. (2006); Spitzer (2006); Takahashi and Tanimoto (2006); La Vacca and Mota (2013). As well as providing the right conditions for the scalar field to play the role of dark energy, the neutrino-scalar coupling can give rise to an attractive fifth-force acting on the neutrinos that is much stronger than gravity. This force gives rise to non-linear “neutrino lumps” on large scales, which have been extensively studied using linear approximation Mota et al. (2008); Pettorino et al. (2010), N-body simulations Baldi et al. (2011); Ayaita et al. (2012, 2013, 2016); Führer and Wetterich (2015); Casas et al. (2016), spherical collapse Wintergerst and Pettorino (2010), and other methods Wintergerst et al. (2010); Nunes et al. (2011); Brouzakis et al. (2011). The effect these neutrino lumps have on the cosmological history depends on the masses of neutrinos. As found in Ref. Casas et al. (2016), for large neutrino masses the neutrino lumps can be stable and can lead to significant backreaction effects. For smaller neutrino masses, however, the neutrino lumps are unstable; they form and dissociate periodically such that backreaction effects are small.

In this work we consider a number of growing neutrino quintessence models, investigating the effects of changing the scalar field kinetic and potential terms and the scalar-neutrino coupling on the background evolution of the scalar field. We study radiation and matter domination analytically, and numerically solve the full system of background equations, in order to gain insight into the robustness of growing neutrino quintessence to variations in the model parameters.

This paper is organised as follows: in Sec. II we describe a particular growing neutrino quintessence model proposed by Wetterich in Ref. Wetterich (2015) and state the equations of motion; in Sec. III we present approximate analytic solutions to the model in Ref. Wetterich (2015) and a related model; in Sec. IV we describe our numerical analysis, present the numerical results and discuss the constraints on model parameters; and finally we summarise in Sec. V.

II The Model

In Ref. Wetterich (2015), Wetterich proposed a unified model of inflation and quintessence. The model consists of a crossover from a past, ultraviolet fixed point to a future, infrared fixed point. Each of these regions is approximately scale invariant and corresponds to inflation and dark energy domination respectively. This is achieved by introducing a scalar field, which in a particular choice of conformal frame acts as a variable Planck mass. The scalar field is coupled to the neutrinos in such a way as to produce growing neutrino quintessence Wetterich (2007); Amendola et al. (2008), with the transition to dark energy domination being driven by the “trigger” event of the neutrinos becoming non-relativistic.

In this work we consider the crossover region, since it can be constrained by presently available cosmological data. According to the model in Ref. Wetterich (2015), radiation domination, matter domination and the present transition to dark energy domination are all part of the crossover. In this period the model effectively reduces to a scalar field minimally coupled to both gravity and the matter sector, with the only exception being the coupling to the neutrino sector. For this reason it is convenient to work in the Einstein frame, in which the model is described by the following action:

[TABLE]

where $M_{\text{P}}$ is the reduced Planck mass, $k^{2}(\varphi)$, $V(\varphi)$ and $C(\varphi)$ are respectively the kinetial, potential and neutrino-scalar coupling function and must be specified in order to choose a particular model. $\Psi_{\text{m}}$, $\Psi_{\gamma}$ and $\Psi_{\nu}$ correspond to the matter, radiation and neutrino fields respectively. The key feature of growing neutrino quintessence models is the function $C(\varphi)$ which couples neutrinos to the scalar field and effectively gives the neutrinos a time-dependent mass given by:

[TABLE]

where $\bar{m}_{\nu}$ is a mass scale. For simplicity we take all neutrino masses to be equal. It is often convenient to work in terms of the dimensionless function:

[TABLE]

By varying the action, Eq. 1, with respect to the metric $g_{\mu\nu}$ one obtains the gravitational field equations:

[TABLE]

and varying with respect to the scalar field $\varphi$ yields the scalar field equation of motion:

[TABLE]

Here $T_{\mu\nu}$ is the total stress-energy-momentum tensor of all species apart from the scalar field and $T_{\mu\nu}^{(\nu)}$ is the stress-energy-momentum tensor for the neutrinos.

If we assume a spatially-flat Friedmann-Lemaître-Robertson-Walker metric of the form

[TABLE]

and assume the scalar field is homogeneous, then Eqs. 4 and 5 become:

[TABLE]

[TABLE]

and

[TABLE]

where $H\equiv\dot{a}/a$ is the Hubble parameter, dots denote differentiation with respect to time, $\rho=\rho_{\text{m}}+\rho_{\nu}+\rho_{\gamma}+\rho_{\varphi}$ and $p=p_{\text{m}}+p_{\nu}+p_{\gamma}+p_{\varphi}$ are the energy density and pressure of all species. The energy density and pressure of the homogeneous scalar field are defined as:

[TABLE]

[TABLE]

Matter and radiation obey the usual conservation equations: $\dot{\rho}_{\text{m}}+3H\rho_{\text{m}}=0$ and $\dot{\rho}_{\gamma}+4H\rho_{\gamma}=0$. However, the neutrinos obey a modified conservation equation due to their interaction with the scalar field:

[TABLE]

[TABLE]

For most of the Universe’s history, the neutrinos are highly relativistic and $\rho_{\nu}-3p_{\nu}\approx 0$ such that the scalar field and the neutrinos are effectively uncoupled and the scalar field energy density tracks that of the dominant species. After the neutrinos become non-relativistic the coupling becomes important and, for large enough $|\beta|$, effectively stops the evolution of the scalar field by providing a force to counter that caused by the gradient of the potential in Eq. 9. As a result, the scalar field’s energy density and pressure are dominated by the potential and the equation of state $w_{\varphi}\equiv p_{\varphi}/\rho_{\varphi}$ approaches $-1$, which is consistent with observations Aghanim et al. (2018).

III Approximate analytic solutions

Under certain simplifying assumptions, it is possible to solve the scalar field equation, Eq. 9, analytically. In this section we consider the behaviour of $\varphi$ before the neutrinos become non-relativistic, both for an exponential and an inverse power law potential. For the exponential case the scalar field evolves linearly with $N=\log(a)$ and there is an approximately constant fraction of early dark energy present. In the inverse power law case we find instead that $\log(\varphi)$ evolves linearly with $N$. For an analytic treatment of a related but distinct model, see Ref. Hossain et al. (2014).

III.1 Exponential potential

First we consider a constant kinetic function $k^{2}(\varphi)=k^{2}_{\text{c}}=\text{const}$ and an exponential potential $V(\varphi)=M_{\text{P}}^{4}\exp(-\alpha\varphi/M_{\text{P}})$, where $\alpha$ is a dimensionless parameter that determines the slope of the potential. Before the neutrinos have become non-relativistic, the model exhibits a scaling solution whereby the energy density of the scalar field tracks that of the dominant species (radiation or matter, depending on the epoch) with the result that the energy density fraction of the scalar field is constant. It is convenient to introduce the energy density of the dominant species as $\rho_{\text{d}}$, which is equal to $\rho_{\gamma}+\rho_{\nu}$ in the radiation-dominated epoch and $\rho_{\text{m}}$ in the matter-dominated epoch. Here we are considering only the epoch in which neutrinos are highly relativistic, so the coupling between the neutrinos and the scalar field is effectively zero and neutrinos can be treated simply as radiation along with the photons.

Sufficiently far from matter-radiation equality one can neglect whichever of matter and radiation is subdominant and write:

[TABLE]

where the energy density of the dominant species evolves as

[TABLE]

where $N=0$ at the present time and $n=4\,(3)$ for radiation (matter) domination. In the scaling solution,

[TABLE]

and the (constant) energy density fraction of the scalar field is given by

[TABLE]

The scalar field itself obeys the following particular solution of Eq. 9

[TABLE]

where $\hat{\varphi}$ is the value $\varphi$ would take at $N=0$ (though note that this bears no relation to realistic present-day values of $\varphi$ since at some point before $N=0$ the neutrinos become important and the scaling solution becomes invalid).

For a slowly-varying kinetic function $k^{2}(\varphi)$ we can expect behaviour that approximates this scaling solution. The procedure for quantifying the deviation from the scaling solution was demonstrated in Ref. Wetterich (2015), which in turn is based on a calculation in Ref. Wetterich (1995). When we carry out this procedure we find disagreement with the results of Ref. Wetterich (2015). The details of our calculation are laid out in Appendix A. Here we simply present the results.

We find that the energy density of the scalar field obeys

[TABLE]

which deviates from the exact scaling solution by a small quantity

[TABLE]

This differs from the corresponding result in Ref. Wetterich (2015) by a factor of $\Omega_{\varphi}$. As an example, we consider the particular kinetic function used in Ref. Wetterich (2015):

[TABLE]

where $\kappa$ is a dimensionless parameter which sets the scale of the function $k_{1}^{2}(\varphi)$. Substituting into Eq. 20, one obtains

[TABLE]

The corresponding result in Ref. Wetterich (2015) is given by

[TABLE]

from which it follows that $\kappa$ must be small compared to $1$, in order to give a small $\bar{u}$ and hence produce behaviour close to the scaling solution. However, since $\Omega_{\varphi}$ is small, we find no such constraint on $\kappa$; $\bar{u}$ is small automatically in Eq. 22.

This has implications for the prospects of constraining the model. A larger value of $\kappa$ gives smaller values of the function $k_{1}^{2}(\varphi)$ and hence smaller values of $\Omega_{\varphi}$ Wetterich (2015). There is a tight upper bound from the Planck experiment on the value of $\Omega_{\varphi}$ at early times. This can translate into a lower bound on $\kappa$, which will be discussed in more detail in Sec. IV. Based on Eq. 23 one would conclude that there are both upper and lower bounds on $\kappa$, which could potentially put a very tight constraint on the model. However, based on our result for $\bar{u}$, which is small irrespective of the magnitude of $\kappa$, one finds no upper bound on $\kappa$. As will be shown in Sec. IV, we can consider values of $\kappa$ much larger than the upper bound found in Ref. Wetterich (2015). Our numerical results match closely our prediction and there is no evidence of any approximation breaking down for large $\kappa$.

III.2 Inverse power law potential

An approximate analytic solution can also be found for models with inverse power law potentials of the form

[TABLE]

where $\tilde{V}$ and $\lambda$ are dimensionless constants.

While the neutrinos are relativistic, Eq. 9 becomes

[TABLE]

Using the same kinetic function as for the exponential potential case, Eq. 21, but with $\alpha=1$ and $\bar{\varphi}=0$ since these parameters relate to the specific model presented in Ref. Wetterich (2015), we have

[TABLE]

For our choice of $k^{2}(\varphi)$ and $V(\varphi)$, Eq. 25 becomes:

[TABLE]

In terms of e-foldings $N$ as the time variable, we have:

[TABLE]

where $\varphi^{\prime}\equiv\text{d}\varphi/\text{d}N$. Finally introducing $\Phi$ via

[TABLE]

Eq. 28 becomes:

[TABLE]

In a radiation (matter) dominated universe, the Hubble parameter evolves according to $H^{2}=\tilde{H}^{2}\exp(-nN)$ where $n=4\,(3)$ and $\tilde{H}$ is a normalising factor. Equation 30 then becomes:

[TABLE]

Motivated by results from numerical simulation (see Sec. IV), which show linear solutions for $\Phi$, we make the following ansatz:

[TABLE]

where $q$ is a dimensionless constant and $\hat{\Phi}$ is the value $\Phi$ would take if this solution were extrapolated to $N=0$.

Under this ansatz Eq. 31 becomes:

[TABLE]

Treating the $N$-dependent and $N$-independent parts of the equation separately, we obtain

[TABLE]

which, substituting into Eq. 33, gives

[TABLE]

Rearranging, we find $\hat{\Phi}$ as

[TABLE]

Thus, in contrast to the previous section, we find that inverse power law potentials admit solutions in which $\log(\varphi)$ evolves linearly with $N$ as opposed to $\varphi$ evolving linearly as in the exponential potential case.

It is also instructive to find an expression for the dark energy density fraction. Substituting our solution for $\varphi$ (Eqs. 29 and 32) into Eq. 10 gives the energy fraction:

[TABLE]

Recalling Eq. 34, we can write

[TABLE]

Thus it turns out that $\Omega_{\varphi}$ is proportional to $\varphi$:

[TABLE]

where $q$ and $\hat{\Phi}$ are given by Eqs. 34 and 36. In contrast to the exponential case, where there is an approximately constant fraction of early dark energy, here the fact the dark energy fraction has an exponential dependence on $N$ implies that at early times (i.e. large negative values of $N$), it automatically makes a negligible contribution to the energy density. These results are confirmed in Sec. IV, with Figs. 5 and 6 showing $\log(\varphi)$ and $\log(\Omega_{\varphi})$ evolving linearly with $N$ with a gradient given by $q$.

IV Numerical evolution

In addition to the analytic approach laid out in Sec. III, we numerically solved the equations of motion. This allows us to confirm the results of Sec. III and to probe the late-universe cosmology that our analytic approach did not capture.

To generate our results we modified the code used by Barreira et al. in Ref. Barreira et al. (2015), which is based on the Code for Anisotropies in the Microwave Background (CAMB) Lewis et al. (2000). We modified the background part of Barreira et al.’s code such that it solved the background equations of motion laid out in Sec. II.

We consider the following choices for the kinetic, potential and coupling functions:

Kinetic function:

• •

$k_{\text{c}}^{2}(\varphi)=\text{const}$ ,

• •

$k_{1}^{2}(\varphi)=\frac{M_{\text{P}}\alpha}{2\kappa(\varphi-\bar{\varphi})}$ .

Coupling function:

• •

$\beta_{\text{c}}(\varphi)=\text{const}$ ,

• •

$\beta_{1}(\varphi)=-\frac{M_{\text{P}}}{\varphi_{\text{c}}-\varphi}$ ,

• •

$\beta_{2}(\varphi)=-\left(\frac{M_{\text{P}}}{\varphi_{\text{c}}-\varphi}\right)^{2}$ ,

• •

$\beta_{3}(\varphi)=-\frac{\gamma M_{\text{P}}}{\varphi}$ .

Potential function:

• •

$V_{\text{exp}}(\varphi)=M_{\text{P}}^{4}\exp(-\alpha\varphi/M_{\text{P}})$ ,

• •

$V_{\text{IPL}}(\varphi)=\tilde{V}M_{\text{P}}^{4}(M_{\text{P}}/\varphi)^{\lambda}$ .

The motivation for choosing these functions is as follows. In Ref. Wetterich (2015), the functions $k^{2}_{1}(\varphi)$ and $V_{\text{exp}}(\varphi)$ are used, with $\beta(\varphi)$ unspecified. We use this as a starting point, and we specify $\beta(\varphi)=\beta_{1}(\varphi)$ as employed in Ref. Casas et al. (2016). We then widen the scope by choosing other functions that could be expected to give rise to growing neutrino quintessence behaviour. Inverse power law potentials have a qualitatively similar ‘decaying’ form to exponential potentials. The couplings $\beta_{\text{c}}$, $\beta_{1}$, $\beta_{2}$, and $\beta_{3}$ each correspond to a function $C(\varphi)$ via Eq. 3, or equivalently:

[TABLE]

The four functions $\beta(\varphi)$ considered here all correspond to a rapidly rising $C(\varphi)$. Thus $V(\varphi)$ and $C(\varphi)$ give rise to an effective potential for the scalar field that has a minimum, which is a necessary condition for growing neutrino quintessence.

In Sec. IV.1 we focus on $k^{2}_{1}(\varphi)$, $\beta_{1}(\varphi)$ and $V_{\text{exp}}(\varphi)$. The scaling solution discussed in Sec. III is verified and a constraint is found on the parameter $\kappa$ in $k_{1}^{2}(\varphi)$ due to its effect on the amount of early dark energy. In Sec. IV.2 we consider $k^{2}_{1}(\varphi)$, $\beta_{1}(\varphi)$ and $V_{\text{IPL}}(\varphi)$, which give rise to qualitatively similar behaviour for the scalar field $\varphi$ but do not produce early dark energy. We discuss the various options for $\beta(\varphi)$ in Sec. IV.3.

IV.1 Exponential potential

In this section we present the results of numerical calculations using $V_{\text{exp}}(\varphi)$, $k^{2}_{1}(\varphi)$ and $\beta_{1}(\varphi)$. During radiation and matter domination we find $\varphi$ evolving linearly with $N=\log(a)$ according to the scaling solution discussed in Sec. III. After the neutrinos become non-relativistic, $\varphi$ starts to oscillate around the minimum of the effective potential formed by $V(\varphi)$ and $\beta(\varphi)$ and comes to a halt to behave as an effective cosmological constant. This behaviour is illustrated in Fig. 1.

Figure 2 shows the evolution of the equation of state of the scalar field. It can be seen that it mimics radiation with a value of $w_{\varphi}=1/3$ when the Universe is radiation dominated, then approaches $w_{\varphi}=0$, mimicking matter when the Universe is matter dominated, and finally tends towards $w_{\varphi}=-1$ after the neutrinos halt the evolution of the scalar field and it mimics a cosmological constant. The first two regimes illustrate the scaling solution, where the energy density of the scalar field tracks that of the dominant species as discussed in Sec. III. This is also illustrated in Fig. 3, in which we have plotted the predictions of the energy density fraction of the scalar field assuming the scaling solution is exactly satisfied both for radiation and matter domination. It can be seen that in the early Universe the numerical result closely follows $\Omega_{\varphi}=4k^{2}(\varphi)/\alpha^{2}$ and at later times it follows $\Omega_{\varphi}=3k^{2}(\varphi)/\alpha^{2}$, with a transition in between, as expected.

Figure 4 shows the effect of varying the model parameter $\kappa$ in $k_{1}^{2}(\varphi)$ on the scalar field evolution and the energy density fraction of the scalar field respectively. Note that the larger the value of $\kappa$ the smaller the amount of early dark energy. This agrees with the scaling solution result, Eq. 17 in Sec. III, since $\kappa$ is effectively a constant that controls the size of the kinetic function $k^{2}_{1}(\varphi)$ as can be seen in Eq. 21.

We find that our numerical results for the evolution of dark energy are well approximated by the early dark energy parametrisation of Doran & Robbers Doran and Robbers (2006), in which the dark energy density fraction is parametrised as follows:

[TABLE]

where $\Omega_{\text{e}}$ (the fraction of early dark energy) and $w_{0}$ (the present day equation of state) are parameters to be fitted and $\Omega_{\text{DE}}^{0}$ and $\Omega_{\text{m}}^{0}$ are the present-day dark energy and matter fractions. For a given value of $\kappa$ we carry out a least-squares fitting of our numerical results to the Doran & Robbers parametrisation to give $w_{0}$ and $\Omega_{\text{e}}$. The Planck Collaboration Ade et al. (2016) finds an upper bound on the parameter $\Omega_{\text{e}}$ of $0.0036$. We find that the value of $\kappa$ required to give rise to this value of $\Omega_{\text{e}}$ is $\kappa=1.8$, with larger values of $\kappa$ resulting in smaller values of $\Omega_{\text{e}}$ and vice versa. We therefore find a lower bound on $\kappa$ of $1.8$.

As discussed in Sec. III, Ref. Wetterich (2015) finds a requirement that $\kappa\ll 1$ in order to ensure that $u$, the deviation of $\Omega_{\varphi}$ from the scaling solution at early times, is small. If this requirement were valid then the model of Ref. Wetterich (2015) would have been ruled out by the constraints on early dark energy. However, due to our finding in Sec. III that $u$ is given by Eq. 22 and not Eq. 23, we find that there is no requirement for $\kappa$ to be small and hence our constraint that $\kappa>1.8$ does not rule out the model.

Varying the parameter $\alpha$ in the potential merely results in a rescaling of $\varphi$ and does not have any effect on the physics. We also studied the case of a constant kinetic function $k^{2}_{\text{c}}$ and found that it made almost no difference to the results. This is as expected, because the varying kinetic function to which we compare it varies very slowly over the relevant part of the Universe’s history.

IV.2 Inverse power law potential

In this section we present the results for models with $V_{\text{IPL}}(\varphi)$, $k^{2}_{1}(\varphi)$ and $\beta_{1}(\varphi)$, with $\kappa=1.8$, $\alpha=1$ and $\bar{\varphi}=0$. We considered several different values of the power $\lambda$ as shown in Figs. 5 and 6. For each value of $\lambda$, an appropriate value of $\tilde{V}$ was chosen to produce the correct dark energy density fraction at the present day. For ease of comparison, the same present-day value of $\varphi$ was chosen for each value of $\lambda$, with $\varphi_{\text{c}}$ being tuned in each case to achieve this.

The choice of $\kappa=1.8$ was made for ease of comparison with the exponential potential, but has no special significance in the inverse power law case. Larger values of $\kappa$ result in an upward shift in $\varphi$ and a corresponding downward shift in $\Omega_{\varphi}$.

Compared to the models with exponential potentials already discussed, the behaviour of models with inverse power law potentials is not drastically different. During radiation and matter domination we find that $\varphi$ evolves exponentially with $N$ as opposed to linearly as it does for models with $V_{\text{exp}}(\varphi)$. However, the qualitative behaviour of the field increasing as long as neutrinos are relativistic and then effectively stopping once they become non-relativistic is still present. Figure 5 shows the evolution of the logarithm of the scalar field against $N$ for different inverse power law potentials. Before the neutrinos become non-relativistic, $\log(\varphi)$ evolves approximately linearly with a gradient of $n/(\lambda+1)$ and an intercept of $\hat{\Phi}$ as predicted in Eqs. 34 and 36.

The evolution of the energy density of the scalar field is shown in Fig. 6. From this it is clear that these models do not give rise to early dark energy; looking back in time, the energy density of the scalar field continues to drop off rapidly. The constraint on $\kappa$ that we found for exponential potentials therefore does not apply to models with inverse power law potentials.

IV.3 Coupling function

In addition to the coupling $\beta_{1}(\varphi)$ already considered, we investigated $\beta_{\text{c}}$, $\beta_{2}(\varphi)$ and $\beta_{3}(\varphi)$. None of these choices led to behaviour significantly different to the $\beta_{1}(\varphi)$ case, provided $|\beta|$ is sufficiently large at the time at which neutrinos become non-relativistic. This requirement is automatically satisfied for $\beta_{1}(\varphi)$ and $\beta_{2}(\varphi)$, since as $\varphi$ approaches $\varphi_{\text{c}}$, $|\beta(\varphi)|$ tends to infinity. The scalar field is never allowed to reach $\varphi_{\text{c}}$, however, because the neutrino coupling term in the scalar field equation, Eq. 9, always acts to decrease the value of $\varphi$. It can be seen that the value of $\varphi_{\text{c}}$ in $\beta_{1}(\varphi)$ and $\beta_{2}(\varphi)$ determines the present-day value of $\varphi$, since the latter will approach ever closer to it but can never exceed it. This is confirmed by our numerical analysis.

For $\beta_{\text{c}}$ and $\beta_{3}(\varphi)$ one does not automatically obtain large $|\beta|$ but it must be set by an appropriate choice of parameters. In the latter case this means choosing a large value of $\gamma$. The requirement on the size of $|\beta|$ is illustrated by the following result from Ref. Wetterich (2014):

[TABLE]

If $|\beta|$ is too small, the coupling term in Eq. 9 will not be large enough to counteract the potential term and the value of $\varphi$ will continue to increase. This will result in both a larger $\Omega_{\nu}$ and a smaller $\Omega_{\varphi}$. Our numerical results are consistent with Eq. 42.

V Conclusions

We have considered various growing neutrino quintessence scenarios inspired by the model proposed by Wetterich in Ref. Wetterich (2015). We studied the early-universe background evolution analytically for both exponential potentials and inverse power law potentials. In the former case we followed the procedure in Ref. Wetterich (2015), finding some disagreement with their results. In the latter case we found an analytic result for the behaviour of such models that we were able to check numerically.

Using a modified version of CAMB, we found the numerical solution to the background equations for several different kinetic, potential and neutrino-scalar coupling functions. We verified our analytic predictions, investigated the circumstances under which growing neutrino quintessence behaviour is obtained, and used early dark energy constraints Ade et al. (2016) to constrain the parameter $\kappa$ that controls the scale of the kinetic function $k^{2}_{1}(\varphi)$.

The following conditions must be met to give rise to growing neutrino quintessence:

• •

$V(\varphi)$ must have a negative gradient in order to cause the value of the scalar field to increase with time. This gradient must be sufficiently steep that $\varphi$ reaches large enough values in the late Universe to act as dark energy. Note that growing neutrino quintessence models such as the ones considered here do not require that $V(\varphi)$ be flat in the late Universe, as other quintessence models often require. The slow evolution of $\varphi$ necessary for it to mimic a cosmological constant is achieved by the presence of the neutrino coupling term, not by slow roll.

• •

$|\beta(\varphi)|$ must be sufficiently large when the neutrinos become non-relativistic that $\beta(\rho_{\nu}-3p_{\nu})$ is able to act as a strong enough restoring force to stop the evolution of $\varphi$ in Eq. 9.

For exponential potentials with slowly varying kinetic functions we found the predicted scaling solution behaviour with an approximately constant early dark energy fraction. Using existing constraints on early dark energy Ade et al. (2016) we were able to constrain the model parameter $\kappa$, which sets the scale of the kinetic function, to be larger than $1.8$, forcing it into a region previously thought excluded Wetterich (2015). However, we also found that there is no upper bound on $\kappa$ and so the model is not ruled out.

As well as the exponential potential considered in Ref. Wetterich (2015), we have also considered inverse power law potentials, since these have a qualitatively similar form to exponential potentials and so could provide the necessary conditions for growing neutrino quintessence. We confirmed that such models can give rise to growing neutrino quintessence and we found that, unlike in the case of exponential potentials, there is no early dark energy present. Planck bounds on early dark energy therefore do not translate into constraints on models with inverse power law potentials.

VI Acknowledgements

FNC is supported by a United Kingdom Science and Technology Facilities Council (STFC) studentship. AA, EJC and AMG acknowledge support from STFC grant ST/P000703/1. BL is supported by the European Research Council via grant ERC-StG-716532-PUNCA, and by STFC Consolidated Grants ST/P000541/1, ST/L00075X/1. We are grateful to Marco Baldi and Christoph Wetterich for useful discussions and to the referee for helpful and detailed comments on an earlier version of this manuscript.

Appendix A

In this Appendix we present our derivation of the results presented in Sec. III. The derivation follows that in Ref. Wetterich (2015), but we find a few disagreements with their results.

In the scaling solution, the energy density of the scalar field is a constant fraction of that of the dominant species:

[TABLE]

To find solutions close to the scaling solution we allow $f$ to vary as a function of $\varphi$:

[TABLE]

and allow a small deviation $\delta(N)$ from the scaling solution result for $\varphi$ (Eq. 18):

[TABLE]

Differentiating Eq. 44, we find:

[TABLE]

where primes denote differentiation with respect to $N$. It will be necessary to employ the $\rho_{\varphi}$ conservation equation, Eq. 13 (with $p_{\nu}=\rho_{\nu}/3$), as well as the definitions of $\rho_{\varphi}$ and $p_{\varphi}$, Eqs. 10 and 11. Using $N$ as the time variable, these are given by:

[TABLE]

[TABLE]

and

[TABLE]

It proves convenient to introduce the constant of proportionality in Eq. 15 as follows:

[TABLE]

Substituting Eqs. 47, 48, 49 and 50 into Eq. 46 yields

[TABLE]

Now, using Eqs. 44, 45 and 50

[TABLE]

Hence Eq. 51 becomes:

[TABLE]

Differentiating Eq. 45 gives

[TABLE]

Rearranging Eq. 48, we write $\varphi^{\prime}$ in terms of $\rho_{\varphi}$ and $V$:

[TABLE]

and hence

[TABLE]

where we have used Eq. 52 again. Substituting into Eq. 54 yields

[TABLE]

For constant $k^{2}$ the scaling solution is recovered, with $f=\text{const}$ and $\delta=0$. Equation 53 then gives

[TABLE]

If, however, $k^{2}$ varies smoothly we expect only a small deviation from this solution. We introduce a function $\zeta(N)$ to quantify the deviation of $f$ from the scaling solution value given by Eq. 58:

[TABLE]

Differentiating Eq. 59 gives

[TABLE]

which, using Eq. 53, gives

[TABLE]

Equations 61 and 57, both contain the term $1/(f\bar{\rho})\exp(-\alpha\delta)$, which using Eq. 59 can be written as

[TABLE]

Equations 61 and 57 can now be written as

[TABLE]

and

[TABLE]

respectively.

Now we recall Eq. 17, but introduce a small deviation $u(N)$, by

[TABLE]

and group the small functions $\delta$ and $\zeta$ through the small function $\Delta$, defined by

[TABLE]

Differentiating Eq. 66, we find

[TABLE]

We can make use of Eqs. 65 and 66 to simplify our equations for $\delta^{\prime}$ and $\zeta^{\prime}$, Eqs. 63 and 64 as follows:

[TABLE]

[TABLE]

Substituting Eqs. 68, 69 and 66 into Eq. 67 gives

[TABLE]

The differential equation for $u$ follows from differentiating Eq. 65 and rearranging as

[TABLE]

which in turn yields

[TABLE]

where we have made use of Eqs. 54 and 64 and the fact that $\Omega_{\varphi}=f/(1+f)$. Equations 72 and 70 can be compared to Eq. (108) in Ref. Wetterich (2015). We find two instances of the factor $(1-u)$ instead of $(1+u)$, and the second term in Eq. 72 differs by a factor of $\Omega_{\varphi}$. This latter difference follows through to give an extra factor of $\Omega_{\varphi}$ in Eq. 22 compared to Ref. Wetterich (2015) which, as discussed in Sec. III, has a crucial impact on the range of possible values for the parameter $\kappa$.

We continue following the procedure of Ref. Wetterich (2015) but with our versions of the $\Delta$ and $u$ equations in order to find an approximate form for $u$. Using Eqs. 68 and 69, Eq. 72 can be rewritten as

[TABLE]

Close to the scaling solution $\Delta$, $u$ and $M_{\text{P}}\text{d}\log k^{2}/\text{d}\varphi$ are all small. Expanding Eqs. 70 and 73 to linear order in small quantities gives

[TABLE]

Setting $\Delta^{\prime}=u^{\prime}=0$, we see that this system of equations admits a constant solution:

[TABLE]

One can then split $u=\bar{u}+\hat{u}$ and $\Delta=\bar{\Delta}+\hat{\Delta}$ into their $N$-independent and $N$-dependent components. The equations of motion for only the $N$-dependent parts are as follows:

[TABLE]

which can be written in the following form:

[TABLE]

where

[TABLE]

The real parts of the eigenvalues of the matrix $A$ are both negative, which implies that the $N$-dependent parts of $\Delta$ and $u$ decay with $N$. Thus the solution with $u=\bar{u}$ and $\Delta=\bar{\Delta}$ is approached. Hence, it is appropriate to use $\bar{u}$ in Eq. 65 which gives rise to Eqs. 19 and 20.