Relaxing cosmological constraints on current neutrino masses

Abstract

We show that a mass-varying neutrino model driven by scalar field dark energy relaxes the existing upper bound on the current neutrino mass to ${\sum m_{\nu}<0.72}$ eV. We extend the standard $\Lambda$ cold dark matter model by introducing two parameters: the rate of change of the scalar field with the number of $e$-folds and the coupling between neutrinos and the field. We investigate how they affect the matter power spectrum, the cosmic microwave background anisotropies and its lensing potential. The model is tested against Planck observations of temperature, polarization, and lensing, combined with baryon acoustic oscillation measurements that constrain the background evolution. The results indicate that small couplings favor a cosmological constant, while larger couplings favor a dynamical dark energy, weakening the upper bound on current neutrino masses.

I Introduction

The standard hot big bang model predicts that the Universe is filled with a background of thermal relic neutrinos, called the cosmic neutrino background, with a temperature and density of the order of the cosmic microwave background (CMB) photons Lesgourgues and Pastor (2006); Gerbino and Lattanzi (2018). Neutrinos are held in thermal equilibrium with the primordial plasma by electroweak interactions until the temperature of the Universe drops to ${T\simeq 1\,\mathrm{MeV}}$. Below this temperature, they decouple from the thermal bath and flow freely along geodesics of spacetime. Since the neutrinos are still ultrarelativistic when they decouple, they retain a relativistic Fermi-Dirac distribution even though they are no longer in thermal equilibrium. Not being subjected to the Boltzmann exponential suppression, we have far more neutrinos than would otherwise be expected. Although relic neutrinos are very abundant, there is still no direct evidence for their background because it is hard to detect at low energy level, given that they have a very small cross section with matter. There is only indirect evidence for the cosmic neutrino background, mainly through gravitational interactions, for which theoretical predictions are in excellent agreement with observations of the CMB and large-scale structure.

The existence of a neutrino mass has been demonstrated on Earth by neutrino flavor oscillation experiments, which measure the difference between the squares of the mass eigenstates Fukuda et al. (1998); Abe et al. (2014). The experiments can constrain the minimum neutrino mass in two different scenarios (normal or inverted hierarchy) for the ordering of the individual masses. The lowest limit is given by the normal ordering of the mass eigenstates, where the total sum, ${\sum m_{\nu}>0.06\,\mathrm{eV}}$, is provided by one massive neutrino and two massless neutrinos Esteban et al. (2019); Roy Choudhury and Hannestad (2020); Gariazzo et al. (2022). This in turn implies a lower limit on the current neutrino energy density in standard cosmology, using ${\Omega_{\nu}h^{2}\simeq\sum m_{\nu}/94\,\mathrm{eV}}$, where $\Omega_{\nu}$ is the neutrino density parameter and $h$ is the current value of the Hubble parameter in units of ${100\,\mathrm{km\,s^{-1}\,Mpc^{-1}}}$. We also know from the Katrin Laboratory’s $\beta$-decay experiment that the mass eigenvalues are below the eV scale Aker et al. (2019). But the strongest constraints come from cosmological observations (see, for example, Refs. Vagnozzi et al. (2017); Tanseri et al. (2022); Vagnozzi et al. (2018); Roy Choudhury and Choubey (2018); Roy Choudhury and Naskar (2019)), which provide an upper limit on the current neutrino density. This can be translated into a maximum value for the neutrino mass. Gerstein and Zeldovich originally derived this limit in the 1960s Gershtein and Zel’dovich (1966). They showed that the requirement that neutrinos do not overclose the Universe, ${\Omega_{\nu}h^{2}<1}$, suggests a cosmological constraint below a hundred eV. This limit was stronger than those of laboratory experiments at the time. Today, it is the Planck 2018 measurements of the temperature and polarization of the CMB Aghanim et al. (2020a) that provide the most robust constraints: ${\sum m_{\nu}<0.26\,\mathrm{eV}}$ at the $95\%$ confidence level (C.L.), for the standard seven-parameter model ${\Lambda\mathrm{CDM}+\sum m_{\nu}}$ (see neutrinos in cosmology Workman et al. (2022)).

Alternative cosmological scenarios that extend the parameter space, including modified gravity theories, may relax the neutrino mass bounds Alam et al. (2021); Valentino et al. (2020); Bellomo et al. (2017); Atayde and Frusciante (2023). For example, fitting a 12-parameter model to Planck and baryon acoustic oscillation (BAO) data increases the limit from ${\sum m_{\nu}<0.13\,\mathrm{eV}}$ to ${\sum m_{\nu}<0.52\,\mathrm{eV}}$ (95% C.L.) Valentino et al. (2020). This extensive scenario includes five additional free parameters: two for the parametrization of the dark energy equation of state, one for the running of the spectral index, one for the effective number of relativistic particles, and one unphysical for the amplitude of the dark matter lensing contribution to the CMB power spectra. In the present work, we also investigate whether the existing cosmological upper limit on the neutrino mass can be relaxed. Our proposed model requires only two extra free parameters describing a possible interaction between the neutrino fluid and the dark energy component Riess et al. (1998); Perlmutter et al. (1999) given by a scalar field.

We consider a mass-varying neutrino (MaVaN) scenario, where the coupling leads to an effective neutrino mass that depends on the value of the field Gu et al. (2003); Fardon et al. (2004); Peccei (2005); Brookfield et al. (2007); Wetterich (2007); Amendola et al. (2008); Ichiki and Keum (2008); Franca et al. (2009); Geng et al. (2016). We employ a minimal parametrization where the scalar field depends linearly on the number of $e$-folds Nunes and Lidsey (2004). It limits the number of additional parameters with respect to the concordance model and alleviates the initial conditions problem Peebles and Ratra (1988) thanks to the scaling behavior of the field at early times. Such parametrization has been used in the context of testing a coupling between quintessence and the electromagnetic sector, and scalar-field-dark-matter interactions da Fonseca et al. (2022a); Barros and da Fonseca (2023); da Fonseca et al. (2022b), but never utilized in the context of neutrino interactions. By testing the model with a particular dataset that combines observations of the CMB, structure growth, and background expansion, we show that the constraint on today’s mass is weakened by neutrinos of growing mass Mota et al. (2008); Nunes et al. (2011) that receive energy from the quintessence component over cosmic time.

A mechanism that couples the scalar field, as early dark energy, to the neutrinos has been proposed Sakstein and Trodden (2020); Carrillo González et al. (2023) to alleviate the Hubble tension, i.e., the discrepancy between the $H_{0}$ determinations of high and low redshift probes Aghanim et al. (2020a); Riess et al. (2019); Wong et al. (2020); Riess et al. (2021); Banerjee et al. (2021); Lee et al. (2022), but it remains to be tested with cosmological observations. Nevertheless, the Hubble tension is not the subject of this paper, since the early dark energy component in our model is insufficient to affect it.

In the next section, we present the MaVaN theory we have chosen to study, together with our specific scalar field parametrization. The phenomenology of the model is analyzed at the background level. In Sec. III, we assess the impact of the interaction on the linear perturbations, as well as the sensitivity of the observables to the coupling. We numerically compute the power spectra of matter, the CMB temperature anisotropies, and the lensing potential with the Einstein-Boltzmann code CLASS Lesgourgues (2011); Blas et al. (2011), which we have modified to compute the theoretical observables of the varying-mass neutrino model. The model is tested against observations in Sec. IV. We estimate the cosmological parameters by performing Bayesian inferences on a dataset that combines Planck measurements of the CMB with the detection of baryon acoustic oscillations. We discuss the results in Sec. V, in particular the constraints we obtain on the current neutrino mass sum $\sum m_{\nu}$ with our chosen dataset in the context of the MaVaN scenario.

II Coupling dark energy to neutrinos

Let us consider a flat Universe with vanishing curvature, spatially homogeneous and isotropic, whose expansion is parametrized by the scale factor $a$ associated with the Friedmann-Lemaître-Roberson-Walker spacetime metric. Further considering that the expansion is sourced by photons ($\gamma$), baryons (b), cold dark matter (c), neutrinos ($\nu$), and a scalar field dark energy ($\phi$) responsible for the current acceleration, the Friedmann equation reads

$$\frac{H^{2}}{H_{0}^{2}}=\Omega_{\gamma}\,a^{-4}+\left(\Omega_{b}+\Omega_{c}\right)a^{-3}+\frac{\rho_{\nu}(a)+\rho_{\phi}(a)}{\rho_{0}}\,,$$ (1)

where $H\equiv\dot{a}/a$ is the Hubble parameter, the dot denoting derivation with respect to cosmic time $t$. $H_{0}$ is the current expansion rate, ${\Omega_{i}=\kappa^{2}\rho_{i,0}/3H_{0}^{2}}$ are the present-day density parameters, $\rho_{i}$ being the energy density of species $i$ and ${\rho_{0}={3H_{0}^{2}/\kappa^{2}}}$ the critical density (${\kappa^{2}\equiv 8\pi G=c=\hbar=k_{B}=1}$ by convention). We use $\rho_{\nu}$ for the neutrino energy density summed over the mass eigenstates, and we choose a dynamical parametrization for the scalar energy density $\rho_{\phi}$ in the following Nunes and Lidsey (2004).

In this study, we want to test a possible interaction between neutrino species and dark energy, in a mass-varying neutrino model where active neutrinos are coupled to the scalar field Gu et al. (2003); Fardon et al. (2004); Peccei (2005); Brookfield et al. (2007); Wetterich (2007); Amendola et al. (2008); Ichiki and Keum (2008); Franca et al. (2009). Because to leading order cosmological data are only sensitive to the total neutrino mass Slosar (2006); Font-Ribera et al. (2014), we assume for practical purposes Di Valentino et al. (2018) two massless neutrinos and a massive neutrino nonminimally coupled to the quintessence component. The coupled neutrino has a varying effective mass, which depends on the value of the scalar field and on a dimensionless and constant parameter $\beta$,

$$m_{\nu}(\phi)=m_{\nu,0}\,e^{\beta(\phi-\phi_{0})}\,,$$ (2)

where ${m_{\nu,0}}=\sum m_{\nu}$ is the current neutrino mass and $\phi_{0}$ is the current field value. This coupling can be interpreted as a conformal transformation in the Einstein frame Damour et al. (1990). The massive neutrinos are free falling along the geodesics given by a different metric, $\tilde{g}^{(\nu)}_{\mu\epsilon}$, which is specific to the neutrino sector only, and is related to the gravitational metric, $g_{\mu\epsilon}$, by a conformal transformation,

$$\tilde{g}^{(\nu)}_{\mu\epsilon}=m_{\nu}^{2}\left(\phi\right)g_{\mu\epsilon}\,.$$ (3)

The stress energy tensors of the neutrino fluid and the scalar field are not conserved separately. We have

$$\nabla^{\mu}T_{\mu\epsilon}^{(\nu)} =-\beta(\rho_{\nu}-3p_{\nu})\nabla_{\epsilon}\phi\,,$$ (4)

$$\nabla^{\mu}T_{\mu\epsilon}^{(\phi)} =\beta(\rho_{\nu}-3p_{\nu})\nabla_{\epsilon}\phi\,,$$ (5)

where the subscript $(\nu)$ stands for the stress-energy tensor of the massive neutrino component and $(\phi)$ for that of dark energy, and $p_{\nu}$ is the pressure in the interacting neutrinos. The time component of Eqs. (4) and (5) give the neutrino and scalar field continuity equations, respectively,

$$\dot{\rho}_{\nu}+3H\left(\rho_{\nu}+p_{\nu}\right) =\beta\left(\rho_{\nu}-3p_{\nu}\right)\dot{\phi}\,.$$ (6)

$$\dot{\rho}_{\phi}+3H\left(\rho_{\phi}+p_{\phi}\right) =-\beta\left(\rho_{\nu}-3p_{\nu}\right)\dot{\phi}\,,$$ (7)

where $p_{\phi}$ is the pressure in the field. The extra source terms vanish without interaction, ${\beta=0}$, or if the massive neutrino particles are ultrarelativistic, behaving as traceless radiation.

Regarding the interacting scalar field, we assume that it is homogeneous and canonical with a potential $V(\phi)$ such that ${\rho_{\phi}=-1/2\,g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+V}$, and the quintessence pressure is ${p_{\phi}=\rho_{\phi}-2\,V}$. Equation (7) can also be written as a modified Klein-Gordon equation, which governs the motion of the scalar field, with an additional term due to the neutrino coupling,

$$\ddot{\phi}+3H\dot{\phi}+\frac{dV}{d\phi}=-\beta\left(\rho_{\nu}-3p_{\nu}\right)\,.$$ (8)

To test the model with observations, we adopt a known phenomenological parametrization, first proposed in Ref. Nunes and Lidsey (2004), where the scalar field depends linearly on the number of $e$-folds, ${N\equiv\ln a}$, throughout the cosmological evolution. We introduce a dimensionless constant $\lambda$ for the slope of the scaling:

$$\phi-\phi_{0}=\lambda\ln a\,.$$ (9)

Given the symmetry ${(\phi,\lambda)\mapsto(-\phi,-\lambda)}$, we will limit the present analysis to ${\lambda>0}$. Also, without loss of generality, we choose to set ${\phi_{0}=0}$, which implies ${\phi<0}$. Accordingly, we illustrate in Fig. 1 the cases ${\beta>0}$ (green dashed line) and ${\beta<0}$ (orange dash-dotted line) which correspond to a growing and to a shrinking neutrino mass, respectively, using Eq. (2), while the mass is constant for ${\beta=0}$ (blue solid line).

Figure 1: Neutrino mass evolution in growing ($\beta>0$) and shrinking ($\beta<0$) scenarios for ${m_{\nu,0}=0.06\,\mathrm{eV}}$.

Figure 2: Left panel: evolution of the equations of state of dark energy ($w_{\phi}$) and of massive neutrinos ($w_{\nu}$) for ${m_{\nu,0}=0.06\,\mathrm{eV}}$ and three illustrative values of $\beta$. Right panel: evolution of energy densities for ${m_{\nu,0}=0.06\,\mathrm{eV}}$. The radiation energy density ($\rho_{r}$) includes photons and ultrarelatistic neutrinos.

This simple approach is a powerful alternative to the popular Chevallier–Polarski–Linder (CPL) parametrization Chevallier and Polarski (2001); Linder (2003), since a large variety of dark energy equation of state evolutions can be captured by just one additional parameter Nunes (2004), thereby limiting the degeneracies in Bayesian inferences. The background evolution of the dark energy component allows for an early dark energy component at high redshift, something that is not present in the CPL option. Our scalar parametrization manages to represent a vast class of dark energy models that rely on scaling or tracking solution. It does not, however, allow phantom models, where the scalar field equation of state is below $-1$. These are known to provide weaker constraints on the neutrino mass Di Valentino et al. (2018). It also does not allow an oscillatory evolution of the field around the minimum of an effective potential, such as the model in Ref. Mota et al. (2008).

An additional advantage is that the scalar field potential can be reconstructed analytically following Refs. Nunes and Lidsey (2004); da Fonseca et al. (2022b, a); Barros and da Fonseca (2023). This is done by solving the first-order differential equation (7) to find $\rho_{\phi}$ using the constraint equation (1) and noting that ${\dot{\phi}=\lambda H}$ according to Eq. (9). The potential happens to be a sum of exponential terms,

$$V\left(\phi\right)=A\,e^{-\frac{4}{\lambda}\phi}+B\,e^{-\frac{3}{\lambda}\phi}+C\,e^{\left(-\frac{3}{\lambda}+\beta\right)\phi}+D\,e^{-\lambda\phi}\,,$$ (10)

where the mass scales are given by the following analytical expressions:

$$A =\frac{\lambda^{2}}{4-\lambda^{2}}\Omega_{r}H_{0}^{2}\,,$$ (11)

$$B =\frac{3}{2}\frac{\lambda^{2}}{3-\lambda^{2}}\left(\Omega_{b}+\Omega_{c}\right)H_{0}^{2}\,,$$ (12)

$$C =\frac{3}{2}\frac{\lambda^{2}+2\beta\lambda}{3-\lambda^{2}-\beta\lambda}\Omega_{\nu}H_{0}^{2}\,,$$ (13)

$$D =3\left(1-\frac{\lambda^{2}}{6}\right)H_{0}^{2}\left[1+\frac{4}{\lambda^{2}-4}\Omega_{r}\right.$$ (14)

$$\left.+\frac{3}{\lambda^{2}-3}\left(\Omega_{b}+\Omega_{c}\right)+\frac{3}{\lambda^{2}-3+\beta\lambda}\Omega_{\nu}\right]\,,$$

$\Omega_{r}$ includes both photons and ultrarelativist neutrinos. As for the density parameter $\Omega_{\nu}$, it corresponds to those neutrinos that are nonrelativistic, and we can write Lesgourgues et al. (2013)

$$\Omega_{\nu}h^{2}=\frac{m_{\nu,0}}{93.14\,\mathrm{eV}}\,.$$ (15)

We specifically modify the CLASS code to evolve the scalar field with the potential (10) and the equation of motion (8). The field parametrization is not implemented as such. Independently of the initial conditions, the potential leads to two successive scaling regimes Barreiro et al. (2000), where the field equation of state, shown in the left panel of Fig. 2, first tracks radiation (${w_{\phi}\sim 1/3}$) and then matter (${w_{\phi}\sim 0}$), before being attracted to the late-time acceleration stage (${w_{\phi}\sim-1}$), always present when ${\lambda^{2}<2}$. At this point, the energy density of the field freezes, mimicking a cosmological constant at late time, as shown in the right panel.

We can see from Fig. 2 that the coupling with neutrinos changes $w_{\phi}$ during the matter-dominated era. For growing masses (${\beta>0}$, dotted line), the field equation of state is smaller compared to the uncoupled case (${\beta=0}$, solid line). On the contrary, $w_{\phi}$ is larger when the energy transfer occurs in the opposite direction, i.e., from neutrinos of shrinking mass (${\beta<0}$, dash-dotted line). Correspondingly, Fig. 3 shows that the nonrelativistic neutrinos which receive energy from the scalar field ($\beta>0$) have lower fractional energy density to reach the same present mass than when they give energy ($\beta<0$).

Figure 3: Evolution of the abundances of massive neutrinos and dark energy energy for ${m_{\nu,0}=0.06\,\mathrm{eV}}$.

In the early Universe, neutrinos are held in equilibrium in the primordial plasma by electroweak interactions with charged leptons tightly coupled to photons by electromagnetic processes. As the temperature of the Universe decreases, the interaction rate drops much faster than the Hubble expansion rate. Around ${1\,\mathrm{MeV}}$, at redshift ${z=1/a-1\sim 10^{10}}$, neutrinos cease to interact with the cosmological plasma and flow freely along geodesics. Unlike other particles, the neutrinos are hot relics that decouple while still ultrarelativistic. Therefore their unperturbed phase-space-distribution function maintains to a good approximation the Fermi-Dirac shape for an ultrarelativistic fermion in thermal equilibrium,

$$f_{0}(q)=\frac{1}{e^{q}+1}\,,$$ (16)

neglecting the chemical potential for neutrinos and antineutrinos (i.e., assuming vanishing lepton asymmetry). ${q=ap/T_{\nu,0}}$ is their normalized comoving momentum, $p$ being the physical momentum, and $T_{\nu,0}$ their current temperature. The unperturbed energy density and pressure of the massive neutrino species are thus given by

$$\rho_{\nu} =\frac{1}{\pi^{2}}\left(\frac{T_{\nu,0}}{a}\right)^{4}\int_{0}^{\infty}f_{0}(q)\,dq\,q^{2}\epsilon\,,$$ (17)

$$p_{\nu} =\frac{1}{3\pi^{2}}\left(\frac{T_{\nu,0}}{a}\right)^{4}\int_{0}^{\infty}f_{0}(q)\,dq\,\frac{q^{4}}{\epsilon}\,,$$ (18)

with

$$\epsilon^{2}=q^{2}+\frac{a^{2}\,m_{\nu}^{2}(\phi)}{T_{\nu,0}^{2}}\,,$$ (19)

where $\epsilon$ is the neutrino comoving energy.

The massive neutrinos in the ultrarelativistic regime, that is up to the nonrelativistic transition, behave as standard radiation with their energy density scaling as $a^{-4}$ and ${p_{\nu}=\rho_{\nu}/3}$. In the nonrelativistic regime, the coupling term sources the evolution of both the scalar field and the neutrino fluid, the latter behaving as pressureless matter, ${p_{\nu}=0}$, with its energy density scaling as ${a^{-3}\,e^{\beta\phi}}$. Deep in this regime, the direct integration of Eq. (6) gives

$$\rho_{\nu} =3H_{0}^{2}\Omega_{\nu}\,e^{-3N+\beta\phi}\,,$$ (20)

$$=3H_{0}^{2}\Omega_{\nu}\,a^{-3+\beta\lambda}\,,$$ (21)

given that ${\phi\simeq\lambda N\simeq\lambda\ln a}$ in the parametrization.

From the evolution of their equation of state, ${w_{\nu}=p_{\nu}/\rho_{\nu}}$, in the left panel of Fig. 2, we see that neutrinos with shrinking mass (${\beta<0}$, dash-dotted line) become nonrelativistic earlier than those with growing mass (${\beta>0}$, dotted line), since the latter are lighter in the past and the transition between the two regimes occurs when the average momentum becomes of the order of the mass, ${\langle p\rangle=3.15\,T_{\nu}=m_{\nu}}$ Gerbino and Lattanzi (2018). The redshift $z_{\mathrm{nr}}$ at which the coupled neutrinos become nonrelativistic is given by Gerbino and Lattanzi (2018)

$$1+z_{\mathrm{nr}}\simeq 1900\left(\frac{m_{\nu,0}}{\mathrm{eV}}\right)\,e^{\beta\phi_{\mathrm{nr}}}\,.$$ (22)

According to our choice for the parametrization, ${\phi_{\mathrm{nr}}\simeq-\lambda\ln{(1+z_{\mathrm{nr}})}}$, we get ${e^{\beta\phi_{\mathrm{nr}}}\simeq(1+z_{\mathrm{nr}})^{-\beta\lambda}}$, and thus

$$1+z_{\mathrm{nr}}\simeq\left[1900\left(\frac{m_{\nu,0}}{\mathrm{eV}}\right)\right]^{\frac{1}{1+\beta\lambda}}\,,$$ (23)

where the transition between the two regimes depends on the two extra parameters $\beta$ and $\lambda$ and, as usual, on the current neutrino mass $m_{\nu,0}$.

III Impact of the coupling on perturbations and observables

III.1 Perturbation equations

We apply the theory of linear perturbations Lifshitz (1946) in the synchronous gauge, adopting the usual conventions of Ref. Ma and Bertschinger (1995). In particular, in this section, the scalars $h$ and $\eta$ represent the metric perturbations, and the energy density fluctuation of the cosmological species $i$ is described by the density contrast ${\delta_{i}\equiv\delta\rho_{i}/\bar{\rho_{i}}}$. The overbar denotes the background quantities. The perturbed energy density and pressure of a given species $i$ are ${\delta\rho_{i}\equiv\rho_{i}(\bm{x},\tau)-\bar{\rho}_{i}(\tau)}$ and ${\delta p_{i}\equiv p_{i}(\bm{x},\tau)-\bar{p}_{i}(\tau)}$, respectively. In this section, the dot denotes derivation with respect to conformal time $d\tau\equiv dt/a$.

Figure 4: Left panel: evolution of the massive neutrino perturbations at the scale ${k=0.1\,h/\mathrm{Mpc}}$ for ${\lambda=0.1}$ and ${m_{\nu,0}=0.06\,\mathrm{eV}}$. Right panel: suppression of the linear matter power spectrum at $z=0$ with respect to the $\Lambda$CDM model.

The perturbed energy density and pressure of the interacting neutrinos have been derived in previous studies (see, e.g., Ichiki and Keum (2008); Brookfield et al. (2007); Franca et al. (2009)):

$$\delta\rho_{\nu} =4\pi\left(\frac{T_{\nu,0}}{a}\right)^{4}\int_{0}^{\infty}f_{0}(q)\,dq\,q^{2}\left(\epsilon\Psi+\beta\delta\phi\frac{m^{2}_{\nu}a^{2}}{\epsilon}\right)\,,$$ (24)

$$\delta p_{\nu} =\frac{4\pi}{3}\left(\frac{T_{\nu,0}}{a}\right)^{4}\int_{0}^{\infty}f_{0}(q)\,dq\,\frac{q^{4}}{\epsilon}\left(\Psi-\beta\delta\phi\frac{m^{2}_{\nu}a^{2}}{\epsilon^{2}}\right)\,,$$ (25)

where ${\delta\phi\equiv\phi(\bm{x},t)-\bar{\phi}(t)}$ denotes the fluctuations of the coupled scalar field, and ${\lvert\Psi\rvert\ll 1}$ is the relative perturbation to the neutrino phase-space distribution,

$$\Psi\left(x^{i},q,n_{j},\tau\right)=\frac{f\left(x^{i},q,n_{j},\tau\right)}{f_{0}\left(q\right)}-1\,,$$ (26)

at first order in the perturbations, $f$ being the perturbed distribution function, and ${n^{i}n_{i}=1}$. The dipole equation for the neutrino hierarchy is affected by the interaction, and the corresponding system of infinite equations becomes the following in Fourier space:

$$\dot{\Psi}_{0} =-\frac{qk}{\epsilon}\Psi_{1}+\frac{\dot{h}}{6}\frac{d\ln f_{0}}{d\ln q}\,,$$ (27)

$$\dot{\Psi}_{1} =\frac{qk}{3\epsilon}\left(\Psi_{0}-2\Psi_{2}\right)-\frac{qk}{3\epsilon}\beta\delta\phi\frac{a^{2}m_{\nu}^{2}}{q^{2}}\frac{d\ln f_{0}}{d\ln q}\,,$$ (28)

$$\dot{\Psi}_{2} =\frac{qk}{5\epsilon}\left(2\Psi_{1}-3\Psi_{3}\right)-\left(\frac{\dot{h}}{15}+\frac{2\dot{\eta}}{5}\right)\frac{d\ln f_{0}}{d\ln q}\,,$$ (29)

$$\dot{\Psi}_{\ell\geq 3} =\frac{qk}{\left(2\ell+1\right)\epsilon}\left[\ell\Psi_{\ell-1}-\left(\ell+1\right)\Psi_{\ell+1}\right]\,,$$ (30)

where $\Psi_{\ell}$ is the $\ell$th Legendre component of the series expansion in multipole space of the perturbation $\Psi$. We modified the noncold dark matter part of the CLASS code Lesgourgues and Tram (2011) to evolve the perturbation equations of the MaVaN model.

In the fluid approximation, on sub-Hubble scales, the Boltzmann hierarchy in the momentum grid is cut at ${l_{\mathrm{max}}=2}$ and the continuity equation reads

$$\dot{\delta}_{\nu}= \,3\left(aH+\beta\dot{\phi}\right)\left(w_{\nu}-\frac{\delta p_{\nu}}{\delta\rho_{\nu}}\right)\delta_{\nu}$$

$$-\left(1+w_{\nu}\right)\left(\theta_{\nu}+\frac{\dot{h}}{2}\right)+\beta\left(1-3w_{\nu}\right)\delta\dot{\phi}\,,$$ (31)

where $\theta_{\nu}$ is the neutrino flux divergence. Moreover, the Euler equation is

$$\dot{\theta}_{\nu}= -aH\left(1-3w_{\nu}\right)\theta_{\nu}-\frac{\dot{w}_{\nu}}{1+w_{\nu}}\theta_{\nu}+\frac{\delta p_{\nu}/\delta\rho_{\nu}}{1+w_{\nu}}k^{2}\delta_{\nu}$$

$$+\beta\frac{1-3w_{\nu}}{1+w_{\nu}}k^{2}\delta\phi-\beta\left(1-3w_{\nu}\right)\dot{\phi}\,\theta_{\nu}-k^{2}\sigma_{\nu}\,,$$ (32)

where the neutrino anisotropic stress $\sigma_{\nu}$ Ma and Bertschinger (1995) is not changed by the coupling. We have adjusted the fluid approximation equations of the noncold dark matter in the CLASS code accordingly.

Deep in the nonrelativistic regime, when ${w_{\nu}=0}$, the ratio ${q/\epsilon}$ vanishes asymptotically and the pressure perturbations in the neutrino fluid, as well as the shear stress, become negligible with respect to density perturbations. The continuity and Euler equations are analogous to those of the coupled cold dark matter model Amendola (2000); da Fonseca et al. (2022b),

$$\dot{\delta}_{\nu}+\theta_{\nu}+\frac{\dot{h}}{2} =\beta\delta\dot{\phi}\,,$$ (33)

$$\dot{\theta}_{\nu}+aH\theta_{\nu} =\beta\left(k^{2}\delta\phi-\dot{\phi}\theta_{\nu}\right)\,.$$ (34)

For the coupled scalar field, the equation of motion of the fluctuations is the following,

$$\delta\ddot{\phi}+2aH\delta\dot{\phi}+\left(k^{2}+a^{2}\frac{d^{2}V}{d\phi^{2}}\right)\delta\phi+\frac{1}{2}\dot{h}\dot{\phi}$$

$$=-a^{2}\beta\left(\delta\rho_{\nu}-3\delta p_{\nu}\right)\,.$$ (35)

As in the background, we evolve the field perturbations with the potential through the above equation in our version of the CLASS code.

III.2 Effects on the matter power spectrum

There are three main stages of the evolution of the neutrino density contrast affected by the coupling. During the radiation-dominated era, when the neutrinos are decoupled from the thermal bath but still relativistic, their perturbations grow as radiation. Later, the neutrinos become nonrelativistic and cluster in the gravitational potential wells of cold dark matter, which is the dominant cosmological component. However, below their free-streaming scale, they do not cluster like cold dark matter Lesgourgues and Pastor (2006). Neutrino free streaming dampens the neutrino fluctuations up to a critical scale depending on the neutrino mass, and giving the oscillatory pattern seen in the left panel of Fig. 4. The free-streaming wave number of Fourier mode reaches a minimum at the nonrelativistic transition, given by Gerbino and Lattanzi (2018)

$$k_{\mathrm{fs}}\simeq 0.018\,\Omega_{m}^{1/2}\left(\frac{m_{\nu,0}}{\mathrm{eV}}\right)^{1/2}e^{\beta\phi_{\mathrm{nr}}/2}\,h\,\mathrm{Mpc}^{-1}\,,$$ (36)

during matter or dark energy domination. Or equivalently, using Eqs. (22) and (23), we get

$$k_{\mathrm{fs}}\simeq 4.1\cdot 10^{-4}\,\Omega_{m}^{1/2}\left(1900\frac{m_{\nu,0}}{\mathrm{eV}}\right)^{\frac{1}{2(1+\beta\lambda)}}\,h\,\mathrm{Mpc}^{-1}\,,$$ (37)

for our particular scalar field parametrization. Above the free-streaming length, the neutrino fluctuations grow unhindered. For growing neutrino masses (${\beta>0}$, green dashed line) the free-streaming scale in Eq.(37) is larger and the growth of the fluctuations is delayed with respect to shrinking neutrino masses (${\beta<0}$, orange dash-dotted line).

Figure 5: Left panel: comparison of the power spectra of the CMB anisotropies with respect to $\Lambda$CDM, for ${\lambda=0.1}$ and ${m_{\nu,0}=0.06\,\mathrm{eV}}$, on large and intermediate scales. Right panel: same as left, on medium and small scales.

Furthermore, the dependence of the neutrino mass on $\beta$ changes the fraction of matter whose fluctuations do not grow like cold dark matter at a given scale. The neutrinos do not contribute to the creation of potential wells below the free streaming scale, and all structure formation is damped because the gravitational wells are not as deep as they would be in the presence of only nonrelativistic matter.

The effect of the coupling on the linear matter power spectrum, which is proportional to the variance of the density fluctuations (which tells how large they are at a given scale), can be seen on small scales, at large wave numbers $k$. The right panel of Fig. 4 shows the residual plot between the MaVaN scenario and the standard $\Lambda$ cold dark matter (CDM) model with massless neutrinos, normalized to the power spectrum of the latter. At scales smaller than the critical scale ($k\gtrsim 10^{-3}$), below which neutrinos do not cluster, the perturbations in the neutrino fluid are completely damped by free streaming and do not contribute to matter perturbations. In this respect, we see with the gray solid line that the ${\nu\Lambda\mathrm{CDM}}$ model suppresses power at these scales compared to the $\Lambda$CDM model with massless neutrinos.

Moreover, the non-negligible fraction of dark energy itself (${\lambda\neq 0}$ and ${\beta=0}$, blue solid line) further reduces the growth of the fluctuations during matter dominance, leading to more power suppression. On the other hand, the matter power spectrum at small scales also depends on how large the neutrino mass was in the past. Growing neutrino masses (${\beta>0}$, green dashed line) reduce the power suppression caused by the scalar field, while shrinking neutrino masses increase the suppression (${\beta<0}$, orange dash-dotted line).

On the contrary, it can be seen that $\beta$ has little influence at large scales ($k\lesssim 10^{-3}$), since the neutrinos cluster independently of their mass.

III.3 Effects on the CMB anisotropies power spectrum

The temperature power spectrum of the unlensed CMB is also affected by the coupling through background and perturbation effects. After their nonrelativistic transition, the massive neutrinos modify the evolution of the cosmological background and the redshift of the matter-dark energy equality, with respect to an uncoupled neutrino universe, affecting the late time integrated Sachs-Wolfe effect. As shown in the left panel of Fig. 5, at large angular scales (lower multipoles ${\ell\lesssim 20}$) reminiscent of the scale-invariant primordial spectrum of the inflationary perturbations, the additional power produced by the time evolution of the gravitational potentials, as dark energy begins to dominate, decreases with negative couplings (orange dash-dotted line).

On the intermediate scales (${20\lesssim\ell\lesssim 500}$), dominated by the imprints of the acoustic oscillations in the baryon-photon fluid, it is again the effect of the coupling on the evolution of the gravitational potential that affects the first peak. In particular, $\beta$ alters the time required for the gravitational potential to stabilize at a nearly constant value after photon decoupling. Keeping the present value of the Hubble parameter $h$ unchanged, a higher neutrino mass at the time of recombination (${\beta<0}$, orange dash-dotted line) increases the height of the first peak given by the early integrated Sachs-Wolfe effect compared to the $\Lambda$CDM cases. Conversely, the amplitude decreases for a lower neutrino mass at recombination (${\beta>0}$, green dashed line).

The amplitude of the first peak is also affected by the uncoupled scalar field alone (${\beta=0}$, blue solid line). It gives a non-negligible contribution to early dark energy, which reduces the fractional energy density of matter at the time of decoupling da Fonseca et al. (2022a, b). The amplitude increases with the value of $\lambda$.

Moreover, for a given physical density of cold dark matter and baryons, the angle subtended by the sound horizon at recombination, $\theta_{s}$, which determines the spacing between the peaks and in particular the position of the first one, is larger for growing neutrino masses ($\beta>0$). The corresponding shift of the position of the CMB peaks towards larger scales can be compensated by the Hubble constant. Indeed, decreasing $H_{0}$ increases the comoving angular diameter distance from the CMB surface and shifts the peaks towards smaller scales.

At small scales (${\ell\gtrsim 500}$), in the right panel of Fig. 5, of the order of the photon mean free path at the time of recombination, a positive coupling (orange dash-dotted line) has the opposite effect of $\lambda$ (${\beta=0}$, blue solid curve) on the exponential damping of the CMB peak structure.

III.4 Effects on the CMB lensing potential

Because the free-streaming neutrinos erase the density perturbations, they affect the CMB light that is distorted by the gravitational lensing caused by the intervening matter distribution between us and the last scattering surface Lesgourgues et al. (2006). The neutrinos reduce the CMB lensing potential, which is a measure of the integral of the gravitational potentials along the line of sight between the recombination time and the present time. The effect of the weak lensing is to smooth the power spectrum of the CMB temperature anisotropies on small scales. Note in Fig. 6 that since the effect is proportional to the energy density of the neutrinos, it can constrain their mass, whose cosmological evolution is controlled by the two parameters $\lambda$ and $\beta$. For example, if the neutrino mass had been too high in the recent past, we would have had less lensing than we observe. The suppression already caused by the scalar field (${\beta=0}$, blue solid curve) is either enhanced by shrinking neutrino masses ($\beta<0$, orange dash-dotted line) or compensated by growing neutrino masses ($\beta>0$, green dashed line).

Figure 6: Upper panel: power spectra of the CMB lensing potential of the $\Lambda$CDM and $\nu\Lambda\mathrm{CDM}$ models, and the coupled scenario with $\lambda=0.1$ and $\sum m_{\nu}=0.06\,\mathrm{eV}$. Lower panel: percentage deviation between models with respect to $\Lambda$CDM.

It is worth noting that, in contrast to the model-independent parametrization for the neutrino mass variation studied in Ref. Franca et al. (2009), we do not find instabilities on large scales in our model Bjaelde et al. (2008), which would be triggered by large coupling values causing the neutrino perturbations to grow rapidly on the largest scales observable.

IV Parameter estimation

IV.1 Methodology

We choose to test the model with the 2018 Planck satellite Planck Collaboration et al. (2011) temperature and polarization measurements of the CMB at the time of last scattering, combined with the CMB lensing potential as a probe of the distribution and evolution of large-scale structure in the late Universe. We also add BAO distance and expansion rate measurements in galaxy clusters ”Eisenstein and others” (2005); Cole et al. (2005) covering different redshift ranges. They constrain the model at the background level Percival et al. (2007) to help break the geometric degeneracy between $H_{0}$ and $\sum m_{\nu}$.

With respect to the Planck observations, the Bayesian inference is performed using the CMB lensing potential power spectrum Aghanim et al. (2020b), as well as the cross-correlation likelihoods of the temperature and polarization anisotropies, $TT$, $EE$, and $TE$, respectively Aghanim et al. (2020c). We reduce the dimensional parameter space by using the lite version of the Planck likelihood on the smaller angular scales, which marginalizes the foreground and instrumental effects, leaving the absolute calibration $A_{\mathrm{planck}}$ as the only nuisance parameter. As for the background probe, we use a compilation of measurements that have detected the BAO peak at different angular separations with galaxy samples at different redshifts. The three data points are from the Sloan Digital Sky Survey DR7 Main Galaxy Sample Ross et al. (2015), the Sloan Digital Sky Survey DR9 release Ahn et al. (2012), and the 6dF Galaxy Survey Beutler et al. (2011), respectively. We will refer to our dataset as Plk18+BAO in the rest of the document.

Our dedicated version of the CLASS code is used to compute the observables confronting the actual data, on the basis of the sampled parameters $\{\lambda,\beta,\sum m_{\nu}\}$ along with the six standard minimal baseline parameters $\{\Omega_{b}h^{2},\Omega_{c}h^{2},\theta_{s},A_{s},n_{s},\tau_{\mathrm{reio}}\}$. We use $\theta_{s}$ instead of $h$ because it is the quantity measured in the CMB observations. $A_{s}$ is the power of the primordial curvature perturbations normalized at the pivot scale ${k=0.05\,\mathrm{Mpc}^{-1}}$, and $n_{s}$ is the power-law index of the scalar spectrum. $\tau_{\mathrm{reio}}$ is the optical depth to reionization, which gives the fraction of photons rescattered by the new population of free electrons from the ionization of neutral hydrogen produced by the light of the first stars. The settings in CLASS are such that the effective number of neutrino families is $N_{\mathrm{eff}}=3.044$ Froustey et al. (2020); Bennett et al. (2021); Akita and Yamaguchi (2020), and the neutrino fractional energy density satisfies Eq. (17).

The likelihoods are minimized by the Monte Carlo code of the MONTEPYTHON parameter estimation package Audren et al. (2013); Brinckmann and Lesgourgues (2019), which samples the parameter space and the posterior probability distributions using a Metropolis-Hastings algorithm with the flat priors specified in Table 1. The size of the prior intervals is sufficient for the corresponding posteriors to fall exponentially within them, except for the parameter $\beta$ as we will discuss later.

Table 1: Flat priors for the sampled model and cosmological parameters.

Parameter	Prior
$\lambda$	$[0,1.4]$
$\beta$	$[0,100]$
$\sum m_{\nu}$	$[0,2]$
$\Omega_{b}h^{2}$	$[0.005,0.1]$
$\Omega_{c}h^{2}$	$[0.01,0.99]$
$100\theta{}_{s}$	$[1,2]$
$\ln 10^{10}A_{s}$	$[2.7,4]$
$n_{s}$	$[0.5,1.5]$
$\tau{}_{\mathrm{reio}}$	$[0.04,0.8]$

We use the GETDIST analyzer Lewis (2019) to obtain the constraints from the Markov chains, and plot the confidence contours and the marginalized posterior distributions. In addition to the sampled parameters, we also infer constraints on derived late-Universe parameters : $H_{0}$, $\Omega_{m}$, $\sigma_{8}$ (which measures the amplitude of matter fluctuation on ${8\,h^{-1}\mathrm{Mpc}}$ comoving scale), and the degeneracy parameter ${S_{8}\equiv\sigma_{8}\sqrt{\Omega_{m}/0.3}}$.

IV.2 Main results

The results of the likelihood analysis are summarized in Table 2, which lists the constraints on the cosmological parameters obtained with the Plk18+BAO dataset for both the uncoupled ($\beta=0$) scalar field model and the coupled MaVaN model. They lead to several remarkable conclusions.

Table 2: Cosmological constraints (95% limits for $\sum m_{\nu}$, mean and $68\%$ limits for the others) for the uncoupled (${\beta=0}$) and MaVaN ($\beta$ free) models.

Parameter	$\beta=0$	$\hskip 28.45274pt\beta$ free
$\lambda$	$<0.0472$	$\hskip 28.45274pt0.059^{+0.029}_{-0.037}$
$\sum m_{\nu}$	$<0.127$	$\hskip 28.45274pt<0.724$
$\Omega_{b}h^{2}$	$0.02241\pm 0.00014$	$\hskip 28.45274pt0.02242\pm 0.00014$
$\Omega_{c}h^{2}$	$0.1199\pm 0.0011$	$\hskip 28.45274pt0.1200\pm 0.0011$
$100\theta{}_{s}$	$1.04186\pm 0.00030$	$\hskip 28.45274pt1.04179\pm 0.00031$
$\ln 10^{10}A_{s}$	$3.050^{+0.014}_{-0.016}$	$\hskip 28.45274pt3.049^{+0.014}_{-0.016}$
$n_{s}$	$0.9667\pm 0.0039$	$\hskip 28.45274pt0.9669\pm 0.0039$
$\tau{}_{\mathrm{reio}}$	$0.0568^{+0.0068}_{-0.0082}$	$\hskip 28.45274pt0.0564^{+0.0068}_{-0.0080}$
$H_{0}$	$67.59^{+0.69}_{-0.57}$	$\hskip 28.45274pt67.61^{+0.65}_{-0.56}$
$\Omega_{m}$	$0.3127^{+0.0074}_{-0.0091}$	$\hskip 28.45274pt0.3176^{+0.0084}_{-0.012}$
$\sigma_{8}$	$0.813^{+0.012}_{-0.0079}$	$\hskip 28.45274pt0.804^{+0.019}_{-0.010}$
$S_{8}$	$0.830\pm 0.012$	$\hskip 28.45274pt0.827\pm 0.015$
$A_{\mathrm{planck}}$	$1.0008\pm 0.0025$	$\hskip 28.45274pt1.0007\pm 0.0025$

First, in the uncoupled scalar field model (${\beta=0}$), the dynamical dark energy component is constrained to be close to a cosmological constant, ${\lambda<0.05}$ (68% C.L.) and ${\lambda<0.09}$ (95% C.L.), with a clear preference for a vanishing scalar field parameter. This is the result of the extreme constraining power of the CMB data on the fraction of early dark energy at the last scattering surface Gómez-Valent et al. (2021, 2022). In our scalar field parametrization, the fraction of dark energy (during the scaling regime with matter) is given by $\Omega_{\phi}=\lambda^{2}/3$ da Fonseca et al. (2022b), which implies $\Omega_{\phi}<0.27\%$ ($95\%\,\mathrm{CL}$) at recombination. It is therefore not surprising that the corresponding upper limit ${\sum m_{\nu}<0.127\,\mathrm{eV}}$ (95% C.L.) that we find is in agreement with recent literature findings for the concordance model based on comparable datasets. For example, including the latest release of BAO eBOSS data in combination with Planck temperature, polarization, and lensing measurements gives an upper limit of ${0.129\,\mathrm{eV}}$ (95% C.L.) for the uncoupled ${\nu\Lambda\mathrm{CDM}}$ model Alam et al. (2021).

Second, we find that the parameter $\beta$ of the MaVaN model, with our choice of scalar field parametrization, is not constrained by the cosmological observations Plk18+BAO, as shown in Fig. 7 (and the Appendix for additional contour and probability distribution plots). However, the interaction does relax the upper bound on the current neutrino mass sum, raising it significantly by almost a factor of six to ${\sum m_{\nu}<0.724\,\mathrm{eV}}$ (95% C.L.). Positive couplings allow for the possibility of heavier neutrinos today. Despite their current higher mass, the influence of neutrinos with growing mass on the cosmological expansion and perturbations over time is tempered by the fact that they were lighter in the past and gradually gained energy from the scalar field. Moreover, for the MaVaN model, the data favor a nonvanishing scalar field parameter, ${\lambda=0.059^{+0.029}_{-0.037}}$ (68% C.L.), at the one-sigma level.

Figure 7: Constraints obtained with the Plk18+BAO dataset on the additional parameters. Probability distributions and 2D marginalized contours (68% and 95% C.L.).

Third, with respect to the standard cosmological parameters, the main differences between the MaVaN model and the uncoupled one are in the matter fluctuation amplitude, $\sigma_{8}$, and the matter density parameter, $\Omega_{m}$. Including the coupling parameter degrades the error bars while virtually preserving the central values, which are consistent with the best fit of the base-$\Lambda$CDM cosmology predicted by the 2018 Planck data Aghanim et al. (2020a). This is particularly the case for the derived parameter $\sigma_{8}$, whose precision decreases by about $50\%$, from $0.813^{+0.012}_{-0.0079}$ to $0.804^{+0.019}_{-0.010}$ (68% C.L.), mainly towards lower values, which could tend to reduce the tension with late cosmological probes Amon et al. (2022); Abbott et al. (2023); Secco et al. (2022). However, the confidence interval of $\Omega_{m}$ increases by $22\%$ in the opposite direction due to the added fractional energy density of heavier neutrinos, from $0.3126^{+0.0074}_{-0.0090}$ to $0.3176^{+0.0084}_{-0.012}$ (68% C.L.). As shown in Appendix Appendix: Parameter constraints for the MaVaN and the uncoupled scalar field models, the posterior distribution for the degeneracy parameter $S_{8}$ is only slightly changed by the coupling parameter. As for $H_{0}$, its posterior distribution is hardly modified by the interaction, confirming the inability of the model to resolve the Hubble tension. The fractional dark energy in the early times is far less than the minimum $10\%$ required by neutrino-assisted early dark energy models to improve the $H_{0}$ tension Poulin et al. (2019); de Souza and Rosenfeld (2023).

Finally, to assess the influence of the coupling on the posterior $\sum m_{\nu}$, we repeat the statistical analysis fixing different values for $\beta$. The results are summarized in Table 3 and Fig. 8. For a mass-shrinking neutrino scenario (${\beta=-1}$), the constraint on today’s mass is tighter (${\sum m_{\nu}<0.102\,\mathrm{eV}}$) than without interaction. Large masses in the past are disfavored by the cosmological data. On the contrary, the limit on today’s neutrino masses relaxes as the value of the positive coupling increases. The highest limit is for $\beta=100$, i.e., for the edge of the prior, with $\sum m_{\nu}<1.130\,\mathrm{eV}$ already exceeding the limit obtained by local experiments Aker et al. (2019). Moreover, there seems to be a plateau around $\beta=20$ where the present mass stabilizes in the ($\beta,\sum m_{\nu}$) plane of the parameter space (see also Fig. 7). In addition, as seen earlier, the energy transferred by the scalar field causes the model to deviate slightly from a constant dark energy component. Dark energy moves away from a cosmological constant as the strength of the coupling increases. However, the value of the scalar field parameter is tightly constrained by the CMB data, which does not allow $\lambda$ to be too large, while it could theoretically be as large as $\sqrt{2}$.

Table 3: Constraints on present neutrino mass in eV and on the dark energy parameter $\lambda$, for different coupling values $\beta$.

$\beta$	$\sum m_{\nu}\,(95\%\,\mathrm{C.L.})$	$\lambda\,(68\%\,\mathrm{C.L.})$
-1	$<0.102$	$<0.0497$
0	$<0.127$	$<0.0472$
1	$<0.147$	$<0.0493$
5	$<0.355$	$0.046^{+0.016}_{-0.043}$
10	$<0.580$	$0.050^{+0.022}_{-0.038}$
20	$<0.542$	$0.052^{+0.024}_{-0.038}$
30	$<0.513$	$0.053^{+0.025}_{-0.039}$
50	$<0.598$	$0.059^{+0.028}_{-0.037}$
75	$<0.822$	$0.064^{+0.031}_{-0.035}$
100	$<1.130$	$0.067\pm 0.030$
free	$<0.724$	$0.059^{+0.029}_{-0.037}$

Figure 8: Constraints obtained for fixed coupling parameters. Probability distributions and 2D marginalized contours (68% and 95% C.L.).

We note in Fig. 8 that the $\sum m_{\nu}$ positive and negative correlations with $\Omega_{m}$ and $\sigma_{8}$, respectively, are preserved by the coupling. The errors on $S_{8}$ are degraded to lower values as the strength of the interaction increases. It is also worth noting that the degeneracy in the ($\sum m_{\nu},H_{0}$) plane, along which the position of the first CMB peak is held unchanged, decreases significantly with increasing value of the parameter $\beta$, weakening the otherwise strong correlation between the neutrino mass and the Hubble constant. This is probably due to the fact that the posterior for $\lambda$ becomes larger with increasing $\beta$. The corresponding non-negligible fraction of dark energy at the time of recombination reduces the sound horizon at decoupling, thus competing with the geometric degeneracy $\sum m_{\nu}-H_{0}$.

V Discussion

In this study, we investigated a model of mass-varying neutrinos in which the mass depends on the value of a scalar field representing the dark energy component. The originality of our work lies in the choice of the quintessence parametrization, which limits the number of free parameters with respect to other approaches such as the CPL parametrization or arbitrary choices of scalar field potentials. While the parametrization is purely phenomenological, it leads to the analytical reconstruction of the potential as a sum of exponential terms, which conveniently endows dark energy with scaling properties. We have introduced two additional free parameters with respect to $\Lambda$CDM, $\beta$ for the coupling strength between the two sectors in Eq. (2) and $\lambda$ for the linear evolution of the field in Eq. (9).

We confirmed the conjecture that a growing neutrino mass scenario relaxes the upper bound on the current neutrino mass derived in the framework of the baseline $\nu\Lambda$CDM model. This is a result that we expect can be replicated for any quintessence model that manifests scaling or tracking behavior. Our goal was to complement previous work, such as the study in Ref. Brookfield et al. (2007). In the latter, one of the models considered is a mass-varying neutrino theory with an identical coupling to the scalar field, which is a classical choice made throughout the literature, and a different potential in the form of a single exponential. In contrast to our case, their quintessence model shows no tracking behavior and allows only shrinking mass scenarios. The authors also found that astronomical observations do not provide strong constraints on the coupling parameter. To constrain the coupling, they fixed the neutrino mass to values higher than ${0.1\,\mathrm{eV}}$, assuming that such a significant mass could be independently confirmed. Alternately, in our likelihood analysis, we set different values of the coupling to constrain the current neutrino mass.

Following the analysis of the model at background level, we evaluated the sensitivity of several observables (matter and CMB power spectra, and CMB lensing potential) to the coupling, using a version of the Boltzmann code CLASS that we adapted to the considered model. We found that growing neutrino masses lead to less matter power suppression than predicted by the presence of a noninteracting scalar field. The coupling also affects the shape of the CMB power spectrum at different scales, in particular through the integrated Sachs-Wolfe effect, in parallel with the effect of the scalar field parameter itself. The CMB lensing potential is sensitive to the interaction. Growing neutrino masses can compensate for the reduction in lensing potential caused by the quintessence fluid. It is therefore theoretically possible to obtain constraints on a putative interaction between the neutrino sector and a dynamical dark energy component.

We performed Bayesian inferences to estimate the parameters of the model using observations of CMB temperature and polarization anisotropies and lensing from the Planck satellite, together with BAO measurements. This dataset allows us to combine cosmological probes covering the early Universe at the time of the last scattering, the formation of large-scale structures in the late universe, and the cosmic expansion history. The Markov chains were generated by the MONTEPYTHON package using a Metropolis-Hastings algorithm to explore the parameter space, which contains nine free parameters, $\{\lambda,\beta,\sum m_{\nu},\Omega_{b}h^{2},\Omega_{c}h^{2},\theta_{s},A_{s},n_{s},\tau_{\mathrm{reio}}\}$, plus one nuisance parameter $A_{\mathrm{planck}}$. We also derived constraints on four late Universe parameters: $H_{0}$, $\sigma_{8}$, $\Omega_{m}$ and $S_{8}$.

The main conclusion is that the model relaxes the upper bound on the neutrino mass, from ${\sum m_{\nu}<0.13\,\mathrm{eV}}$ in $\nu\Lambda$CDM to ${\sum m_{\nu}<0.72\,\mathrm{eV}}$ in our model (95% C.L.). The dynamical property of dark energy is enhanced because the scalar field looses energy to the neutrino, whose mass grows with time. With a non-vanishing scalar field parameter, ${\lambda=0.059^{+0.028}_{-0.036}}$ (68% C.L.), the data disfavors a static field. This is in contrast with previous work where CMB data constraining a scalar field interacting with dark matter model favors $\Lambda$CDM da Fonseca et al. (2022b). The coupling $\beta$ is, however, poorly constrained by the observations of the selected dataset. As for the cosmological parameters of the standard model, it is mostly the confidence intervals for $\sigma_{8}$ and $\Omega_{m}$ that are broadened by the interaction, but in a way that weakly affects the posterior distribution of the degeneracy parameter $S_{8}$. The posterior distribution of $H_{0}$ is also almost unchanged, as expected. The scalar field parameter is extremely constrained and the fraction of dark energy at the last scattering is too small to make any impact in resolving the Hubble tension.

As regards forthcoming work, it would be worthwhile to complement the CMB constraints that the Planck data place on the parameters of our model with precise measurements of the anisotropies at larger multipoles ($l\gtrsim 3000$). These small angular scales measured with sufficient accuracy can reveal coupling signatures, especially when the strength of the interaction is small Brax et al. (2023a). For example, an indication of nonvanishing coupling between neutrinos and dark matter at the one-sigma level has been found using the Atacama Cosmology Telescope’s high-multipole observations of CMB temperature and polarization anisotropies Brax et al. (2023b). In addition, alternative CMB data on the lensing power spectrum could also be used to further constrain the neutrino-scalar field interaction scenario affecting structure growth Qu et al. (2023). As for the late-time Universe probes of large-scale structure, it would be adequate to use the Kilo-Degree Survey (KiDS) weak lensing observations Hildebrandt et al. (2017) to test the MaVaN model, like in the coupled dark matter case da Fonseca et al. (2022b).

In the future, combinations of CMB experiments (such as Simons Observatory Ade et al. (2019), CMB-S4 Abazajian et al. (2016), or LiteBird Hazumi et al. (2019)) with large-scale structure experiments [for example Dark Energy Spectroscopic Instrument (DESI) Collaboration et al. (2016), Euclid Laureijs et al. (2011), Large Synoptic Survey Telescope (LSST) Collaboration et al. (2009), or Square Kilometre Array (SKA) Bacon et al. (2020)] are expected to significantly improve the robustness of the cosmological neutrino mass measurements Brinckmann et al. (2019); Archidiacono et al. (2017), possibly allowing us to discriminate between models of interacting dark energy.

Appendix: Parameter constraints for the MaVaN and the uncoupled scalar field models

In this appendix we provide the triangular plots (Fig. 9) of the analysis made with the Plk18+BAO data for the MaVaN model ($\beta$ free) and the uncoupled case ($\beta=0$), i.e., no interaction between the quintessence component and the neutrino fluid.

Figure 9: Probability distributions and 2D marginalized contours (68% and 95% C.L.) obtained with the Plk18+BAO dataset.