Origin of cosmological neutrino mass bounds: background versus perturbations

Toni Bertólez-Martı́nez ,1, ∗ Ivan Esteban ,2, 3, † Rasmi Hajjar ,4, ‡ Olga Mena ,4, § and Jordi Salvado 1, ¶
1
Departament de Fı́sica Quàntica i Astrofı́sica and Institut de Ciències del Cosmos,
Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Spain
2
Department of Physics, University of the Basque Country UPV/EHU, PO Box 644, 48080 Bilbao, Spain
3
EHU Quantum Center, University of the Basque Country UPV/EHU
4
Instituto de Fı́sica Corpuscular (IFIC), University of Valencia-CSIC,
Parc Cientı́fic UV, c/ Catedrático José Beltrán 2, E-46980 Paterna, Spain
(Dated: November 25, 2024)
The cosmological upper bound on the total neutrino mass is the dominant limit on this fundamen-
tal parameter. Recent observations—soon to be improved—have strongly tightened it, approaching
arXiv:2411.14524v1 [astro-ph.CO] 21 Nov 2024

the lower limit set by oscillation data. Understanding its physical origin, robustness, and model-
independence becomes pressing. Here, we explicitly separate for the first time the two distinct cos-
mological neutrino-mass effects: the impact on background evolution, related to the energy in neu-
trino masses; and the “kinematic” impact on perturbations, related to neutrino free-streaming. We
scrutinize how they affect CMB P anisotropies, introducing two effective masses enclosing background
( mBackg. ) and perturbations ( mPert.
P
ν ν ) effects. We analyze CMB data, finding that the neutrino-
mass bound is mostly a background measurement,P i.e., how the neutrino energyP density evolves
with time. The bound on the “kinematic” variable mPert.
ν is largely relaxed, mPert.
ν < 0.8 eV.
This work thus adds clarity to the physical origin of the cosmological neutrino-mass bound, which is
mostly a measurement of the neutrino equation of state, providing also hints to evade such a bound.

I. INTRODUCTION energy density. At temperatures ≲ MeV, their weak-

interaction rate with the primordial plasma falls below
Neutrinos are hot thermal relics, permeating the Uni- the expansion rate of the Universe. Then, they decouple
verse and being the most abundant particles after Cos- and stream freely, with their momentum being redshifted
mic Microwave Background (CMB) photons. In particle by the expansion. When neutrino momenta fall below
physics, neutrino masses have meant the first departure their mass, neutrinos move at non-relativistic speeds,
from the Standard Model (SM), where no gauge-invariant contributing to the Universe’s matter energy density.
renormalizable neutrino mass term can be written. Yet Due to their abundance, neutrinos leave a measurable
neutrino oscillation experiments have robustly measured imprint on cosmological observables. This leads to up-
two squared-mass differences, |∆m231 | ≡ |m23 − m21 | ≃ per
P bounds on the total neutrino mass that are around
2.5 · 10−3 eV2 and ∆m221 ≡ m22 − m21 ≃ 7.5 · 10−5 eV2 [1– mν ≲ 0.05–0.3 eV [10–13], potentially raising a ten-
3], confirming that neutrinos are massive and that the sion between oscillation results and cosmological observa-
SM must be extended. Two mass orderings are possible, tions. As cosmological datasets grow, the neutrino mass
normal (NO, ∆m231 > 0) and inverted (IO, ∆m231 < 0), will soon be either measured or a bound in clear conflict
and oscillation experiments aim to determine it [4–8]. with oscillation experiments will be placed [14–17].
The absolute neutrino mass scale, however, remains Cosmological inferences, however, are indirect. The
unknown, as oscillations are only sensitive
P to mass differ- cosmological neutrino-mass bound is known to be corre-
ences. Oscillation
P results constrain m ν > 0.06 eV for lated with, among others, the equation of state of dark
NO and mν > 0.1 eV for IO [2],Pand direct searches energy [18–39]; the Hubble constant H0 [40–43] and the
in the PKATRIN experiment imply mν < 1.35 eV [9]; amplitude parameter σ8 [40, 44, 45], that are in tension
where mν ≡ m1 + m2 + m3 is the total neutrino mass. among different datasets [46–51]; CMB lensing [26, 52–
Even if neutrino masses are tiny, determining them is a 57], where some observations are anomalous [40, 58–60];
mandatory first step to understand not only the origin or new physics in the neutrino sector [61–82].
of the particle masses and hierarchies, but also the Uni- To clarify the robustness of cosmological neutrino-mass
verse’s evolution and structure formation—neutrinos are, determinations, and the complementarity with direct
so far, the only known form of (hot) dark matter. searches, it is key to understand the physical effects of
In the early Universe, neutrinos are relativistic and neutrino masses that cosmology is most sensitive to [83].
behave as radiation, constituting up to 40% of the total Is current and upcoming data only sensitive to the back-
ground energy in neutrino masses, degenerate with other
sources of energy and cosmological unknowns? Or is it
∗ sensitive to a characteristic scale-dependent imprint, in-
[email protected]
† [email protected]
duced by neutrinos not moving at the speed of light?
‡ [email protected] In this paper, we explicitly separate background
§ [email protected] neutrino-mass effects, which capture the evolution of the
¶ [email protected] average neutrino energy density; from perturbations ef-
 2

fects, which capture “kinematic” scale-dependent effects tons, baryons, cold dark matter, neutrinos and a cosmo-
directly related to the neutrino speed. As a first step, logical constant). This is the only equation where neu-
we focus on their impact on CMB anisotropies, leaving trinos modify the background evolution of the Universe.
other observables for future work [84]. To gain intuition, The time dependence of the neutrino energy density
we study in detail the physical origin of detectable ef- ρν is, in turn, governed by the covariant conservation of
fects. Finally, we carry out an analysis of current CMB the neutrino stress-energy tensor
data, finding out that it mostly constraints background 1 dρν
neutrino-mass effects. In turn, the limit on perturbations = −3[1 + wν (a)] , (2)
ρν d ln a
effects is significantly weaker than the standard bounds.
Our results provide insight into the physical origin with wν ≡ Pν /ρν the neutrino equation of state and
of cosmological neutrino-mass bounds, and serve as a Pν the neutrino pressure. That is, the cosmological im-
benchmark of extended models that could affect them. pact of neutrinos at the background level is fully deter-
As one direct consequence, we find that any change on mined by their energy density at an initial time (usually
the background expansion history that is degenerate with parametrized in terms of Neff ) and their equation of state,
the impact of neutrino masses can strongly affect the cos- which controls how fast ρν dilutes.
mological neutrino-mass determination. Our definitions Differences among particle physics models enter when
of background and perturbations effects are not tied to specifying wν . For massive, non-interacting neutrinos
particular observables (e.g., CMB anisotropies are per- that decoupled while being relativistic,
turbations, but their evolution is sensitive to the back- 2
d p √ 2p 2 f0 (ap)
R 3
ground expansion rate), providing a generic framework to 1 p +mν
understand the cosmological effects of neutrino masses. wν (a, mν ) = R p , (3)
3 d3 p p2 + m2ν f0 (ap)
The structure of the manuscript is as follows. In Sec-
tion II, we describe our formalism to separate background with p linear momentum, mν each neutrino mass and
and perturbations effects. In Section III, we examine and 1 1
illustrate their physical effects on CMB anisotropies at f0 (ap) = 3 ap/T ν (4)
(2π) e 0 + 1
different angular scales. In Section IV, we present our
statistical analysis together with the cosmological con- a redshifted Fermi-Dirac distribution, with T0ν ≃ 1.9 K
straints on the key parameters of this study. Finally, we the current neutrino temperature. Equation (3) encloses
draw our conclusions in Section V. In the Appendices, the only background effect of neutrino masses.
we provide further explanations and results.
Background effect
II. FORMALISM

After decoupling from the primordial plasma, neutri- ( oscillations)

nos only affect cosmological observables via gravity, i.e.,
(CMB)
via their stress-energy tensor Tνµ that enters the Einstein
(KATRIN)
equations. In this Section, we set up the formalism to
separate background and perturbations effects. We start
from the most general parametrization, and compute the
relevant quantities for the specific case of massive, non-
interacting neutrinos. On top of establishing our frame-
work, this identifies the neutrino-mass effects that cos-
mology can be sensitive to, enabling a better understand-
ing of the cosmological neutrino-mass determination.

A. Evolution of the background Recombination

At the background level, the Einstein equations lead

to the first Friedmann equation
2
a′

8πG 2
P
FIG. 1. Impact of mν on the neutrino equation of state wν ,
= a (ργ + ρb + ρcdm + ρν + ρΛ ) , (1)
a 3 that controls how the neutrino energy density dilutes. When
the neutrino temperature drops below their mass, wν varies
with primes derivatives with respect to conformal time, from 1/3 (radiation) to 0 (matter). The only background effect
a the scale factor, G Newton’s constant, and ρi the back- of neutrino masses is to change how fast the neutrino energy
ground energy density of the species i (we include pho- dilutes, which modifies the expansion rate of the Universe.
 3

PFigure 1 shows how, by changing the equation of state, where δ ≡ δρ/ρν , ϕ and ψ are the gravitational potentials
mν determines how fast the neutrino energy density in the conformal Newtonian gauge, c2s ≡ δP/δρ is the so-
dilutes. We plot the equation of state of non-interacting called squared sound-speed, and
massive neutrinos for different total neutrino masses.
The chosen values are the smallest Pν′ wν′
P value allowed by neu- c2ad ≡ = w ν − ′ , (11)
trino oscillation experiments, mν > 0.06 eV [1, 2, 85]; ρ′ν a
3 a (1 + wν )
the current cosmological limit from Planck CMB data
is the so-called adiabatic squared sound-speed. Physi-
P
alone, mν < 0.24 eV [40]; and the current limit from
cally, the first three terms in Eq. (9) correspond to energy
P
the KATRIN experiment, m ν < 1.35 eV [9].
dilution due to bulk motions, gravitational redshift, and
P
As the figure shows, mν controls at which time wν
switches from 1/3 (radiation, ρν ∝ a−4 ) to 0 (matter, the expansion of the Universe; respectively. In Eq. (10),
ρν ∝ a−3 ). The change happens earlier for higher the first term corresponds to drag due to the expansion
P
mν ,
leading to a slower dilution of ρν in Eq. (1) due to the en- of the Universe; the second and third terms to isotropic
ergy in neutrino masses. Hence, the only background ef- and anisotropic momentum flow, respectively; and the
fect of increasing neutrino masses, with other parameters fourth term to gravitational forces.
fixed, is to increase the expansion rate of the Universe. To solve these equations, c2s and σ (together with the
The current CMB limit on background neutrino effects background quantity wν ) have to be provided. An in-
can thus be fully recast as enforcing that the neutrino determinacy arises because c2s is gauge-dependent, that
energy density dilutes as radiation approximately until is, its value depends on the coordinate system used to
recombination. While in the minimal scenario only neu- separate background from perturbations. Physically, by
trino masses control this effect, it is rather indirect, and performing a coordinate transformation one can transfer
non-minimal cosmological extensions —either a differ- pressure and energy density from perturbations to the
ent neutrino equation of state or additional background background and vice versa. A related issue is that, for
components— may mimic it [61, 62, 65–82]. adiabatic perturbations that behave as background as
Below, we separate this effect from perturbations ef- k → 0, the evolution equations for small k should only
fects. To do so, we parametrize depend on background quantities. This may enforce a
P wν as in the massive, non-
interacting case wν = wν (a, mBackg. consistency relation between c2s , σ, and wν . These issues
P Backg.
ν ); with mν
a parameter that we name background neutrino mass. can be overcome if, instead of c2s , the equations are writ-
ten in terms of the so-called effective sound-speed [87–89]
′
B. Evolution of perturbations k 2 c2s δ + 3 aa (1 + wν )c2ad θ
c2eff ≡ ′ . (12)
k 2 δ + 3 aa (1 + wν )θ
At the perturbations level, for scalar perturbations
the neutrino stress-energy tensor Tνµ can be written in As we show in Appendix B, c2eff is gauge-invariant and,
terms of the energy density perturbation δρ, the pres- when expressing the evolution equations in terms of it,
sure perturbation δP , the velocity divergence θ, and the adiabatic perturbations behave as background as k → 0
anisotropic stress σ. In Fourier space, [86] regardless of the values of c2eff , σ, and wν . Physically, ceff
is the sound speed in a frame comoving with neutrinos.
δρ ≡ −T00 − ρν , (5)
The equations above always hold. Differences among
1 particle physics models (including neutrino-mass effects)
δP ≡ Tii − Pν , (6)
3 enter when specifying c2eff and σ. For massive, non-
ik i Ti0 interacting neutrinos; they can be computed from the
θ≡ , (7) perturbed distribution function
ρν + P ν
k̂i k̂j [Tji − (Pν + δP )δji ] h i
σ≡− , (8) f (⃗k, p⃗, η) = f0 (ap) 1 + Ψ(⃗k, p⃗, η) , (13)
ρ ν + Pν
with k the comoving wavenumber. Physically, θ repre- with η conformal time. The explicit expression for the
sents the divergence of the bulk velocity of energy per- stress-energy tensor leads to [86]
turbations; and δP and σ the isotropic and anisotropic Z
components of the linear momentum flux, respectively. 1 p
δ= d3 p p2 + m2ν f0 Ψ , (14)
These parameters source the perturbed Einstein equa- ρν
Z
tions [86]. They evolve following energy-momentum con- 1
servation, which in the conformal Newtonian gauge reads θ= d3 p (i⃗k · p⃗)f0 Ψ , (15)
ρ ν + Pν
Z
a′ 2 1
δ ′ = −(1 + wν )θ + 3(1 + wν )ϕ′ − 3 δP = d3 p pvf0 Ψ , (16)


cs − wν δ , (9) 3
a
a′
 
c2s
Z
1 1
θ′ = − (1 − 3c2ad )θ + k2 δ − k2 σ + k2 ψ , (10) σ= d3 p pv − (k̂ · p̂)2 f0 Ψ , (17)
a 1 + wν ρ ν + Pν 3
 4

Perturbations effects
( oscillations) (CMB) (KATRIN)

kFS kFS kFS

kFS(a)>k

kFS(a)>k
kFS(a)<k

kFS(a)<k

kFS(a)<k
Radiation

Radiation

Radiation

mν on the comoving neutrino squared-sound-speed c2eff [left] and anisotropic stress σ [center, for fixed
P
FIG. 2. Impact of
scale and varying time; right, for fixed time and varying scale]. These parameters control all perturbations effects of neutrinos
(see text). c2eff is essentially scale-independent and falls when the neutrino temperature drops below their mass. In turn, σ gets
suppressed above a characteristic scale, the neutrino free-streaming length, that depends on the neutrino mass.

p
where v ≡ p/ p2 + m2ν is the neutrino velocity. These tions get steadily damped with time, as neutrino free-
expressions explicitly show the physical meaning of θ, δP , streaming steadily suppresses perturbations. However,
and σ for massive, non-interacting neutrinos. The evolu- when the mode becomes larger than the free-streaming
tion of Ψ follows the perturbed Boltzmann equation [86]. scale (which shrinks with time as neutrinos become non-
Figure 2 shows that the anisotropic
P stress contains the relativistic), i.e., when k < kFS ; neutrinos cluster instead
leading “kinematic” effects of mν . We plot the time- of free-streaming, the momentum flux diminishes, and σ
and scale-dependence of c2eff and σ (normalized to an ini- decays much faster. The right panel also shows the two
tial comoving curvature perturbation R = 1) of non- distinct behaviors as a function of scale at fixed time. At
interacting massive neutrinos, for the same total neutrino scales below the free-streaming scale, k > kFS , σ is larger
masses as in Fig. 1. We fix other cosmological parameters at large scales, which had less time to evolve and are less
to the best fit of the Planck 2018 CMB analysis [40]. damped by free-streaming. On the contrary, at scales
The scale dependence can be understood in terms of above the free-streaming scale, k < kFS , σ is smaller at
the neutrino free-streaming wavenumber [83] large scales, where neutrino clustering reduces the mo-
mentum flux. These distinct behaviors are a direct con-
a2 H(a)  mν  sequence of neutrinos not moving at the speed of light.
kFS (a) ≃ 0.776 h Mpc−1 , (18)
H0 1 eV To separate the background and perturbations effects
with H the Hubble parameter and h the reduced Hub- of neutrino masses, below we parametrize c2eff and σ as
ble constant. Physically, if k < kFS , perturbations tend in the massive,
P non-interacting case; computing them for
to collapse gravitationally; whereas if k > kFS , velocity a valuePof mν that we name perturbations neutrino
dispersion inhibits gravitational collapse [83].PThus, kFS mass, mPert.
ν (see Appendix A for details).
directly encodes the “kinematic” impact of mν ̸= 0,
i.e., that neutrinos do not move at the speed of light.
As the left panel of Fig. 2 shows, the time dependence
of c2eff resembles that of the equation of state, falling from III. PHENOMENOLOGY
1/3 to 0 when neutrinos become non-relativistic. The
scale dependence introduced by kFS is subleading (see
Ref. [90] for a discussion). Physically, super-horizon adi- The formalism described above allows to implement
abatic perturbations behave as background, c2eff → c2ad at neutrino background effects, parametrized by wν ; and
all scales, and scale dependence only appears after sub- perturbations effects, parametrized by c2eff and σ; in a
leading sub-horizon evolution [90, 91]. model-independent way. In this Section, we use this for-
The center and right panels of Fig. 2 show that malism to separate the background and perturbations ef-
the evolution of σ is much more scale-dependent, with fects of neutrino masses. We discuss the physical effects
a characteristic feature at k = kFS . As the center on CMB anisotropies, their origin, and how they can be
panel shows, at early times σ oscillates. These oscilla- observationally identified.
 5

A. Implementation with cs,γ (z) the sound speed of the baryon-photon fluid
and zrec the redshift of recombination. Since θs is very-
To implement background effects, we solve Eq. (2) for well measured, changes in it are compensated by mod-
an equation of state ifying H0 , which modifies the cosmological constant Λ.
P ofBackg.
massive, non-interacting neutrinos
This changes the Λ-induced late-time boosting of large-
with total mass mν ,
 Z a scale anisotropies through the late integrated Sachs-
1 + wν (a, mνBackg. )
P 
ρν (a) = ρν (a1 ) exp −3 da , Wolfe (LISW) effect [95]. In addition, the amplitude of
a1 a super-horizon CMB perturbations depends on the expan-
(19) sion rate of the Universe aroundP recombination (see Ap-
where a1 is an initial scale factor. For simplicity, we pendix C), which depends on mBackg. ν . This introduces
assume the standard initial energy density for three light a subleading depletion of low-ℓ CMB anisotropies.
neutrino species, [83] Secondly, high multipoles (ℓ ≳ 500) correspond to
7π 2 ν 4 modes that enter the horizon much before recombina-
ρν (a1 ) = (T ) , (20) tion. There are two main ways in which neutrino masses
40a41 0 affect these modes, both of which are mainly sensitive to
with the initial condition evaluated P when neutrinos are background effects.
ultrarelativistic, i.e., T0ν /a1 ≫ mBackg.
ν . The hereby On the one hand, these modes are damped below a
computed neutrino energy density affects the expansion characteristic angular scale θD due to the finite mean
of the Universe via the Friedmann equation, Eq. (1). free path of photons, where [83]
To implement perturbations effects, we solve Eqs. (9) sZ ,Z
and (10) for values ofPc2eff and σ corresponding to neutri- ∞
1 + z dz zrec
dz
nos with total mass mPert. , together with the pertur- θD ∼ , (22)
ν
zrec n e (z)σT H(z) 0 H(z)
bations equations of other species. We set initial condi-
tions corresponding to adiabatic perturbations of an ul- with ne the electronPdensity and σT the Thomson scatter-
trarelativistic relic [86]. We provide the technical details ing cross section. mBackg. slows down the dilution of
ν
of our computation in Appendix A. The hereby computed the neutrino energy density, increasing H(z) both in the
neutrino energy density and pressure perturbations δρν numerator and denominator. Overall, θD gets reduced,
and δPν , velocity divergence θ, and anisotropic stress σ which is visible as an excess at high ℓ in Fig. 3. As the
affect the evolution of the Universe via the perturbed
Einstein equations. We implement our modified cosmo-
logical evolution in the public code CLASS [91–94].

B. Physical effects on the CMB

Mass only in background

mBackg. mPert.
P P
Figure 3 shows that ν and ν leave dis-
tinct signatures on CMB anisotropies at different scales.
To better represent the observable effects, in the fig-
ure we fix the well-measured cosmological parameters
Mass both in background
{100θs , ωb , ωcdm , As , ns , τreio } to the best fit of the and perturbations
Planck 2018 CMB analysis [40]. Below, we explore in
detail the physical origin of the different effects. In Ap- Mass only in perturbations
pendix C, we provide explicit checks and figures that fur-
ther illustrate our results.
Firstly, low multipoles (ℓ ≲ 10) correspond to modes
that enter the horizon at very late times, when neutrinos
constitute a subleading component of the energy den-
sity of the Universe. The main neutrino-mass effects are
thus induced when the modes are larger than the hori-
zon. Since, for adiabatic perturbations, the evolution of
super-horizon modes Pdepends only on background quan-
tities, the effect of PmPert.ν at these scales is negligible. FIG.
P 3. Background- and perturbations-induced impact of
The main effect of m Backg.
is indirect: an increased mν on CMB anisotropies. As detailed in the main text;
ν
P Backg.
mν slows down the dilution of the neutrino energy at low and high ℓ, the background affects the LISW effect,
Silk damping, and lensing; at intermediate ℓ, both back-
density, increasing H(z) and, in principle, modifying the
ground and perturbations directly couple to photon-baryon
observed angular scale of CMB peaks [83] oscillations via gravity. Most effects are background-induced.
Z ∞  Z zrec
cs,γ (z) dz dz Perturbations-induced effects are mainly relevant at interme-
θs = , (21) diate ℓ, where their effect is opposite to that of the background.
zrec H(z) 0 H(z)
 6
P Backg.
damping scale only depends P onPert.
background quantities, logical evolution. In short, mν encodes the equa-
this effect is insensitive to mν (see Appendix C). tion of state,P
i.e., how fast the neutrino energy density di-
On the other hand, high multipoles are affected by lutes; while mPert.
ν contains a more direct “kinematic”
weak gravitational lensing, i.e., by the random gravita- effect related to the free-streaming nature of neutrinos.
tional deflection of CMB photons due to the large-scale In this Section, we carry out anPanalysis of CMB P data to
structure of the Universe. Lensing smooths out the power quantify the allowed values of mBackg.
ν and m Pert.
ν .
spectrum and transfers power from low multipoles to We analyze the Planck 2018 temperature, polarization,
high multipoles [96]. Neutrino masses affect CMB lens- and lensing power spectra (TT, TE, EE+lowE+lensing
ing via both background P and perturbations effects (see in Ref. [40]). To do so, we modify the public code
Appendix C for details). mBackg.
ν reduces CMB lens- CLASS [91–94], as mentioned above and detailed in Ap-
ing by accelerating the expansion of the Universe, which pendix A, to solve the evolution of cosmological pertur-
suppresses structure formation. This is visible in Fig. 3 bations; and we explore the parameter space with the
as wiggles thatP are in phase with the CMB power spec- public Markov Chain Monte Carlo code COBAYA [97, 98].
trum. In turn, mPert. ν enhances neutrino clustering and In Appendix D, we summarize our priors and convergence
structure growth, which enhances CMB lensing in a scale-
P Backg. criteria, and we provide the full results of our analysis.
dependent way. This effect is relevant if mν ̸= 0, Fig. 4 shows that splitting neutrino mass effects among
because otherwise the neutrino energy density is too di- background and perturbations strongly increases the al-
luted at late times and neutrino perturbations
P Pert.effects are lowed perturbations effects. We show in solid the 1D pos-
negligible. Overall, the main impact of mP ν on CMB terior
P probabilities
P and 2D credible regions of our analysis
lensing is to partially reduce the effect of mBackg.
ν as for mBackg.
ν , mPert.
ν , the Hubble parameter H0 , and
can be seen in Fig. 3. In our data analysis below, this the amplitude parameter σ8P . Dashed lines
Pcorrespond to
leads to a partial degeneracy among both parameters. the standard scenario, i.e., mBackg.
ν = m Pert.
ν .
Finally, intermediate multipoles (10 ≲ ℓ ≲ 500) cor- The 95% CL limits on the total neutrino mass read
respond to modes that enter the horizon around recom- X
bination. Anisotropies at these scales are largely influ- mν < 0.24 eV (Planck 2018 [40]) ,
enced by the gravitational potentials around recombina- X
tion, that act as a driving term for photon-baryon acous- mBackg.
ν < 0.29 eV (This paper) ,
tic oscillations. More precisely, decaying gravitational
X
Pert.
mν < 0.79 eV (This paper) .
potentials
P Backg. increase the amplitude of oscillations [89].
mν boosts this decay by increasing the expan- That is, the standard limit on
P
mν is mostly a limit on
sion rate of the Universe, thus boosting the amplitude of P Backg.
mν (the limit on
anisotropies as can be seen in Fig. 3. In turn,
P Pert.
mν PthePert.
latter is slightly weaker due
to a degeneracy with mν that we discuss below).
enhances neutrino clustering and structure growth in a As the figure shows, the posterior
P
scale-dependent way (see Appendix C), slowing down the P probabilities of mν ,
in the standard analysis; and mBackg.
ν , in our analysis;
decay of the gravitational potentials and decreasing the almost
P
amplitude of anisotropies. On top of these effects, decay- P Pert.match. The correspondence among mν and
mν , however, is null.
ing gravitational potentials further increase the ampli-
In other words, CMB data tightly constraints the neu-
tude of CMB anisotropies through the integrated Sachs-
trino equation of state; but, compared to the standard
Wolfe effect, boosting the aforementioned P effects (see Ap-
scenario, the limit on “kinematic” effects of neutrino
pendix C). The fact that the effects of mBackg.
ν and
P Pert. masses is relaxed by about a factor of 3.
mν are opposite at intermediate ℓ leads to a partial
The background neutrino mass is correlated with H0
degeneracy among both parameters in our data analysis.
and σ8 . These correlations are also present in the stan-
We conclude that,P while background effects are present dard scenario, and they are due to the neutrino contribu-
at all multipoles, mPert.
ν on its own is mainly relevant tion to the total energy density of the Universe. Larger ρν
at intermediate ℓ, where anisotropies are suppressed due increases the expansion rate, which suppresses structure
to the direct gravitational impact of neutrino perturba- formation (i.e., σ8 ); and, as described in the previous Sec-
tions. EvenP ifBackg.
this effect can be partially hidden by a tion, modifies the angular scale of the CMB peaks that is
non-zero mν , unique background effects at high ℓ degenerate
allow disentangling both. In the next section, we carry P Pert.with H0 . These are both background effects,
so mν is not strongly correlated with H0 and σ8 .
out a data analysis to quantify the effects allowed by cur- Figure 5 shows that intermediate-ℓ data constrains
rent CMB data. The higher statistics at large ℓ forecasts both perturbations and background mass effects, whereas
a stronger limit on background effects. high-ℓ data mainly constrains background-mass effects.
In both cases, perturbations-mass effects can partially
compensate background mass effects; we discuss this be-
IV. DATA ANALYSIS AND RESULTS low. As in Fig. 3, we fix {100θs , ωb , ωcdm , As , ns , τreio }
to their Planck 2018 best-fit values [40]. The blue and
In the Sections above, we have discussed how the differ- purple lines are close to the 2σ allowed region in our
ent neutrino masses that we introduce affect the cosmo- analysis, whereas the red line is more strongly excluded.
 7

Mass in background

Mass in perturbations

Mass in background and perturbations

FIG. 5. Impact on CMB anisotropies of parameters ex-

FIG. 4. CMB limits on the separate neutrino-massP effects due cluded by our analysis. Background effects are most excluded
Pert.
to background and perturbations. The P limit on m ν P gets by high-ℓ data (CMB lensing), and perturbations effects by
strongly relaxed compared to that on mν . The CMB mν intermediate-ℓ data (direct coupling via gravity). Opposed
limit is mostly a limit on background neutrino-mass effects. physical effects induce a degeneracy between
P Pert.
m ν and
P Backg.
mν (see text).

As discussed in the previous Section, the sensitiv- P

ity at intermediate multipoles is due to the direct im- bound on mν is mainly a limit on background effects,
pact of neutrino masses on the gravitational potentials. i.e., on how the neutrino energy density evolves with
The high-ℓ sensitivity to background mass effects can time. Once we separate neutrino-mass effects into back-
be traced back to the impact on the damping tail and ground and perturbations (the latter containing more
CMB lensing. This is particularly relevant in light of the direct “kinematic” signatures related to neutrino free-
lensing anomaly: Planck 2018 data prefers more CMB streaming), the limit on perturbations effects gets largely
lensing than what is present in the standard ΛCDM relaxed. This could serve as an insight for theories that
scenario [40, 58, 99–101] (see, however, Refs.P[59, 60]). evade cosmological neutrino mass bounds, as our results
Since, as described in the previous Section, mBackg.
ν show that such theories should have a radiation-like dilu-
reduces CMB lensing, the anomaly enhances the high-ℓ tion of the neutrino energy density until recombination.
constraints, and the resulting limit is somewhat stronger
than what would be expected. This is visible in the bot-
tom panel of Fig. 5 as wiggles in the data that are out of
P Backg. V. CONCLUSIONS AND OUTLOOK
phase with the effect of mν .
An interesting feature of Fig. 4, also visible in Fig. 5, The absolute neutrino mass scale remains unknown.
is that when splittingP neutrino masses P into background Cosmological surveys are close to either measuring it or
and perturbations, PmBackg.ν and mPert.
ν are partly excluding the lower limit set by oscillations. In this pa-
Pert.
degenerate.
P Backg. When m ν increases, the limit over per, we explore the origin of the cosmological neutrino-
mν can relax by almost a factor of two. This mass bound from CMB data. Our results show that the
is driven Pby opposite effects of background and pertur- driving constraint arises from the contribution of the en-
bations: mνBackg. suppresses structure formation due ergy in neutrino masses to the expansion of the Universe.
to
P the increased expansion rate of the Universe, whereas We focus on the macroscopic quantities that capture
mPert.
ν enhances structure formation due to neutrinos the cosmological impact of neutrinos. This enables us to
clustering as matter. As described in the previous Sec- explicitly discriminate, for the first time, between back-
tion, this introduces opposite effects on the amplitude of ground and perturbations effects of neutrino masses. The
CMB anisotropies, at intermediate ℓ; and on CMB lens- separation can be understood in terms of standard “fluid”
ing, at high ℓ. TheP degeneracy is not perfect, partly be- variables: the equation of state wν , which governs how
cause increasing mBackg.
ν also
P affects the CMB damp- fast the background neutrino energy density dilutes; the
ing tail that is insensitive to mPert.
ν . sound speed in the frame comoving with neutrinos c2eff ,
Overall, our results explicitly show that the CMB which captures isotropic neutrino momentum flow; and
 8

the anisotropic stress σ, which captures anisotropic neu- fects of neutrino masses. Since current and near-future
trino momentum flow. wν encloses all background effects, observations of the matter power spectrum (particularly
directly impacting the P expansion rate of the Universe— Baryon Acoustic Oscillation, BAO, measurements) have
which increases as mν increases due to the energy in a strong impact on neutrino-mass determinations [10, 11],
neutrino masses. In turn, c2eff and σ enclose the perturba- we will explore them in detail in incoming work [84]. For
tions effects. c2eff is almost scale-independent, and σ con- BAOs, which measure quantities that can be expressed in
tains the main “kinematic” impact of neutrino masses, terms of background neutrino
P Backg. properties P
[61], we foresee
related to the free-streaming scale set by neutrinos not a stronger impact on mν than on mPert.
ν .
moving at the speed of light. Our separation among background and perturbation
Since our goal is to disentangle among background and neutrino-mass effects opens many research avenues. The
perturbations effects, we explore two types of neutrino observed
P Backg.CMB-lensing excess drives the P strong limit on
masses: the one that P governs wν , that we name back- mν and the best-fit for nonzero mPert.ν in Fig. 4.
ground neutrino mass mBackg.
ν ; and the one that gov- In light of this, our framework could be explored with
2 state-of-the-art CMB likelihoods where this anomaly is
erns c
P Pert.eff and σ, that we name perturbations neutrino mass P Backg.
mν . Although these parameters are phenomeno- not present
P Pert.[56, 59, 60]. The opposite effects of m ν
logical, they encode the distinct physical implications of and mν on CMB lensing, together with the scale-
neutrino masses in cosmology. Hence, they serve as a dependence of perturbations effects, could also be lever-
benchmark to understand the effects that cosmology is aged to separate both effects in future high-precision de-
most sensitive to, and to shed light P on potential degen- terminations of CMB lensing [16, 17, 101–103]. More-
eracies of extended models with mν . over, some cosmological tensions are correlated with neu-
The effects on the CMB temperature anisotropies trino masses, motivating new studies that address them
can be split into different multipole regions (see Sec- in the context of our separation of neutrino-mass effects.
tion III B). The low-ℓ region has a minor impact This may shed light on the physics that can alleviate
P on the
constraints due to cosmic variance. At high ℓ, mBackg.
ν
these tensions. Beyond pure cosmology studies, our re-
modifies the CMB damping tail, and both background sults provide a benchmark to build models that evade the
and perturbations effects impact P CMB lensing in op- cosmological neutrino-mass bound. In short, the main
posite directions. While a larger mBackg.
ν suppresses way to relax the cosmological bound is by modifying the
structure formation (and hence CMB lensing) by in- expansion history of the Universe, with the “kinematic”
creasing properties of the model being less important.
P Pert. the expansion rate of the Universe, a larger
mν enhances structure formation above the neu- While we wait for a positive neutrino-mass signal from
trino free-streaming scale. At intermediate ℓ, neutrino terrestrial experiments, cosmological measurements lead
perturbations directly affect photon-baryon current limits. As they improve, we find ourselves in an
P oscillations
via gravity,
P Pert. where again the effects of mBackg.
ν and era where cosmological limits are approaching values dis-
mν are opposite. P As a consequence, there is a favored by oscillation experiments. If this tension grows,
slight degeneracy among mBackg.
P Pert.
ν and m ν . Over- the solution may rely on a non-standard neutrino sector.
all, high-ℓ data is mostly sensitive to background effects, If, in turn, a neutrino-mass signal is found, scrutinizing
whereas intermediate-ℓ data determines both background the robustness of this determination will be mandatory.
and perturbations effects. In both cases, understanding the involved physical effects
We then carry out an analysis of Planck 2018 CMB and degeneracies is key for solid progress in cosmology
data to shed light on the effects that observations con- and particle physics. Future observations will guide the
strain. We conclude that the Planck 2018 neutrino-mass next steps for the physics of neutrinos, the first particle
bound is a bound on the background effects, i.e., on the whose mass may be first measured outside laboratories.
evolution of the neutrino energy density. This provides
a rule-of-thumb to understand if CMB data excludes a
model with new physics in the neutrino sector: if its equa- ACKNOWLEDGMENTS
tion of state significantly deviates from wν = 1/3 around
recombination, the model is probably excluded. This work has been supported by the Span-
The perturbations limit on “kinematic” effects of neu-
P Pert. ish MCIN/AEI/10.13039/501100011033 grants PID2020-
trino masses is consequently relaxed, mν < 0.8 eV. 113644GB-I00, PID2023-148162NB-C22 (RH and OM),
The limit is still competitive—models that dramati- PID2021-123703NB-C21, PID2022-136510NB-C33 (IE),
cally affect free-streaming properties of neutrinos are still and PID2022-126224NB-C21 (TB and JS); and by
excluded—, and it is similar to the projected reach of KA- the European Union’s Horizon 2020 research and
TRIN [9]; yet in the standard scenario such high neutrino innovation program under the Marie Sklodowska-
masses are excluded within ∼ 7σ. This result underscores Curie grants HORIZON-MSCA-2021-SE-01/101086085-
the complementarity among laboratory and cosmological ASYMMETRY and H2020-MSCA-ITN-2019/860881-
determinations of the neutrino mass. HIDDeN. RH is supported by the Spanish grant FPU19/
In this first paper, we have focused on the consequences 03348 of MU, and TB by the Spanish grant PRE2020-
for the CMB of “energy-dilution” versus “kinematic” ef- 091896. The authors also acknowledge support from the
 9

Generalitat Valenciana grants PROMETEO/2019/083 the Institut de Ciències del Cosmos of the Universitat
and CIPROM/2022/69 (RH and OM). IE acknowledges de Barcelona (RH), and Institut de Fı́sica Corpuscular
support from the Basque Government (IT1628-22). The (IE and TB) for their hospitality during the comple-
authors acknowledge the Galileo Galilei Institute (GGI) tion of this work. OM acknowledges the financial sup-
for Theoretical Physics (RH and OM), the Center for port from the MCIU with funding from the European
Cosmology and AstroParticle Physics (CCAPP) at the Union NextGenerationEU (PRTR-C17.I01) and Gener-
Ohio State University (RH), the Fermi National Accel- alitat Valenciana (ASFAE/2022/020). TB and JS ac-
erator Laboratory FERMILAB (TB, RH, OM and JS), knowledge support from the “Unit of Excellence Maria
the University of the Basque Country UPV/EHU (RH), de Maeztu 2020-2023” award to the ICC-UB CEX2019-
000918-M.

[1] P. F. de Salas, D. V. Forero, S. Gariazzo, [14] Euclid Collaboration, M. Archidiacono et al., Euclid
P. Martı́nez-Miravé, O. Mena, C. A. Ternes, preparation. Sensitivity to neutrino parameters,
M. Tórtola, and J. W. F. Valle, 2020 global 2405.06047.
reassessment of the neutrino oscillation picture, JHEP [15] DESI Collaboration, A. Aghamousa et al., The
02 (2021) 071, [2006.11237]. DESI Experiment Part I: Science,Targeting, and
[2] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, Survey Design, 1611.00036.
I. Martinez-Soler, J. a. P. Pinheiro, and T. Schwetz, [16] Simons Observatory Collaboration, P. Ade et al.,
NuFit-6.0: Updated global analysis of three-flavor The Simons Observatory: Science goals and forecasts,
neutrino oscillations, 2410.05380. JCAP 02 (2019) 056, [1808.07445].
[3] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, [17] LiteBIRD Collaboration, E. Allys et al., Probing
A. Melchiorri, and A. Palazzo, Unfinished fabric of the Cosmic Inflation with the LiteBIRD Cosmic
three neutrino paradigm, Phys. Rev. D 104 (2021), Microwave Background Polarization Survey, PTEP
no. 8 083031, [2107.00532]. 2023 (2023), no. 4 042F01, [2202.02773].
[4] T2K Collaboration, K. Abe et al., Measurements of [18] S. Hannestad, Neutrino masses and the dark energy
neutrino oscillation parameters from the T2K equation of state - Relaxing the cosmological neutrino
experiment using 3.6 × 1021 protons on target, Eur. mass bound, Phys. Rev. Lett. 95 (2005) 221301,
Phys. J. C 83 (2023), no. 9 782, [2303.03222]. [astro-ph/0505551].
[5] NOvA Collaboration, M. A. Acero et al., Improved [19] W. Yang, E. Di Valentino, S. Pan, and O. Mena,
measurement of neutrino oscillation parameters by the Emergent Dark Energy, neutrinos and cosmological
NOvA experiment, Phys. Rev. D 106 (2022), no. 3 tensions, Phys. Dark Univ. 31 (2021) 100762,
032004, [2108.08219]. [2007.02927].
[6] JUNO Collaboration, F. An et al., Neutrino [20] S. Vagnozzi, S. Dhawan, M. Gerbino, K. Freese,
Physics with JUNO, J. Phys. G 43 (2016), no. 3 A. Goobar, and O. Mena, Constraints on the sum of
030401, [1507.05613]. the neutrino masses in dynamical dark energy models
[7] Hyper-Kamiokande Collaboration, K. Abe et al., with w(z) ≥ −1 are tighter than those obtained in
Hyper-Kamiokande Design Report, 1805.04163. ΛCDM, Phys. Rev. D 98 (2018), no. 8 083501,
[8] DUNE Collaboration, B. Abi et al., Deep [1801.08553].
Underground Neutrino Experiment (DUNE), Far [21] S. Roy Choudhury and A. Naskar, Strong Bounds on
Detector Technical Design Report, Volume II: DUNE Sum of Neutrino Masses in a 12 Parameter Extended
Physics, 2002.03005. Scenario with Non-Phantom Dynamical Dark Energy
[9] Katrin Collaboration, M. Aker et al., Direct (w(z) ≥ −1) with CPL parameterization, Eur. Phys. J.
neutrino-mass measurement based on 259 days of C 79 (2019), no. 3 262, [1807.02860].
KATRIN data, 2406.13516. [22] E. di Valentino, S. Gariazzo, and O. Mena, Model
[10] DESI Collaboration, A. G. Adame et al., DESI 2024 marginalized constraints on neutrino properties from
VI: Cosmological Constraints from the Measurements cosmology, Phys. Rev. D 106 (2022), no. 4 043540,
of Baryon Acoustic Oscillations, 2404.03002. [2207.05167].
[11] DESI Collaboration, A. G. Adame et al., DESI 2024 [23] M.-M. Zhao, Y.-H. Li, J.-F. Zhang, and X. Zhang,
VII: Cosmological Constraints from the Full-Shape Constraining neutrino mass and extra relativistic
Modeling of Clustering Measurements, 2411.12022. degrees of freedom in dynamical dark energy models
[12] J.-Q. Jiang, W. Giarè, S. Gariazzo, M. G. Dainotti, using Planck 2015 data in combination with
E. Di Valentino, O. Mena, D. Pedrotti, S. S. da Costa, low-redshift cosmological probes: basic extensions to
and S. Vagnozzi, Neutrino cosmology after DESI: ΛCDM cosmology, Mon. Not. Roy. Astron. Soc. 469
tightest mass upper limits, preference for the normal (2017), no. 2 1713–1724, [1608.01219].
ordering, and tension with terrestrial observations, [24] X. Zhang, Impacts of dark energy on weighing
2407.18047. neutrinos after Planck 2015, Phys. Rev. D 93 (2016),
[13] D. Wang, O. Mena, E. Di Valentino, and S. Gariazzo, no. 8 083011, [1511.02651].
Updating neutrino mass constraints with Background [25] R.-Y. Guo, J.-F. Zhang, and X. Zhang, Exploring
measurements, 2405.03368. neutrino mass and mass hierarchy in the scenario of
 10

vacuum energy interacting with cold dark matte, Chin. cosmological data, Phys. Rev. D 87 (2013), no. 4
Phys. C 42 (2018), no. 9 095103, [1803.06910]. 043515, [1211.2154].
[26] S. Roy Choudhury and S. Hannestad, Updated results [42] J.-P. Hu and F.-Y. Wang, Hubble Tension: The
on neutrino mass and mass hierarchy from cosmology Evidence of New Physics, Universe 9 (2023), no. 2 94,
with Planck 2018 likelihoods, JCAP 07 (2020) 037, [2302.05709].
[1907.12598]. [43] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez,
[27] H. Li and X. Zhang, Constraining dynamical dark S. J. Witte, V. Poulin, and J. Lesgourgues, The H0
energy with a divergence-free parametrization in the Olympics: A fair ranking of proposed models, Phys.
presence of spatial curvature and massive neutrinos, Rept. 984 (2022) 1–55, [2107.10291].
Phys. Lett. B 713 (2012) 160–164, [1202.4071]. [44] I. G. Mccarthy, S. Bird, J. Schaye, J. Harnois-Deraps,
[28] Y.-H. Li, S. Wang, X.-D. Li, and X. Zhang, A. S. Font, and L. Van Waerbeke, The BAHAMAS
Holographic dark energy in a Universe with spatial project: the CMB–large-scale structure tension and the
curvature and massive neutrinos: a full Markov Chain roles of massive neutrinos and galaxy formation, Mon.
Monte Carlo exploration, JCAP 02 (2013) 033, Not. Roy. Astron. Soc. 476 (2018), no. 3 2999–3030,
[1207.6679]. [1712.02411].
[29] J.-F. Zhang, Y.-H. Li, and X. Zhang, Cosmological [45] K. N. Abazajian and A. Kusenko, Hidden treasures:
constraints on neutrinos after BICEP2, Eur. Phys. J. Sterile neutrinos as dark matter with miraculous
C 74 (2014) 2954, [1404.3598]. abundance, structure formation for different production
[30] J.-F. Zhang, M.-M. Zhao, Y.-H. Li, and X. Zhang, mechanisms, and a solution to the σ8 problem, Phys.
Neutrinos in the holographic dark energy model: Rev. D 100 (2019), no. 10 103513, [1907.11696].
constraints from latest measurements of expansion [46] E. Abdalla et al., Cosmology intertwined: A review of
history and growth of structure, JCAP 04 (2015) 038, the particle physics, astrophysics, and cosmology
[1502.04028]. associated with the cosmological tensions and
[31] C.-Q. Geng, C.-C. Lee, R. Myrzakulov, M. Sami, and anomalies, JHEAp 34 (2022) 49–211, [2203.06142].
E. N. Saridakis, Observational constraints on varying [47] E. Di Valentino, O. Mena, S. Pan, L. Visinelli,
neutrino-mass cosmology, JCAP 01 (2016) 049, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and
[1504.08141]. J. Silk, In the realm of the Hubble tension—a review of
[32] Y. Chen and L. Xu, Galaxy clustering, CMB and solutions, Class. Quant. Grav. 38 (2021), no. 15
supernova data constraints on ϕ CDM model with 153001, [2103.01183].
massive neutrinos, Phys. Lett. B 752 (2016) 66–75, [48] L. Verde, T. Treu, and A. G. Riess, Tensions between
[1507.02008]. the Early and the Late Universe, Nature Astron. 3
[33] A. Loureiro et al., On The Upper Bound of Neutrino (2019) 891, [1907.10625].
Masses from Combined Cosmological Observations and [49] L. Verde, N. Schöneberg, and H. Gil-Marı́n, A tale of
Particle Physics Experiments, Phys. Rev. Lett. 123 many H0 , 2311.13305.
(2019), no. 8 081301, [1811.02578]. [50] E. Di Valentino et al., Cosmology Intertwined III: f σ8
[34] S. Wang, Y.-F. Wang, D.-M. Xia, and X. Zhang, and S8 , Astropart. Phys. 131 (2021) 102604,
Impacts of dark energy on weighing neutrinos: mass [2008.11285].
hierarchies considered, Phys. Rev. D 94 (2016), no. 8 [51] S. Roy Choudhury and S. Choubey, Updated Bounds
083519, [1608.00672]. on Sum of Neutrino Masses in Various Cosmological
[35] W. Yang, R. C. Nunes, S. Pan, and D. F. Mota, Scenarios, JCAP 09 (2018) 017, [1806.10832].
Effects of neutrino mass hierarchies on dynamical dark [52] E. Di Valentino and A. Melchiorri, Neutrino Mass
energy models, Phys. Rev. D 95 (2017), no. 10 103522, Bounds in the Era of Tension Cosmology, Astrophys.
[1703.02556]. J. Lett. 931 (2022), no. 2 L18, [2112.02993].
[36] Q.-G. Huang, K. Wang, and S. Wang, Constraints on [53] W. Giarè, O. Mena, and E. Di Valentino, Lensing
the neutrino mass and mass hierarchy from impact on cosmic relics and tensions, Phys. Rev. D
cosmological observations, Eur. Phys. J. C 76 (2016), 108 (2023), no. 10 103539, [2307.14204].
no. 9 489, [1512.05899]. [54] N. Craig, D. Green, J. Meyers, and S. Rajendran, No
[37] R. K. Sharma, K. L. Pandey, and S. Das, Implications νs is Good News, 2405.00836.
of an Extended Dark Energy Model with Massive [55] D. Green and J. Meyers, The Cosmological Preference
Neutrinos, Astrophys. J. 934 (2022), no. 2 113, for Negative Neutrino Mass, 2407.07878.
[2202.01749]. [56] D. Naredo-Tuero, M. Escudero,
[38] M. Zhang, J.-F. Zhang, and X. Zhang, Impacts of dark E. Fernández-Martı́nez, X. Marcano, and V. Poulin,
energy on constraining neutrino mass after Planck Living at the Edge: A Critical Look at the
2018, Commun. Theor. Phys. 72 (2020), no. 12 Cosmological Neutrino Mass Bound, 2407.13831.
125402, [2005.04647]. [57] I. J. Allali and A. Notari, Neutrino mass bounds from
[39] A. R. Khalifeh and R. Jimenez, Distinguishing Dark DESI 2024 are relaxed by Planck PR4 and
Energy models with neutrino oscillations, Phys. Dark cosmological supernovae, 2406.14554.
Univ. 34 (2021) 100897, [2105.07973]. [58] P. Motloch and W. Hu, Tensions between direct
[40] Planck Collaboration, N. Aghanim et al., Planck measurements of the lens power spectrum from Planck
2018 results. VI. Cosmological parameters, Astron. data, Phys. Rev. D 97 (2018), no. 10 103536,
Astrophys. 641 (2020) A6, [1807.06209]. [Erratum: [1803.11526].
Astron.Astrophys. 652, C4 (2021)]. [59] M. Tristram et al., Cosmological parameters derived
[41] E. Giusarma, R. De Putter, and O. Mena, Testing from the final Planck data release (PR4), Astron.
standard and nonstandard neutrino physics with Astrophys. 682 (2024) A37, [2309.10034].
 11

[60] E. Rosenberg, S. Gratton, and G. Efstathiou, CMB large-neutrino-mass cosmologies, Phys. Rev. D 105
power spectra and cosmological parameters from Planck (2022), no. 6 063501, [2111.14870].
PR4 with CamSpec, Mon. Not. Roy. Astron. Soc. 517 [78] I. M. Oldengott, G. Barenboim, S. Kahlen, J. Salvado,
(2022), no. 3 4620–4636, [2205.10869]. and D. J. Schwarz, How to relax the cosmological
[61] I. Esteban and J. Salvado, Long Range Interactions in neutrino mass bound, JCAP 04 (2019) 049,
Cosmology: Implications for Neutrinos, JCAP 05 [1901.04352].
(2021) 036, [2101.05804]. [79] J. Alvey, M. Escudero, and N. Sabti, What can CMB
[62] I. Esteban, O. Mena, and J. Salvado, Nonstandard observations tell us about the neutrino distribution
neutrino cosmology dilutes the lensing anomaly, Phys. function?, JCAP 02 (2022), no. 02 037, [2111.12726].
Rev. D 106 (2022), no. 8 083516, [2202.04656]. [80] A. Cuoco, J. Lesgourgues, G. Mangano, and S. Pastor,
[63] C. D. Kreisch, F.-Y. Cyr-Racine, and O. Doré, Do observations prove that cosmological neutrinos are
Neutrino puzzle: Anomalies, interactions, and thermally distributed?, Phys. Rev. D 71 (2005) 123501,
cosmological tensions, Phys. Rev. D 101 (2020), no. 12 [astro-ph/0502465].
123505, [1902.00534]. [81] I. J. Allali, D. Aloni, and N. Schöneberg, Cosmological
[64] M. Lattanzi and M. Gerbino, Status of neutrino probes of Dark Radiation from Neutrino Mixing, JCAP
properties and future prospects - Cosmological and 09 (2024) 019, [2404.16822].
astrophysical constraints, Front. in Phys. 5 (2018) 70, [82] C. Benso, T. Schwetz, and D. Vatsyayan, Large
[1712.07109]. neutrino mass in cosmology and keV sterile neutrino
[65] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, dark matter from a dark sector, 2410.23926.
Cosmological Limits on the Neutrino Mass and [83] J. Lesgourgues, G. Mangano, G. Miele, and S. Pastor,
Lifetime, JHEP 04 (2020) 020, [1909.05275]. Neutrino Cosmology. Cambridge University Press, 2,
[66] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, 2013.
Determining the Neutrino Lifetime from Cosmology, [84] T. Bertólez-Martı́nez, I. Esteban, R. Hajjar, O. Mena,
Phys. Rev. D 103 (2021), no. 4 043519, [2002.08401]. and J. Salvado. To appear (2024).
[67] M. Escudero, J. Lopez-Pavon, N. Rius, and [85] F. Capozzi, E. Lisi, A. Marrone, and A. Palazzo,
S. Sandner, Relaxing Cosmological Neutrino Mass Current unknowns in the three neutrino framework,
Bounds with Unstable Neutrinos, JHEP 12 (2020) 119, Prog. Part. Nucl. Phys. 102 (2018) 48–72,
[2007.04994]. [1804.09678].
[68] G. Franco Abellán, Z. Chacko, A. Dev, P. Du, [86] C.-P. Ma and E. Bertschinger, Cosmological
V. Poulin, and Y. Tsai, Improved cosmological perturbation theory in the synchronous and conformal
constraints on the neutrino mass and lifetime, JHEP Newtonian gauges, Astrophys. J. 455 (1995) 7–25,
08 (2022) 076, [2112.13862]. [astro-ph/9506072].
[69] N. Bellomo, E. Bellini, B. Hu, R. Jimenez, [87] W. Hu, Crossing the phantom divide: Dark energy
C. Pena-Garay, and L. Verde, Hiding neutrino mass in internal degrees of freedom, Phys. Rev. D 71 (2005)
modified gravity cosmologies, JCAP 02 (2017) 043, 047301, [astro-ph/0410680].
[1612.02598]. [88] J. Garriga and V. F. Mukhanov, Perturbations in
[70] G. Dvali and L. Funcke, Small neutrino masses from k-inflation, Phys. Lett. B 458 (1999) 219–225,
gravitational θ-term, Phys. Rev. D 93 (2016), no. 11 [hep-th/9904176].
113002, [1602.03191]. [89] W. Hu, Structure formation with generalized dark
[71] C. S. Lorenz, L. Funcke, E. Calabrese, and matter, Astrophys. J. 506 (1998) 485–494,
S. Hannestad, Time-varying neutrino mass from a [astro-ph/9801234].
supercooled phase transition: current cosmological [90] C. Nascimento, Accurate fluid approximation for
constraints and impact on the Ωm -σ8 plane, Phys. Rev. massive neutrinos in cosmology, Phys. Rev. D 108
D 99 (2019), no. 2 023501, [1811.01991]. (2023), no. 2 023505, [2303.09580].
[72] G. Dvali, L. Funcke, and T. Vachaspati, Time- and [91] J. Lesgourgues and T. Tram, The Cosmic Linear
Space-Varying Neutrino Mass Matrix from Soft Anisotropy Solving System (CLASS) IV: efficient
Topological Defects, Phys. Rev. Lett. 130 (2023), no. 9 implementation of non-cold relics, JCAP 09 (2011)
091601, [2112.02107]. 032, [1104.2935].
[73] C. S. Lorenz, L. Funcke, M. Löffler, and E. Calabrese, [92] J. Lesgourgues, The Cosmic Linear Anisotropy Solving
Reconstruction of the neutrino mass as a function of System (CLASS) I: Overview, 1104.2932.
redshift, Phys. Rev. D 104 (2021), no. 12 123518, [93] D. Blas, J. Lesgourgues, and T. Tram, The Cosmic
[2102.13618]. Linear Anisotropy Solving System (CLASS) II:
[74] M. Escudero, T. Schwetz, and J. Terol-Calvo, A Approximation schemes, JCAP 07 (2011) 034,
seesaw model for large neutrino masses in concordance [1104.2933].
with cosmology, JHEP 02 (2023) 142, [2211.01729]. [94] J. Lesgourgues, The Cosmic Linear Anisotropy Solving
[Addendum: JHEP 06, 119 (2024)]. System (CLASS) III: Comparision with CAMB for
[75] M. Sen and A. Y. Smirnov, Refractive neutrino LambdaCDM, 1104.2934.
masses, ultralight dark matter and cosmology, JCAP [95] R. K. Sachs and A. M. Wolfe, Perturbations of a
01 (2024) 040, [2306.15718]. cosmological model and angular variations of the
[76] Y. Farzan and S. Hannestad, Neutrinos secretly microwave background, Astrophys. J. 147 (1967)
converting to lighter particles to please both KATRIN 73–90.
and the cosmos, JCAP 02 (2016) 058, [1510.02201]. [96] S. Dodelson and F. Schmidt, Modern Cosmology, 2nd
[77] J. Alvey, M. Escudero, N. Sabti, and T. Schwetz, Edition. Academic Press, 4, 2020.
Cosmic neutrino background detection in [97] J. Torrado and A. Lewis, Cobaya: Code for Bayesian
 12

Analysis of hierarchical physical models, JCAP 05

(2021) 057, [2005.05290].
[98] J. Torrado and A. Lewis, “Cobaya: Bayesian analysis
in cosmology.” Astrophysics Source Code Library,
record ascl:1910.019, Oct., 2019.
[99] P. Motloch and W. Hu, Lensinglike tensions in the
P lanck legacy release, Phys. Rev. D 101 (2020), no. 8
083515, [1912.06601].
[100] G. Efstathiou and S. Gratton, A Detailed Description
of the CamSpec Likelihood Pipeline and a Reanalysis of
the Planck High Frequency Maps, 1910.00483.
[101] SPT-3G Collaboration, F. Ge et al., Cosmology
From CMB Lensing and Delensed EE Power Spectra
Using 2019-2020 SPT-3G Polarization Data,
2411.06000.
[102] SPT-3G Collaboration, B. A. Benson et al.,
SPT-3G: A Next-Generation Cosmic Microwave
Background Polarization Experiment on the South
Pole Telescope, Proc. SPIE Int. Soc. Opt. Eng. 9153
(2014) 91531P, [1407.2973].
[103] LiteBIRD Collaboration, A. I. Lonappan et al.,
LiteBIRD science goals and forecasts: a full-sky
measurement of gravitational lensing of the CMB,
JCAP 06 (2024) 009, [2312.05184].
[104] M. Shoji and E. Komatsu, Massive Neutrinos in
Cosmology: Analytic Solutions and Fluid
Approximation, Phys. Rev. D 81 (2010) 123516,
[1003.0942]. [Erratum: Phys.Rev.D 82, 089901 (2010)].
[105] M. Zaldarriaga and D. D. Harari, Analytic approach to
the polarization of the cosmic microwave background in
flat and open universes, Phys. Rev. D 52 (1995)
3276–3287, [astro-ph/9504085].
[106] J. M. Bardeen, Gauge-invariant cosmological
perturbations, Phys. Rev. D 22 (Oct, 1980) 1882–1905.
[107] H. Kodama and M. Sasaki, Cosmological Perturbation
Theory, Prog. Theor. Phys. Suppl. 78 (1984) 1–166.
[108] V. Mukhanov, Physical Foundations of Cosmology.
Cambridge University Press, 2005.
[109] W. Hu, N. Sugiyama, and J. Silk, The Physics of
microwave background anisotropies, Nature 386 (1997)
37–43, [astro-ph/9604166].
[110] A. Lewis and A. Challinor, Weak gravitational lensing
of the CMB, Phys. Rept. 429 (2006) 1–65,
[astro-ph/0601594].
[111] A. Lewis, Efficient sampling of fast and slow
cosmological parameters, Phys. Rev. D 87 (2013),
no. 10 103529, [1304.4473].
[112] A. Gelman and D. B. Rubin, Inference from Iterative
Simulation Using Multiple Sequences, Statistical
Science 7 (1992), no. 4 457 – 472.
 13

Appendix A: Evolution of neutrino perturbations

In this Appendix, we detail how we compute the effective neutrino sound-speed c2eff and the neutrino anisotropic
stress σ, as functions of the perturbations neutrino mass mPert.
P
ν . We then substitute
P their values in Eqs. (9) and (10),
that we solve numerically within the P public code CLASS to compute the impact of mPert.
ν on CMB observables.
Technically, c2eff and σ depend on mPert.
ν as well as on the gravitational potentials ϕ and ψ (in this Appendix,
as in the entire manuscript, we work in the Newtonian gauge). Physically, this captures the backreaction of gravity
onto neutrinos, in the same way that the background equation of state wν depends both on the neutrino mass and
on the scale factor a (see Eq. (3)). To compute this dependence for massive, non-interacting neutrinos; we start from
the collisionless Boltzmann equation for the perturbed neutrino distribution Ψ [86],
 
∂Ψ q d ln f0 ϵ dϕ
+ iµ Ψ = iµ ψ − . (A1)
∂x ϵ d ln q q dx

Here, x ≡ kη is the product p of Fourier wavenumber and conformal time, q ≡ ap is the comoving 3-momentum with
direction n̂, µ ≡ ⃗k · n̂, ϵ ≡ q 2 + m2 a2 is the comoving energy, and f0 is given by Eq. (4) in the main text. Following
Refs. [90, 104], this equation can be implicitly solved in terms of ϕ and ψ. After expanding in Legendre polynomials
and integrating by parts,
 Z x   
d ln f0 ′ ϵ ′ q ′ ′ ′
Ψℓ (x) = Ψ(0)jℓ (y(0, x)) + ϕ(0)jℓ (y(0, x)) − ϕ(x)jℓ (0) − dx ψ(x ) + ϕ(x ) jℓ (y(x , x)) , (A2)
d ln q 0 q ϵ

where
∞
X
Ψ(k, q, µ, x) ≡ (−i)ℓ (2ℓ + 1)Ψℓ (k, q, x)Pℓ (µ) , (A3)
ℓ=0

with Pℓ the Legendre polynomials; jℓ (x) the spherical Bessel functions; and
Z x2
q
y(x1 , x2 ) ≡ ′)
dx′ , (A4)
x1 ϵ(x

the distance traveled by neutrinos between times x1 /k and x2 /k, divided by the mode size. The super-horizon
adiabatic initial condition for Ψ, assuming radiation domination and ultrarelativistic neutrinos, is [86]

1 d ln f0
Ψ(x = 0) = ψ(η = 0) . (A5)
2 d ln q

These equations allow to explicitly compute Ψ(k, q, µ, η) as a function of the gravitational potentials
pand the neutrino
P Pert.
dispersion relation ϵ(q). The latter depends on the total perturbations mass mν , since ϵ(q) = q 2 + a2 (mPert.
ν )2
with mPert.
ν each individual perturbations neutrino mass. Once Ψ is known, c2
eff and σ can be computed from Eqs. (12)
and (14) to (17) in the main text.
 14

Appendix B: Consistency with the background and gauge-invariance of the effective sound speed

As discussed in Section II B in the main text, the fact that adiabatic perturbations behave as background as k → 0
may enforce a consistency relation between c2s , σ and wν . On top of that, an indeterminacy is introduced by c2s being
gauge-dependent. In this Appendix, we show that writing the perturbations equations Eqs. (9) and (10) in terms of
the effective sound-speed c2eff overcomes both issues.

1. Gauge invariance

Under a gauge transformation

η → η + α(⃗x, η) , (B1)
⃗ x, η) + ⃗ϵ (⃗x, η) ,
⃗x → ⃗x + ∇β(⃗ (B2)

matter perturbations transform as [86]

δ → δ − αρ′ /ρ , (B3)
2
θ → θ − αk , (B4)
δP → δP − αP ′ , (B5)
σ → σ. (B6)

As a first consequence, σ is gauge-invariant. Straightforward substitution of the gauge transformation in Eq. (12)
shows, using Eq. (11), that c2eff is also gauge-invariant.

2. Consistency with background

For adiabatic perturbations, super-horizon long-wavelength perturbations must behave as background at a slightly
different time [83, 86, 96]. That is,

lim (P + δP )(η) = P (η + δη) = P (η) + P ′ (η)δη , (B7)

k→0
′
lim (ρ + δρ)(η) = ρ(η + δη) = ρ(η) + ρ (η)δη . (B8)
k→0

This immediately implies that, for our perturbations formalism to be consistent with the background,

δPν P′ w′
lim c2s ≡ lim = ′ν = wν − a′ ν ≡ c2ad (B9)
k→0 k→0 δρν ρν 3 a (1 + wν )

must be fulfilled. It is straightforward to check that, if this condition holds, enforcing the background equation Eq. (2)
at a conformal time η + δη leads to Eq. (9) (assuming θ → 0, which follows from Eq. (15) as k → 0).
To write this condition as a condition on ceff , we invert Eq. (12)
a′
2 3 a (1 + wν )θ
c2s = c2eff + (c2eff − cad ) , (B10)
k2 δ
It follows from this equation that, if
′
3 a (1 + wν )θ
lim a 2 = −1 (B11)
k→0 k δ

consistency with background is satisfied for all values of c2eff and wν . One can check (see Eqs. (B3) to (B6)) that this
condition is gauge-invariant.
We now turn to explicitly checking that Eq. (B11) is indeed satisfied for adiabatic perturbations. In the conformal
Newtonian gauge, super-horizon adiabatic initial conditions in the radiation-domination regime for ultrarelativistic
neutrinos correspond to δ = −2ψ, θ = k 2 ηψ/2, wν = 1/3, a′ /a = η −1 [86]; which satisfy Eq. (B11). Equation (B11)
 15

is also maintained by super-horizon evolution: using Eqs. (9) and (10) with k → 0 (where θ → 0 and σ → 0, see
Eqs. (15) and (17)), as well as the Friedmann and perturbed Einstein equations [86], we obtain
" ′ #  " a′ #
d 3 aa (1 + wν )θ a′
  
3 Ptot 3 a (1 + wν )θ
+δ = − 1+ + 3wν +δ , (B12)
dη k2 a 2 ρtot k2

with Ptot and ρtot the total pressure and energy density of the Universe. Since the right-hand side vanishes for the
′
initial conditions, super-horizon evolution keeps the identity true at all times, even if aa or wν change with time.
Finally, our definition of σ is also consistent with the background. Adiabatic perturbations behave as background
at super-horizon scales, which implies a diagonal stress-energy tensor, i.e., σ → 0 as k → 0. Equations (17) and (A2)
trivially fulfill this, because, j2 (0) = 0. That is, the anisotropic stress is non-zero only when the distance travelled by
neutrinos is of the order of the mode size. Since σ is a gauge-invariant quantity, this is true in all gauges.
 16

Appendix C: Detailed physical effects on the CMB anisotropies

P Backg. P Pert.
In this Appendix, we provide more details and explicit checks of the different effects of mν and mν on
the CMB anisotropies; with an emphasis on the underlying physical processes.
Unless stated otherwise, all the CMB power spectra shown throughout this manuscript are obtained by fixing the
standard ΛCDM parameters to their best-fit values of the Planck 2018 TT, TE, EE + lowE + lensing analysis [40].
The values of these parameters
P are 100θs = 1.04172, ωb = 0.02237, ωcdm = 0.1200, log(1010 As ) = 3.044, ns =
0.9649, τreio = 0.0544, and mν = 0.06 eV.

1. Modified Sachs-Wolfe separation

The amplitude of the ℓ-th multipole of CMB temperature anisotropies, without including CMB lensing, is given by
Z
1 dk
Cℓ = PR (k)[Θγℓ (η0 , k)]2 , (C1)
2π 2 k
with PR (k) the primordial curvature spectrum—in the standard scenario As = (k/k0 )ns −1 [83]—, and Θγℓ (η0 , k) the
present photon temperature transfer function, normalized so that the curvature perturbation at initial times R = 1.
To gain physical intuition, it is customary to write Θγℓ in terms of a line-of-sight integral [105]
Z η0 n o
′
g(η) (Θγ0 + ψ) + g(η)k −2 θB + e−τ (ϕ′ + ψ ′ ) jℓ [k (η0 − η)] dη ,

Θγℓ (η0 , k) = (C2)
ηin

with ηin the initial conformal time, θB the baryon velocity divergence,
Z η0
τ (η) ≡ a(η ′ )ne (η ′ )σT dη ′ (C3)
η

the photon optical depth, which changes from 0 well before recombination to 1 after recombination; and g(η) ≡ [e−τ (η) ]′
the visibility function, which is sharply peaked around recombination [83]. The Bessel function jℓ sets ℓ ∼ k(η0 − η).
Equation (C2) provides a simple physical interpretation of the observed CMB anisotropies. The first term can be
understood as the temperature anisotropy at recombination, redshifted by the local gravitational potential (Sachs-
Wolfe effect, SW); the second term as a Doppler shift due to baryon velocities at recombination; and the third term
as the accumulated gravitational redshift as CMB photons travels from recombination to us (Integrated Sachs-Wolfe
effect, ISW). The ISW effect can be further understood as caused by the depth of gravitational potential wells changing
while photons are inside them, due to the combined action of structure growth and the expansion of the Universe.
This separation, however, is not particularly appropriate to study the impact of neutrinos. Both (Θγ0 + ψ) and
(ϕ′ + ψ ′ ) depend on the neutrino anisotropic stress σ even at super-horizon scales [86]. Thus, the SW contribution
to modes that are much larger than the horizon would depend on σ; and, as σ changes with time, there would be
an ISW effect even for modes that are well outside the horizon. This would contradict the expected super-horizon
behavior of adiabatic perturbations: as k → 0, these perturbations behave as background and their evolution cannot
depend on perturbation-related quantities such as σ.
Indeed, one can numerically check that the ℓ → 0 limit of Eq. (C2) is independent of the anisotropic stress around
recombination, even if the SW and ISW terms separately are not. Technically, the SW term can be integrated by
parts, leading to terms proportional to e−τ that may as well be interpreted as an ISW contribution. Physically, as
argued by Bardeen in his seminal paper on gauge-invariant cosmological perturbations [106], there is an ambiguity of
what one means by a temperature or metric perturbation at scales comparable to or bigger than the horizon.
To facilitate physical understanding, below we redefine the separation between SW (understood as the temperature
fluctuation at recombination, redshifted by the local gravitational potential) and ISW (understood as the accumulated
gravitational redshift from recombination until today) contributions, so that they both depend only on background
quantities in the k → 0 limit.
To such purpose, we split the last term in Eq. (C2) as
   
6(1 + wtot ) 2  a  1 + 3wtot 2a ′ 6
ϕ+ψ = ϕ+ ψ + ′ ϕ′ − ϕ−ψ+ ϕ , (C4)
5 + 3wtot 3(1 + wtot ) a 5 + 3wtot 3 a′ 1 + 3wtot
with wtot ≡ Ptot /ρtot the total equation of state of the Universe. We identify the comoving curvature perturbation
a′ θtot 2  a ′
R≡ϕ+ = ϕ + ψ + ϕ , (C5)
a k2 3(1 + wtot ) a′
 17

with θtot the total velocity divergence of all the matter content in the Universe. Here, we have used the Einstein
equations to express θtot in terms of ϕ and ψ. Thus, we can rewrite Eq. (C4) as
 
6(1 + wtot ) 1 + 3wtot 2a ′ 6
ϕ+ψ = R− ϕ−ψ+ ϕ . (C6)
5 + 3wtot 5 + 3wtot 3 a′ 1 + 3wtot
It immediately follows from the Einstein equations that R is conserved on super-horizon scales.
Substituting Eq. (C6) in the third term in Eq. (C2), and integrating by parts the second term in Eq. (C6), we
obtain
Z η0 n o
−2
′
g(η)ΘSW θB jℓ [k (η0 − η)] + e−τ ΘISW

Θγℓ (η0 , k) = γℓ + g(η)k γℓ dη , (C7)
ηin

with
  
1 + 3wtot 2a ′ 6
ΘSW
γℓ ≡ Θγ0 + ψ + ϕ−ψ+ ϕ jℓ [k (η0 − η)] =
5 + 3wtot 3 a′ 1 + 3wtot
(C8)
a′ θtot
  
1 + 3wtot
= R + Θγ0 + jℓ [k (η0 − η)] ,
5 + 3wtot a k2
where we have used the Einstein equations to go from the first to the second line; and
 ′  
6(1 + wtot ) 1 + 3wtot 4 a ′
ΘISW
γℓ ≡ R jℓ [k(η 0 − η)] − ϕ − ψ + ϕ k jℓ′ [k(η0 − η)] . (C9)
5 + 3wtot 5 + 3wtot 1 + 3wtot a′
The first term in Eq. (C8) is, on super-horizon scales, explicitly σ-independent (we normalize the transfer function
to R = 1 at initial times). Since on such scales R is conserved, it only depends on the background quantity wtot —
the amplitude of super-horizon perturbations does depend on the background equation of state [107, 108]—, and it
provides the correct ℓ → 0 SW plateau of the CMB power spectrum [95, 96] (see Eq. (C10) below). The second
term in Eq. (C8) is zero for super-horizon scales and adiabatic perturbations (see Appendix B 2). It corresponds to
one of the gauge-invariant density perturbations defined by Bardeen [106], and it represents the photon temperature
perturbations in a gauge comoving with matter.
The first term in Eq. (C9) is also, on super-horizon scales, explicitly σ-independent. It only introduces an ISW
effect if wtot changes, which would affect the super-horizon gravitational potentials. The second term in Eq. (C9) is
explicitly suppressed as k → 0.
Figure C1 showsPour separationP among SW and ISW contributions to the CMB temperature anisotropies, for
mBackg. Pert.
P
different values of ν and mν . The blue line has the same background as m ν = 0, and the SW
contribution is consequently suppressed as ℓ → 0. The purple and red lines also share the same background evolution,
but they have sub-percent differences as ℓ → 0. This is reasonable, as these multipoles correspond to scales where
(kℓ=2 ηrec )2 ∼ 0.004, so they are not completely out of the horizon around recombination.

2. Physical effects on the SW contribution

The SW effect accounts for temperature anisotropies atPrecombination,Predshifted by the local gravitational poten-
tial. We now break down the different physical effects of mBackg.
ν and mPert.
ν on the SW contribution in Fig. C1.

a. Low multipoles

Low multipoles (ℓ ≲ 10) correspond to modes that are super-horizon at recombination. Their evolution only
′
depends on background quantities, and it can be analytically estimated. Setting k ≪ aa in Eq. (C8), the ℓ → 0
plateau is given by [96]
 2
SW 1 + 3wtot (ηrec )
lim ℓ(ℓ + 1)Cℓ ≈ 8 As , (C10)
ℓ→0 5 + 3wtot (ηrec )
which coincides with the standard result ℓ(ℓ + 1)CℓSW ≈ 8As /25 for matter domination, wtot (ηrec ) = 0 [95, 96].
P Backg.
Increasing mν reduces the neutrino contribution to wtot , reducing the ℓ → 0 limit of the SW contribution.
P Pert.
There is also a sub-percent effect of mν that, as explained above, is due to the corresponding modes not being
completely out of the horizon at recombination.
 18

P
FIG. C1. Background- and perturbations-induced impact of mν on the Sachs-Wolfe (SW, left) and Integrated Sachs-Wolfe
(ISW, right) effects as defined in Eqs. (C8) and (C9), without including CMB lensing. The SW contribution can be understood
as due to the temperature fluctuations at recombination, redshifted by the local gravitational potential. In turn, the ISW
contribution can be understood as the accumulated gravitational redshift from recombination until today. Top plots show the
contribution to the total CMB anisotropy power spectrum, and bottom plots the relative difference with respect to massless
neutrinos. The SW contribution, as defined in this work, only depends on background quantities as k → 0.

b. Intermediate multipoles

Intermediate multipoles (10 ≲ ℓ ≲ 500) correspond to modes that are comparable to the horizon at recombination.
Their evolution is largely influenced by the gravitational potentials ϕ and ψ. It can be approximated as [109]
 
Θγ0 ψ
meff Θ′′γ0 + k 2 ≃ −meff k2 − ϕ′′ , (C11)
3 3

where meff ≡ 1 + ρργb +P

+Pγ ≃ 1 + 3ρb /(4ργ ); with ρb , Pb , ργ , and Pγ the baryon energy density, baryon pressure, photon
b

energy density, and photon pressure, respectively. Equation (C11) is the equation of a forced harmonic oscillator,
with the gravitational potentials playing the role of an external force.
Decaying gravitational potentials enhance the amplitude of oscillations via Eq. (C11). Physically, large initial
potentials force the fluid into a highly compressed state. If they then decay, photon pressure overcomes gravity and
the photon-baryon fluid oscillates with a larger amplitude. On top of that, decaying gravitational potentials reduce
the redshift experienced by CMB photons as they leave the last-scattering surface, further increasing the anisotropies.
For a detailed explanation of these effects, we refer to the work of Ref. [109]. Gravitational potentials decay in the
radiation-dominated era [96], and since recombination happens soon after matter-radiation equality, they are still
decaying when SW anisotropies get frozen (see Fig. C2).
P Backg.
The effect of neutrino masses is then straightforward. A large mν increases the expansion rate of the
Universe, boosting the decay of gravitational potentials, P Pert.as shown in Fig. C2; and increasing the SW contribution
to CMB anisotropies, as shown in Fig. C1. A large mν enhances neutrino clustering above the free-streaming
scale, slowing down the decay of gravitational potentials, as shown in Fig. C2; and decreasing the SW contribution
to CMB anisotropies,
P Backg.as shown
P inPert.
Fig. C1. These effects are opposite, and thus the net effect for standard massive
neutrinos with mν = mν is smaller, as shown in Figs. C1 and C2.
P Pert.
Figure C3 shows the scale dependence of the clustering effect induced by mν . As described in Section II B in
the main text, neutrino clustering is a scale-dependent “kinematic”P effect that reflects that neutrinos do not move at
thePspeed of light. The characteristic scale, kFS , is proportional to mPert.
ν , see Eq. (18). As can be seen in Fig. C3,
Pert.
if mν increases, the aforementioned depletion of CMB anisotropies affects higher multipoles.
 19

Recombination Recombination
Radiation domination is still close,
Mass in perturbations:
so potentials are decaying.
neutrino clustering
slows down the decay.
Mass only in
perturbations
Mass only in
background

Mass in background:
a faster expansion rate
boosts the decay.

(Mpc) (Mpc)

FIG. C2. Gravitational potential ψ (left) and its derivative with respect to conformal time (right) for a mode entering the
horizon (kηrec ∼ 1) at the time of recombination. We normalize to an initial comoving curvature perturbation R = 1. This mode
corresponds to ℓ ∼ 100. Background and perturbations neutrino masses affect intermediate-ℓ CMB anisotropies by boosting and
slowing down the decay of ψ, respectively.

P Pert. P Pert.
FIG. C3. Scale-dependent impact of increasing mν on the SW contribution to CMB anisotropies. Larger mν
slow down the decay of gravitational potentials through neutrino clustering (see Fig. C2), reducing the amplitude of acoustic
oscillations at recombination. Since neutrino clustering is a scale-dependent effect controlled by the free-streaming length, kFS ,
Eq. (18), larger perturbations masses increase kFS and propagate the impact to higher multipoles.

c. High multipoles

High multipoles (ℓ ≳ 500) correspond to modes which enter the horizon much before recombination. As discussed
in Section III B in the main text, these modes are affected by diffusion damping before recombination, an SW effect
controlled by the damping scale
P θBackg.
D , Eq. (22). SinceP
θD only depends on background
P Backg.quantities, the high-ℓ contribution
in Fig. C1 depends only on mν and not on mPert.
ν . An increased mν boosts the expansion rate of the
Universe, reducing θD and enhancing the high-ℓ anisotropies in Fig. C1. We have checked that all the enhancement
is due to a modified θD , because it can be completely removed by artificially keeping θD fixed (technically, this is
achieved by changing the primordial Helium fraction Yp , which rescales ne in Eq. (22) [83]).
 20

3. Physical effects on the ISW contribution

The ISW contribution accounts for the accumulated gravitational redshift of CMB photons on their way to Earth.
Physically, when photons enter gravitational potential wells they get redshifted, and when they exit them they get
blueshifted. If the depth of the gravitational potential wells is constant, both effects compensate each other. A
nonzero ISW effect thus requires time-dependent gravitational potentials. This can only be the case if the Universe
is not matter-dominated [96].
The ISW effect is only relevant at multipoles ℓ ≲ 200, as shown in the top right panel of Fig. C1. At smaller scales,
gravitational potentials experience a long period of radiation domination inside the horizon, and they are suppressed.
It can be split into an early and a late contribution.
The early contribution is due to the residual radiation left after recombination, and it affects modes that were inside
the horizon at that time, 10 ≲ ℓ ≲ 200. As discussed above (see Fig. C2), background neutrino P
P Backg. masses boost the
decay of gravitational potentials and perturbations masses slow it down. Therefore, mν and mPert.
ν produce
an enhancement and a depletion of the ISW effect, respectively; as can be seen in Fig. C1.
The late contribution is due to the late-time epoch of accelerated expansion caused by the cosmological constant Λ,
and it affects modes that entered the horizon recently, ℓ ≲ 10. Λ makes the gravitational potentials decay, inducing a
late ISW effect. As explained in Section III B in the main text,Pthe main neutrino-mass effect is indirect: the angular
scale of CMB peaks θs , Eq. (21), is very well measured, and mBackg.
ν would modify it by changing the expansion
rate of the Universe. To compensate for this effect, H0 gets P reduced, reducing Λ and the late ISW effect as can be seen
in Fig. C1. There is also a subleading reduction caused by mPert.
ν , as neutrino clustering
P Backg.slows down the Λ-induced
decay of potentials. This effect is small, see Figs. 3 and C1, and it is only present if mν ̸= 0, because otherwise
the neutrino energy density is too diluted at late times and neutrinos do not affect gravitational potentials.

4. CMB lensing

On their way to Earth, CMB photons get randomly deflected by the gravitational pull of the large-scale structure
of the Universe [96, 110]. This weak gravitational lensing is a second-order non-linear effect, since it is a perturbative
deflection of perturbative anisotropies, but with the precision of current data it is detectable. Lensing has two
observable features on CMB anisotropies, that are more prominent at high multipoles. First, it smoothes out the
power spectrum, since it mixes photons coming from different points in the last scattering surface. Second, it transfers
power from large scales to small scales, due to its non-linear nature. This leads to increased anisotropies at high ℓ.
Neutrino masses affect the evolution of gravitational potentials, as explained above and shown in Fig. C2. Con-
sequently, they modify CMB lensing. On the one hand, background neutrino masses increase the expansion rate
of the Universe and boost the decay of gravitational potentials. This reduces CMB lensing. On the other hand,
perturbations masses make neutrinos cluster below the free-streaming scale and slow down the decay of potentials.
This enhances CMB lensing P in Backg.
a scale-dependent
P Pert. way.
Figure C4 shows how mν and mν control these effects. The left panel shows CMB anisotropies
without lensing,
P Backg. while the right panel includes lensing. The high-ℓ tail of the unlensed power spectrumP
P Backg. only depends
on mν through diffusion damping, as explained in Appendix C 2. In turn, since m
P Backg.ν and mPert.
ν have
opposite effects on lensing, they produce opposite features at lensed high-ℓ. As mν reduces lensing, CMB
anisotropies are less smoothed, which is visible as wiggles in phase with the CMB peaks. Weaker P Pert. lensing also implies
less power transferred to high multipoles, which leads to a slight depletion at ℓ ≳ 2000. As mν increases lensing,
its effect is opposite, reducing the amplitude of the wiggles
P Backg. and slightly enhancing anisotropies in the ℓ ≳ 2000 region.
Perturbations mass effects are only relevant if mν ̸= 0, because otherwise the neutrino energy density is too
diluted at late times and neutrinos doPnot affect gravitational potentials.
As
P mentioned above, the impact of mPert. ν is scale-dependent, with the characteristic scale kFS being proportional
Pert.
P Pert.
to mν (see Eq. (18)). Figure C5 shows this effect and how mν controls it. We plot the power spectrum
of the lensing potential, which is directly related to gravitational potentials [96, 110] (the scale dependence in CMB
anisotropies is less evident,
P Pert.because lensing is non-linear and thePtranslation between k and ℓ is not straightforward).
As the figure shows, mν enhances lensing, and increasing mPert.
ν shifts its strongest impact to higher ℓ.
Finally, for completeness and to understand thePrelevance of lensing in setting CMB limits, Fig. C6 shows the
predictions for Fig. 5 without lensing. A value of mBackg.
ν that is excluded by our analysis predicts
P Backg.an unlensed
damping tail well within error bars. Therefore, the constraining power of Planck 2018 data on mν at high ℓ
comes from its effect on CMB lensing.
 21

Lensed

P Backg. P Pert.
FIG. C4. Impact of mν and mν on CMB anisotropies,
P Backg. with (left) and without (right) the effect of
Pgravitational
lensing.
P Pert. At high ℓ, unlensed anisotropies only depend on m ν through diffusion damping. However,
P Backg. both mBackg.
ν and
mν P have an effect on the lensed anisotropies, and they do so in opposite directions. mν reduces CMB lensing,
while mPert.
ν enhances it.

P Pert. P Pert.
FIG. C5. Scale-dependent impact of mν on the CMB lensing power spectrum. mν increases gravitational potentials
through neutrino clustering, P
increasing CMB lensing. Neutrino clustering is a scale-dependent effect, controlled by the scale
kFS (see Eq. (18)). A larger mPert.
ν increases kFS , shifting the impact of neutrino clustering to higher multipoles.
 22

FIG. C6. Impact on CMB anisotropies of parameters excluded by our analysis, without gravitational lensing. This Figure is
analogous to Fig. 5 in the main text. Since the predictions are unlensed,
P Backg. but data is lensed, here the Planck18 data only shows
visually the uncertainty of the measurements. CMB bounds on mν are mainly due to the effect of lensing at high-ℓ.
 23

Appendix D: Full results of the statistical analysis

In this Appendix, we provide the full results of our statistical analysis. We carry out a Bayesian analysis with the
Markov Chain Monte Carlo code COBAYA [97, 98], as mentioned in the main text.
Table I contains the priors on the cosmological parameters over which we scan.

log(1010 As ) mBackg. mPert.

P P
Parameter ns 100θs ωb ωcdm τreio ν [eV] ν [eV]

Prior U[1.61, 3.91] U[0.8, 1.2] U[0.5, 10] U[0.005, 0.1] U[0.001, 0.99] U[0.01, 0.8] U[0.0, 3.0] U[0.0, 3.0]

TABLE I. Cosmological parameters that we scan over, and their corresponding priors. U[a, b] denotes a uniform distribution
with lower limit a and upper limit b.

Figure D7 shows the 1D posterior probabilities and 2D credible regions for all parameters in our analysis. We also
include the posteriors on the derived parameters H0 and σ8 , as they show the largest degeneracies with neutrino
masses. All posteriors are well-contained within their priors, and the convergence of the MCMC run is determined
by an R − 1 < 0.02 Gelman-Rubin test [111, 112].
PTable II summarizes
P Pert. the 1–3σ credible intervals on the total neutrino mass, both for the standard scenario
( mBackg.
ν = mν ) and for our analysis that splits among background and perturbations neutrino-mass effects.
P Backg.
Both in PFig. D7 and Table II, there is a high degree of consistency when comparing mν with Planck 2018
results on mν . This is one of our main results, i.e., that the CMB is much more sensitive to the P Backg. dilution of
temporal
the neutrino energy density
P than to the “kinematic” effects of neutrino masses. The limits on
P Pert. mν are slightly
weaker than those on mν due to the degeneracy P Pert.with m ν discussed in
P Pert. the main text. We also note from Figs. 4
and D7 and Table II that the best fit for m ν is non-zero, with m ν = 0 being excluded at ∼ 1σ. This is
driven by the lensing anomaly ( mPert.
P
ν increases CMB lensing, see Section III in the main text and Appendix C 4
above), but the result is not statistically significant.

Planck 2018; TT, TE, EE+lowE+lensing

Mass type 68% CL 95% CL 99% CL

P
mν < 0.11 eV < 0.24 eV < 0.35 eV

mBackg.
P
ν < 0.13 eV < 0.29 eV < 0.40 eV

0.40+0.19
P Pert.
mν −0.29 eV < 0.79 eV < 0.97 eV

TABLE II. Neutrino-mass limits at different confidence

P Backg.levels.
P These results correspond to the posteriors shown in Fig. D7.
mPert.
P
mν refers to the standard scenario in which mν = ν .
 24

FIG. D7. CMB 1σ and 2σ credible regions for all parameters in our analysis, together with individual posterior probabilities.
Dashed lines correspond to the standard results, where neutrino-mass effects are not split between background and perturba-
tions. In each subfigure, unshown parameters are marginalized over.