# Reconstruction of the neutrino mass as a function of redshift

Christiane S. Lorenz<EMAIL_ADDRESS>Institute for Particle Physics
and Astrophysics, ETH Zürich, Wolfgang-Pauli-Strasse 27, CH-8093 Zürich,
Switzerland  Lena Funcke<EMAIL_ADDRESS>Perimeter Institute
for Theoretical Physics, 31 Caroline Street North, Waterloo, Ontario, N2L 2Y5,
Canada  Matthias Löffler  Seminar for Statistics, Department of Mathematics,
ETH Zürich, Rämistrasse 101, CH-8092 Zürich, Switzerland  Erminia Calabrese
School of Physics and Astronomy, Cardiff University, The Parade, Cardiff, CF24
3AA, United Kingdom

(Received ; published – 00, 0000)

###### Abstract

We reconstruct the neutrino mass as a function of redshift, $z$, from current
cosmological data using both standard binned priors and linear spline priors
with variable knots. Using cosmic microwave background temperature,
polarization and lensing data, in combination with distance measurements from
baryonic acoustic oscillations and supernovae, we find that the neutrino mass
is consistent with $\sum m_{\nu}(z)=$ const. We obtain a larger bound on the
neutrino mass at low redshifts coinciding with the onset of dark energy
domination, $\sum m_{\nu}(z=0)<1.46$ eV (95% CL). This result can be explained
either by the well-known degeneracy between $\sum m_{\nu}$ and
$\Omega_{\Lambda}$ at low redshifts, or by models in which neutrino masses are
generated very late in the Universe. We finally convert our results into
cosmological limits for models with non-relativistic neutrino decay and find
$\sum m_{\nu}<0.21$ eV (95% CL), which would be out of reach for the KATRIN
experiment.

##  I Introduction

Cosmological surveys and particle physics experiments are independent and
complementary probes of neutrino properties. Neutrino oscillation experiments
have measured the squared mass differences between neutrino mass eigenstates,
giving a lower bound of 59 meV for the total sum of the neutrino masses, $\sum
m_{\nu}$ Group  _et al._ (2020). In addition, the Karlsruhe Tritium Neutrino
Experiment (KATRIN) has constrained the electron neutrino mass to be lower
than $m_{\nu,e}<0.8$ eV at 90% CL Aker  _et al._ (2021a, b). Independent from
these constraints, measurements of the cosmic microwave background (CMB) with
the Planck satellite mission, combined with distance measurements from
baryonic acoustic oscillations (BAO) from the Sloan Digital Sky Survey (SDSS)
and 6dF have provided a tight upper bound on the total sum of neutrino masses,
$\sum m_{\nu}<120$ meV at 95% CL Aghanim  _et al._ (2018).

Contrary to neutrino mass direct detection limits, the cosmological neutrino
mass bound assumes a specific cosmological model, usually the $\Lambda$Cold
Dark Matter ($\Lambda$CDM) model or its single-parameter extensions. In more
extended cosmological models, the cosmological neutrino mass limits can become
less stringent due to opening wider parameter spaces and/or covering more
complex physics scenarios.

On the one hand, the lower bound of $\sum m_{\nu}=59$ meV imposed by neutrino
oscillation experiments can be relaxed in cosmological analyses to $\sum
m_{\nu}=0$ meV, for example, if cosmological neutrinos disappear in the late
Universe. Although the original “neutrinoless Universe” proposal Beacom  _et
al._ (2004) has been ruled out111By free-streaming of the cosmic neutrino
background before photon decoupling Hannestad (2005); Lancaster  _et al._
(2017), by the resulting phase shift in the CMB peaks Follin  _et al._ (2015),
and by precise CMB measurements of the effective number of species in the
early Universe Aghanim  _et al._ (2018)., other scenarios like models
predicting postrecombination neutrino mass generation and subsequent relic
neutrino annihilation Dvali and Funcke (2016a) are still possible. In such
models, the neutrino mass parameter cannot be captured using cosmological data
and can only be measured using terrestrial and astrophysical experiments such
as KATRIN. The possible cosmological disappearance of the neutrino mass
parameter has also been proposed in the context of modified gravity theories
Bellomo  _et al._ (2017); Hagstotz  _et al._ (2019a, b); these should soon be
tested with surveys like Euclid Hagstotz  _et al._ (2019b). If we allow for
strong fine tuning, another possibility to eliminate the cosmological neutrino
mass bounds would be to postulate a new light scalar particle that couples to
neutrinos, with a coupling constant that needs to be smaller than $g\sim
10^{-7}$ to avoid laboratory constraints Esteban and Salvado (2021).

On the other hand, the upper cosmological neutrino mass bound is sensitive to
a number of model assumptions and can be slightly relaxed, for example when
the dark energy equation of state is allowed to vary in time Joudaki and
Kaplinghat (2012); Joudaki (2013); Allison  _et al._ (2015); Lorenz  _et al._
(2017); Yang  _et al._ (2017); Vagnozzi  _et al._ (2018a); Roy Choudhury and
Choubey (2018); Roy Choudhury and Hannestad (2019), when the curvature of the
Universe is not fixed Roy Choudhury and Hannestad (2019), when considering
additional relativistic degrees of freedom Aghanim  _et al._ (2018), or when
assuming non-standard momentum distributions of the cosmic neutrinos Oldengott
_et al._ (2019). Moreover, it has been shown that the cosmological neutrino
mass bound can be substantially weakened when neutrinos are unstable and thus
their lifetime is smaller than the age of the Universe Escudero and Fairbairn
(2019); Chacko  _et al._ (2020a, b); Escudero  _et al._ (2020); Barenboim  _et
al._ (2020) or when neutrino masses are varying in time Fardon  _et al._
(2004); Dvali and Funcke (2016a); Koksbang and Hannestad (2017); Lorenz  _et
al._ (2019). Note, however, that the strongly relaxed neutrino mass bounds of
$\sum m_{\nu}$$<0.9$ eV (95% CL) with neutrino decays Chacko  _et al._ (2020a)
or $\sum m_{\nu}$$<4.8$ eV (95% CL) with time-varying neutrino masses Lorenz
_et al._ (2019) have been derived from previous releases of cosmological data,
including the Planck 2015 CMB data. The neutrino mass bound has been
substantially tightened with the Planck 2018 release, $\sum m_{\nu}$=$0.12$ eV
(95% CL, Planck 2018 TTTEEE+lowE+lensing+BAO Aghanim  _et al._ (2018)) and
therefore, we expect the above-mentioned bounds also to change.

In a previous publication Lorenz  _et al._ (2019), some of us investigated how
the standard cosmological neutrino mass bound would be affected if neutrino
masses were to be generated late in the Universe. This analysis suggested that
the combination of current CMB temperature, polarization and lensing data, as
well as BAO and supernovae (SN) data prefers neutrino masses to be generated
at low redshifts, allowing a significantly larger cosmological neutrino mass
upper bound. An alternative interpretation of the results is that current
cosmological datasets do not necessarily require neutrinos to be massive,
pushing the scale factor of the phase transition in the model to very low
redshifts. The trend of larger neutrino mass bounds in the late Universe had
already been noted before (see, e.g., Refs. Brookfield  _et al._ (2006a, b);
Koksbang and Hannestad (2017); Battye and Moss (2014); Wyman  _et al._ (2014);
Beutler  _et al._ (2014); Poulin  _et al._ (2018) and references therein).
This could either arise due to the well-known degeneracy between neutrino
masses and dark energy, or be explained by new physics beyond the Standard
Model of Cosmology. In particular, it seems striking that the energy scales of
dark energy and neutrino masses are numerically very close,
$\sqrt[4]{\rho_{\Lambda}}\sim m_{\nu}\sim{\rm meV}$. If cosmological data
permit larger neutrino masses during dark energy domination, there could be an
intriguing theoretical connection between these two phenomena (see, e.g.,
Refs. Fardon  _et al._ (2004); Bjaelde  _et al._ (2008); Ayaita  _et al._
(2016); Mandal  _et al._ (2019); D’Amico  _et al._ (2019); Dvali and Funcke
(2016a) and references therein).

The tight neutrino mass limit from cosmology is making a direct detection from
Tritium $\beta$-decay experiments like KATRIN very challenging. Indeed,
KATRIN’s target sensitivity of 200 meV at 90% CL for the $\nu_{e}$ mass
Drexlin  _et al._ (2013) will unlikely hit the cosmological limit obtained in
either $\Lambda$CDM or its simple extensions Aghanim  _et al._ (2018); Roy
Choudhury and Hannestad (2019). Assuming all observations and analysis
assumptions are correct and free of systematic effects, we are entering a
regime where a direct detection of neutrino masses with particle detectors
such as KATRIN could become a strong hint for nonstandard neutrino physics. At
the same time, the next generation of cosmological surveys aim to improve
their reach in neutrino mass sensitivity and to make the first detection
Abazajian  _et al._ (2015); Mishra-Sharma  _et al._ (2018).

To work toward these future goals, we present here a model-independent
approach to investigate the possibility that neutrino masses change on
cosmological timescales, reconstructing the neutrino mass as a function of
redshift.

A similar methodology has been applied to reconstruct possible variations in
other cosmological parameters, in particular dark energy parameters Holsclaw
_et al._ (2010, 2011); Seikel  _et al._ (2012a); Seikel and Clarkson (2013);
Alberto Vazquez  _et al._ (2012); Montiel  _et al._ (2014); Hee  _et al._
(2016, 2017); Poulin  _et al._ (2018); Gerardi  _et al._ (2019); Joudaki  _et
al._ (2018); Keeley  _et al._ (2019); Park  _et al._ (2021); Bonilla  _et al._
(2021) and the Hubble parameter Shafieloo  _et al._ (2012); Seikel  _et al._
(2012b); Bernal  _et al._ (2016); Poulin  _et al._ (2018); Gómez-Valent and
Amendola (2018); Renzi and Silvestri (2020); Bernal  _et al._ (2021) as a
function of redshift, the shape of the primordial power spectrum as a function
of wave number Bridle  _et al._ (2003); Guo  _et al._ (2011); Hazra  _et al._
(2014); Hunt and Sarkar (2014); Lévy  _et al._ (2020), and parameters
describing deviations from general relativity as a function of redshift and
wave number Garcia-Quintero  _et al._ (2020). Exploring variations of
cosmological parameters in time also gains in importance in light of tensions
between parameters inferred from current high- and low-redshift data (see,
e.g., Ref. Verde  _et al._ (2019)).

In general, reconstructions fully accounting for degeneracies and the
interplay between cosmological parameters, such as the dark matter density and
dark energy parameters, are especially hard to achieve Kunz (2009); Amendola
_et al._ (2013); Busti and Clarkson (2016). In the case of the neutrino mass,
a model-independent reconstruction is particularly challenging. A
comprehensive model for time-varying neutrino masses would require an
interaction with other energy sectors, such as dark radiation or dark energy,
to satisfy energy conservation laws. Additionally, including such interactions
permits to fully capture and exploit the physics signatures of the model.
While in the case of dark radiation this interaction is negligible for many
cosmological scenarios (e.g., the dark radiation resulting from neutrino
decays has negligible cosmological impact Chacko  _et al._ (2020a)), in the
case of dark energy this interaction could alter $w_{\rm de}(z)$ and thus
generate multiple signatures that allow us to place stronger limits on the
model.

In this paper, we take a conservative approach and use a generic function to
model the neutrino mass sum. We also ignore potential additional constraining
power coming from the inclusion of the coupling between the neutrino and the
dark sector which is theory specific, and perform a model-independent analysis
that only focuses on neutrinos without assuming any specific interactions with
other dark sectors. This yields a conservative approach in the modelling but
not necessarily the most conservative neutrino mass constraints (see the
discussion in Sec. V). The method used here expands on the work done in Ref.
Lorenz  _et al._ (2019), which assumed a specific time-varying neutrino mass
model, and also spans models with neutrino decays for which we will set new
limits.

The paper is structured as follows. In Sec. II we explain the theoretical
background of cosmological neutrino mass constraints and neutrino mass models.
In Sec. III we present our methodology, in particular the different datasets
and reconstruction methods. We present our results in Sec. IV and summarize
and discuss them in Sec. V.

##  II Theoretical background

###  II.1 Neutrino mass constraints from cosmological probes

At the beginning of their cosmic journey, neutrinos are relativistic particles
and behave as a radiation component in the early Universe. When their kinetic
energy term drops below the mass term due to cooling in an expanding Universe,
neutrinos start to behave as nonrelativistic, massive particles. The redshift
of this transition is inversely proportional to the neutrino mass Ichikawa
_et al._ (2005): neutrinos with a smaller mass become nonrelativistic at a
later time compared to neutrinos with a higher mass. Neutrino masses thus
affect cosmological observations during different cosmic epochs, which in turn
allows us to probe $\sum m_{\nu}$ at different redshifts.

_The cosmic microwave background–_ Constraining neutrino properties with the
CMB has been a rich research area with extensive literature (see, e.g., Refs.
Archidiacono  _et al._ (2017); Lesgourgues  _et al._ (2013); Vagnozzi  _et
al._ (2017)). To summarize, the key effects we look for in CMB probes are: (i)
Effects on the background evolution of the Universe via changes in the angular
diameter distance at recombination, $D_{A}(z_{rec})$ and the Hubble parameter
Archidiacono  _et al._ (2017); Lesgourgues  _et al._ (2013). In general, there
is a strong degeneracy between neutrino parameters and the Hubble constant
Archidiacono  _et al._ (2017); Howlett  _et al._ (2012); Sutherland (2018);
Moresco  _et al._ (2012); Hall and Challinor (2012).

(ii) Effects on the evolution of perturbations. In particular, the integrated
Sachs-Wolfe effect is affected, both at early times during radiation
domination (eISW) Archidiacono  _et al._ (2017); Lesgourgues  _et al._ (2013),
as well at late times when dark energy starts to dominate the evolution of the
Universe ($\ell$ISW) Cabass  _et al._ (2015). The eISW depends on the
relativistic degrees of freedom, $N_{\mathrm{eff}}$, as well on neutrino
masses, but in a different way. Whereas $N_{\mathrm{eff}}$ mostly changes the
amplitude of the first peak in the CMB anisotropy power spectrum, neutrino
masses affect the amplitude of the eISW on a large range of multipoles
depending on the neutrino mass and the corresponding free-streaming scale
Archidiacono  _et al._ (2017); Lesgourgues  _et al._ (2013); Hou  _et al._
(2014). The $\ell$ISW effect dominates at $\ell<30$ and is therefore elusive
because of cosmic variance in CMB data alone. However, when cross-correlated
with galaxy number counts, this effect is a promising avenue for measuring
neutrino masses Lesgourgues  _et al._ (2008).

(iii) Effects on the matter distribution deflecting the CMB photons by
gravitational lensing Lewis and Challinor (2006). This will leave both an
imprint on the CMB temperature and polarization anisotropies (in the
high-$\ell$ region of the spectra), as well as generate a CMB lensing
convergence signal, $C_{\ell}^{\phi\phi}$. The latter probe will capture the
small-scale suppression of the matter power spectrum due to large neutrino
thermal velocities and corresponding free-streaming out of density
fluctuations Lesgourgues  _et al._ (2013); Lesgourgues and Pastor (2006).
Depending on whether the neutrino wavelength is above or below the free-
streaming wavelength, neutrinos cluster as cold dark matter and baryons, or
free-stream out of gravitational wells. This slows down the clustering of
matter, leading to a suppression of the matter power spectrum on the
corresponding scales which is more pronounced for larger neutrino masses Hu
_et al._ (1998). The CMB damping tail is measured with high precision with
current data Henning  _et al._ (2018); Aghanim  _et al._ (2020a); Choi  _et
al._ (2020) and CMB lensing is now in a high signal-to-noise regime Bianchini
_et al._ (2020); Darwish  _et al._ (2020).

(iv) The optical depth to reionization $\tau$ is degenerate with the amplitude
of scalar primordial fluctuations $A_{s}$. Since the amount of clustering of
cosmic structures is tightly linked to $\sum m_{\nu}$, this in turn becomes a
strong degeneracy between $\tau$ and $\sum m_{\nu}$. CMB polarization
measurements of $\tau$ enable us to obtain a tighter constraint on $\sum
m_{\nu}$ Allison  _et al._ (2015); Calabrese  _et al._ (2017).

_Baryonic acoustic oscillations–_ Before CMB decoupling, baryons and photons
are tightly coupled to each other. At the time of recombination, the CMB
photons decouple from the baryons, and the oscillations of the baryon-photon
fluid are frozen in the CMB anisotropies. In addition to the CMB, these
oscillations also leave a characteristic imprint in the large scale structure
of the Universe, both transverse as well as along the line of sight Eisenstein
_et al._ (2005); Cole  _et al._ (2005). In particular, the size of the BAO is
known, and therefore BAO can be used as standard rulers to infer either
$H(z)r_{s}(z^{*})$ or $D_{A}/r_{s}(z^{*})$, where $D_{A}(z)$ is the angular
diameter distance (see, e.g., Refs. Weinberg  _et al._ (2013); Bassett and
Hlozek (2009) for reviews) and $r_{s}$ the sound horizon at decoupling
$z^{*}$. These quantities are in particular sensitive to the matter density
$\Omega_{m}h^{2}$. BAO measurements allow to constrain the energy contribution
of massive neutrinos to the matter density. The BAO feature has been detected
in the clustering of galaxies Eisenstein  _et al._ (2005); Cole  _et al._
(2005); Alam  _et al._ (2020), as well as in the clustering of low-redshift
quasars de Carvalho  _et al._ (2018); Ata  _et al._ (2018); Hou  _et al._
(2021), and in the correlations of Lyman-$\alpha$ systems McDonald and
Eisenstein (2007); Bautista  _et al._ (2017) (see below). This latter probe
provides an additional measurement to standard BAO data and extends BAO
observations toward high redshifts.

_Lyman- $\alpha$ forest–_ The absorption lines of neutral hydrogen in the
intergalactic medium (IGM) in quasar spectra are sensitive to cosmological
parameters, and probe cosmic structure formation at redshifts between $z\sim
2-6$. In particular, Lyman-$\alpha$ forest measurements characterize small
structures on the scale of sub-Mpc to Mpc. This anchors the level of the
matter power spectrum on scales between $k\sim 0.1-2$ and therefore probes the
regime where the suppression of the matter power spectrum due to massive
neutrinos is the most pronounced (see, e.g., Refs. Seljak  _et al._ (2005,
2006); Palanque-Delabrouille  _et al._ (2015a, b, 2020)). The relevant summary
statistic is the 1D flux power spectrum, which is related to the matter power
spectrum using a nonlinear transformation Croft  _et al._ (1998); Zhan  _et
al._ (2005). The level of the 1D flux power spectrum depends both on the
amplitude of the linear matter power spectrum and on the sum of neutrino
masses, resulting in a degeneracy between these two parameters Pedersen  _et
al._ (2020). This limits the ability of current Lyman-$\alpha$ forest
measurements to constrain neutrino masses. The current upper bound from
Lyman-$\alpha$ measurements alone is $\Sigma m_{\nu}<0.71$ eV (95% CL)
Palanque-Delabrouille  _et al._ (2020). Probing cosmological parameters with
the Lyman-alpha forest is also challenging because of observational and
astrophysical systematics, which need to be modeled in the analyses Hui and
Gnedin (1997); Walther  _et al._ (2020).

_Supernovae–_ Measurements of the luminosity distance $D_{L}(z)$ from
supernovae explosions are not directly sensitive to neutrino properties.
However, they can be used to constrain dark energy parameters, such as the
dark energy density $\Omega_{\Lambda}$ and the dark energy equation of state
$w_{\mathrm{de}}(z)$ (see, e.g., Ref. Goobar and Leibundgut (2011) for a
review), which helps significantly to break degeneracies with neutrino
parameters.

Combining the different cosmological probes described above enhances
significantly the constraining power coming from only one of the datasets
(see, e.g., Refs. Archidiacono  _et al._ (2017); Abazajian  _et al._ (2015);
Vagnozzi  _et al._ (2017); Boyle (2019); Aghanim  _et al._ (2018)). This is
due to the fact that (i) different probes are capturing the physics of the
Universe at different redshifts, and (ii) different probes depend on different
parameter combinations, which follow different degeneracy directions. For
example, Ref. Archidiacono  _et al._ (2017) describes in detail how the
combination of CMB and BAO helps to break the $H_{0}-\sum m_{\nu}$ degeneracy.

###  II.2 Neutrino mass models

Most of the possible neutrino mass models, in particular the ones arising from
the famous seesaw mechanism Minkowski (1977); M. Gell-Mann, P. Ramond, and R.
Slansky (edited by P. Van Nieuwenhuizen and D. Freedman, Elsevier, Amsterdam,
1979); Yanagida (1980); Mohapatra and Senjanović (1980); Schechter and Valle
(1980, 1982), cannot be tested using cosmological data. From a cosmological
perspective, more focus is then naturally placed on studying neutrino mass
mechanisms that make cosmologically testable predictions, such as models
providing neutrino mass variations on cosmologically interesting timescales.

As mentioned earlier, when allowing for the neutrino mass to vary as a
function of redshift, an interaction with the dark energy or dark radiation
sector is required to satisfy energy conservation laws. The most popular
models that couple the neutrino and dark energy sectors are mass varying
neutrino (MaVaN) scenarios (see, e.g., Ref. Fardon  _et al._ (2004)). Here,
the neutrinos couple to a light scalar field, which slowly rolls in a flat
potential. Originally, this direct link between neutrino masses and dark
energy was proposed to explain the similarity between the energy scales of
neutrino masses and quintessence-like dark energy ($E\sim 10^{-3}$ eV).
However, the coupling to a light scalar mediates an attractive force between
the neutrinos and leads to bound state formation Afshordi  _et al._ (2005),
such that the light scalar field can only explain dark energy under rather
special circumstances (see, e.g., Refs. Bjaelde  _et al._ (2008); Ayaita  _et
al._ (2016); Mandal  _et al._ (2019)). For example, the Growing Neutrino
Quintessence model studied in Ref. Ayaita  _et al._ (2016) yields time-
dependent neutrino masses, which vanish in the early Universe and raise from
$m_{\nu}=0$ eV at $a\lesssim 0.2$ to $m_{\nu}\sim 0.7$ eV at $a=1$.

Coupling the mass-varying neutrino sector to dark radiation is less trivial
and usually occurs with simultaneously coupling both sectors to the dark
energy sector. For example, it has been proposed that a simultaneous variation
of neutrino masses, a light sterile neutrino fraction, and dark energy can be
used to test certain aspects of the Weak Gravity Conjecture D’Amico  _et al._
(2019).

Another gravitational avenue to couple neutrinos to dark energy and dark
radiation is the gravitational neutrino mass model proposed in Ref. Dvali and
Funcke (2016a), which served as a motivation for the previous cosmological
study in Ref. Lorenz  _et al._ (2019). The key cosmological prediction of this
model is that the cosmological neutrino mass parameter vanishes. Indeed, the
model predicts massless neutrinos in the early Universe and the generation of
neutrino masses in a late-time cosmological phase transition at $T\lesssim
m_{\nu}$. In the phase transition, neutrino masses are generated, followed by
neutrino decay into the lightest mass eigenstate and rapid annihilation into
massless Goldstone bosons (see Refs. Dvali and Funcke (2016a, b); Funcke  _et
al._ (2020) for more details). Thus, the model predicts massless neutrinos in
the early Universe and massless dark radiation in the late Universe. The
intermediate regime with massive neutrinos, which exists directly after the
phase transition, exists only for cosmologically negligible timescales. This
implies that the cosmological neutrino mass parameter would be zero, as
currently preferred by cosmological data Aghanim  _et al._ (2018). This also
implies that experiments aiming at a direct detection of the relic neutrino
background, such as the proposed PTOLEMY experiment Betti  _et al._ (2019),
would not be able to detect this background,222At first sight, one might
expect that PTOLEMY could detect the massless Goldstone bosons, which are
neutrino-composite bosons (similar to the light quark-composite mesons in
QCD). However, the boson’s energy is $E=T_{\nu}<m_{\nu}$, while the neutrino
capture would release the other neutrino of the bound state, requiring an
energy of at least $m_{\nu}$. This would violate energy conservation, unless
one of the neutrinos is almost massless. We thank Pedro Machado for bringing
up this argument. unless we allow for substantial neutrino asymmetries (see
below).

If the model in Ref. Dvali and Funcke (2016a) is extended by allowing for
neutrino asymmetries (i.e., more neutrinos than antineutrinos or vice versa),
not all neutrinos would find an antineutrino partner to annihilate, leaving
behind a fraction of relic neutrinos with a nonzero cosmological mass
parameter. As Ref. Lorenz  _et al._ (2019) demonstrated, even in this case,
the cosmological neutrino mass bound would be substantially weakened to $\sum
m_{\nu}<4.8$ eV (95% CL). However, we note that this weakened bound was
obtained with Planck 2015 data and is expected to become more stringent when
including the Planck 2018 dataset. This is because the $\Lambda$CDM constraint
on $\sum m_{\nu}$ strengthened by a factor of $\mathcal{O}(2)$ from $\sum
m_{\nu}=0.21$ eV (95% CL, Planck 2015 without polarization Ade  _et al._
(2016)) to $\sum m_{\nu}=0.12$ eV (95% CL, Planck 2018 Aghanim  _et al._
(2018)).

Finally, neutrinos can couple to the dark radiation sector if they are
unstable, with lifetimes shorter than the age of the Universe. In particular,
2-body decays of neutrinos into BSM particle species, such as massless or very
light sterile neutrinos and Goldstone bosons, have been proposed to relax
cosmological neutrino mass bounds. Early cosmological studies based on Planck
2015 data implied a relaxation of the cosmological neutrino mass bound up to
$\sum m_{\nu}\lesssim 0.9$ eV Chacko  _et al._ (2020a, b); Escudero  _et al._
(2020). As above, this bound is expected to become tighter when including
Planck 2018 data (see Sec. IV.3). Thus, neutrino decay can only slightly
alleviate cosmological bounds, similar to other extensions of the $\Lambda$CDM
model, such as dynamical dark energy. However, such models leave the
intriguing possibility of observing less (or even completely vanishing)
neutrino mass in the late Universe, similar to the model proposed in Ref.
Dvali and Funcke (2016a), and thus make cosmological observations crucial for
determining the neutrino mass origin and beyond-SM neutrino interactions.

These and other theoretical motivations to study neutrino mass variations on
cosmologically interesting timescales, which typically yield nontrivial
interactions with other cosmological sectors such as dark energy and dark
radiation, are the main rationale for the work that we present here.

##  III Methodology

Dataset | Redshift range
---|---
Planck 2018 CMB TTTEEE Aghanim  _et al._ (2018, 2020a) | mostly $z=1100$
Planck 2018 CMB lowl Aghanim  _et al._ (2018, 2020a) | mostly $z=1100$
Planck 2018 CMB lowE Aghanim  _et al._ (2018, 2020a) | mostly $z=8$ and $z=1100$
Planck 2018 CMB lensing Aghanim  _et al._ (2018, 2020b) | $0\leq z\leq 1100$
BAO (6dF) Beutler  _et al._ (2011) | $z=0.106$
BAO (SDSS DR7 BOSS MGS) Ross  _et al._ (2015) | $z=0.15$
BAO (SDSS DR12 BOSS) Alam  _et al._ (2017) | $z=0.38,0.51,0.61$
BAO (SDSS DR14 eBOSS quasars)Ata  _et al._ (2018) | $z=1.52$
BAO (SDSS DR14 eBOSS Ly-$\alpha$) Blomqvist  _et al._ (2019) | $z=2.34$
BAO (SDSS DR14 eBOSS cross Ly-$\alpha$-QSO) Blomqvist  _et al._ (2019) | $z=2.35$
SN (Pantheon) Scolnic  _et al._ (2018) | $0.01<z<2.3$
Table 1:  List of datasets used in the analyses performed in this work and
their corresponding redshift range.

###  III.1 Data

Motivated by the expected contribution highlighted in Sec. II.1, we include
the following datasets in our analysis:

  1. 1.)

CMB and CMB lensing: We use CMB temperature, polarization and lensing data
from the Planck 2018 data release Aghanim  _et al._ (2018).

  2. 2.)

BAO: We use BAO from 6dF Beutler  _et al._ (2011), the Sloan Digital Sky
Survey (SDSS) DR7 Main Galaxy Sample (MGS) Ross  _et al._ (2015) and the SDSS
Baryon Oscillation Spectroscopic Survey (BOSS) twelfth data release (DR12)
Alam  _et al._ (2017). The mild discrepancies seen between Planck and Lyman-
alpha BAO data have decreased in recent releases of the Lyman-alpha BAO from
SDSS DR14 eBOSS Blomqvist  _et al._ (2019); Schöneberg  _et al._ (2019).
Therefore we include also this additional BAO dataset which provide two data
points at $z\sim 1.5$ and $z\sim 2.3$. We also include the BAO measurement
from quasars from the 14th data release of the extended BOSS (eBOSS) quasar
sample Ata  _et al._ (2018), giving us another data point at $z\sim 1.5$.

  3. 3.)

Supernovae: We additionally add type IA supernova data from from the Pantheon
Supernovae Sample Scolnic  _et al._ (2018) in order to break the degeneracies
between dark energy and neutrino masses.

We demonstrate the impact of these specific datasets later in Sec. IV, and
summarize their respective redshift range in Tab. 1.

###  III.2 Reconstruction method

We modify the publicly available Einstein-Boltzmann code CAMB Lewis  _et al._
(2000) and the corresponding Monte-Carlo Markov chain package CosmoMC Lewis
and Bridle (2002) in order to implement different reconstruction methods for
redshift-dependent neutrino masses.

####  III.2.1 Binned reconstruction

As a first attempt, we parametrize the neutrino mass with a step function for
which it is assumed that the neutrino mass has a constant positive amplitude
in each redshift bin. In the statistics literature this is called a
regressogram Wasserman (2006). To maximally exploit our datasets we use a
parametrization with six redshift intervals:

$\sum m_{\nu}(z)=\begin{cases}\sum m_{\nu,0}&\text{($0\leq z<z_{1}$)}\\\ \sum
m_{\nu,1}&\text{($z_{1}\leq z<z_{2}$)}\\\ \sum m_{\nu,2}&\text{($z_{2}\leq
z<z_{3}$)}\\\ \sum m_{\nu,3}&\text{($z_{3}\leq z<z_{4}$)}\\\ \sum
m_{\nu,4}&\text{($z_{4}\leq z<z_{5}$)}\\\ \sum m_{\nu,5}&\text{($z\geq
z_{5}$)}.\end{cases}$ (1)

We choose the edges of the redshift bins to pinpoint specific transitions in
the composition of the Universe and to highlight the impact of using different
cosmological probes. We set $z_{1}=0.5$, $z_{2}=3$, $z_{3}=10$, $z_{4}=100$
and $z_{5}=1100$. The first bin explores $\sum m_{\nu}$ during dark energy
domination, the second bin spans the BAO interval, the third bin includes most
of the remaining information expected in CMB lensing Carbone  _et al._ (2008),
the fourth bin covers intermediate redshifts, and the second to the last bin
stretches out to the time of CMB decoupling. The last bin for $z>1100$
captures all the integrated information of the pre-recombination Universe.
This choice of bins will also allow us to set constraints on neutrino decay
models Escudero and Fairbairn (2019); Chacko  _et al._ (2020a, b); Escudero
_et al._ (2020), as described in Sec. IV.3. The neutrino masses in the
individual bins are five additional parameters in the model, compared to the
standard $\Lambda$CDM model with massive neutrinos. The discontinuity in the
binned $\sum m_{\nu}(z)$ leads only to small discontinuities in related
quantities, such as the Hubble rate $H(z)$, which however cause no numerical
instabilities in the analysis. A similar parametrization has recently been
used for the reconstruction of the dark energy equation of state
$w_{\mathrm{de}}(z)$ from gamma-ray bursts Muccino  _et al._ (2020).

####  III.2.2 Spline priors

To perform a smoother fit and to potentially identify features in the neutrino
mass sum that are hidden in the binned reconstruction, we also consider
Bayesian regression splines with variable knot points Dimatteo and Kass
(2001).

The knots correspond to the bin margins of the regressogram prior function
described above, and model change points where the trajectory of the neutrino
mass might change its slope. Previous literature in cosmology has often used
splines with fixed knot positions (see, e.g. Guo  _et al._ (2011); Bernal  _et
al._ (2016); Poulin  _et al._ (2018)). This, however, requires to choose the
positions of the knots in advance and can therefore significantly influence or
bias the result of the reconstruction. Here, we estimate the position of the
knots from the data and include them as free parameters in the analysis. Using
variable knot positions yields a more flexible fit than in the case with fixed
knots and allows the fit to adapt to underlying features of the model. This
methodology has previously successfully been applied several times in
cosmology Vázquez  _et al._ (2012); Alberto Vazquez  _et al._ (2012); Millea
and Bouchet (2018), most recently in Ref. Bernal  _et al._ (2021).

Compared to binned priors, the resulting reconstruction is smooth and not
piecewise constant anymore. We model two knots, $z_{1}$ and $z_{2}$, and
linearly interpolate between $\sum m_{\nu,z_{0}}$, $\sum m_{\nu,z_{1}}$, $\sum
m_{\nu,z_{2}}$ and $\sum m_{\nu,z_{3}}$. For our analysis, we choose $z_{0}=0$
and $z_{3}=1100$. In that case, the neutrino mass at redshift $z$ is given by

$\displaystyle\sum m_{\nu}(z)$ (2) $\displaystyle=$
$\displaystyle\begin{cases}\sum m_{\nu,z_{0}}+(\sum m_{\nu,z_{1}}-\sum
m_{\nu,z_{0}})\frac{z}{z_{1}}&\text{($z_{0}\leq z<z_{1}$)}\\\ \sum
m_{\nu,z_{1}}+(\sum m_{\nu,z_{2}}-\sum
m_{\nu,z_{1}})\frac{z-z_{1}}{z_{2}-z_{1}}&\text{($z_{1}\leq z<z_{2}$)}\\\ \sum
m_{\nu,z_{2}}+(\sum m_{\nu,z_{3}}-\sum
m_{\nu,z_{2}})\frac{z-z_{2}}{z_{3}-z_{2}}&\text{($z_{2}\leq z<z_{3}$)}\\\ \sum
m_{\nu,z_{3}}&\text{($z\geq z_{3}$)}.\end{cases}$

We then choose a uniform prior on the logarithm of the corresponding redshift
with $-1\leq\log_{10}(z_{i})_{i=1,2}\leq 3.041$, covering the entire redshift
region until recombination and focusing on low redshifts where most of the
data is located. We further need to impose that $z_{2}$ has to be equal or
larger than $z_{1}$, in accordance to the definition above in Eq. (2). In
addition, we choose $z_{1,2}>0.1$, corresponding to the lowest redshift of the
BAO dataset (see Table 1).

With this parametrization, we have five additional free parameters compared to
the standard $\Lambda$CDM with massive neutrinos, three for $\sum m_{\nu}$ at
$z_{1}$, $z_{2}$ and $z_{3}$, and two for the positions of $z_{1}$ and
$z_{2}$. In order to obtain pointwise credible bands for $\sum m_{\nu}(z)$, we
compute the neutrino mass sum for two hundred points in $z$ as derived
parameters with GetDist Lewis (2019) to sample well the entire redshift range.
We then compute pointwise the mean and credible intervals for each point in
$z$.

##  IV Constraints from current data

For our neutrino mass reconstruction, we run MCMC chains with our modified
version of CosmoMC. Since both our methods only touch the neutrino mass
modeling, we vary the standard $\Lambda$CDM parameters alongside with the new
parameters of the reconstruction: the cold dark matter density
$\Omega_{c}h^{2}$, the baryon density $\Omega_{b}h^{2}$, the scalar spectral
index $n_{s}$, the amplitude of primordial fluctuations $A_{s}$ and the
optical depth to reionization $\tau$. In addition, we have six additional
parameters for the reconstruction with fixed bins ($\sum m_{\nu,0}$, $\sum
m_{\nu,1}$, $\sum m_{\nu,2}$, $\sum m_{\nu,3}$, $\sum m_{\nu,4}$ and $\sum
m_{\nu,5}$), and six for the reconstruction with linear splines and variable
knots ($\sum m_{\nu,z_{0}}$, $\sum m_{\nu,z_{1}}$, $\sum m_{\nu,z_{2}}$, $\sum
m_{\nu,z_{3}}$ and the position of the knots $\log_{10}(z_{1})$ and
$\log_{10}(z_{2})$). We choose a uniform prior between [0:5] for $\sum
m_{\nu}$ in the individual redshift bins (similarly to the Planck 2018
analysis Aghanim  _et al._ (2018)), and we do not assume a specific
correlation between the values of $\sum m_{\nu}$ in the different bins,
leaving them to vary independently. In addition, we assume a normal neutrino
mass hierarchy in line with current cosmological constraints Aghanim  _et al._
(2018); Gariazzo  _et al._ (2018); Roy Choudhury and Hannestad (2019). The
choice of the mass hierarchy should not significantly affect our final results
presented in Sec. IV, as current cosmological data cannot (yet) distinguish
between the different neutrino mass hierarchies Jimenez  _et al._ (2010);
Gerbino  _et al._ (2017); Vagnozzi  _et al._ (2017); Lattanzi and Gerbino
(2018); Gariazzo  _et al._ (2018); De Salas  _et al._ (2018); Roy Choudhury
and Hannestad (2019); Mahony  _et al._ (2020); Archidiacono  _et al._ (2020);
Xu  _et al._ (2020); Stöcker  _et al._ (2020); Hergt  _et al._ (2021).

###  IV.1 Binned reconstruction

We show the results of the binned neutrino mass reconstruction in Fig. 1, and
in Table 2. We plot both the 95% and the 68% CL upper limits for each mass
parameter, as well as the 95% limit for $\sum m_{\nu}(z)={\rm const.}$ (i.e.,
the standard single parameter extension of $\Lambda$CDM for massive neutrinos)
obtained from the same data combinations. We always plot the 95% CL and 68% CL
one-tail upper limits, e.g., we have removed the 32% and 5% highest samples
for these limits. For consistency we choose the same type of limits for all
the individual redshift bins, and can therefore compare better the results for
the individual bins and data combinations. For the full data combination
including CMB, CMB lensing, BAO and SN, we find

$\displaystyle\sum m_{\nu}(0\leq z<0.5)$ $\displaystyle<1.13\text{ eV}$
$\displaystyle\sum m_{\nu}(0.5\leq z<3)$ $\displaystyle<0.42\text{ eV}$
$\displaystyle\sum m_{\nu}(3\leq z<10)$ $\displaystyle<0.37\text{ eV}$
$\displaystyle\sum m_{\nu}(10\leq z<100)$ $\displaystyle<0.19\text{ eV}$
$\displaystyle\sum m_{\nu}(100\leq z<1100)$ $\displaystyle<0.32\text{ eV}$
$\displaystyle\sum m_{\nu}(z\geq 1100)$ $\displaystyle<0.40\text{ eV}\hskip
10.00002pt\text{(95\% CL)}.$

These results are consistent with neutrino masses constant in time. In
addition, we observe that the neutrino mass bound becomes less stringent at
low redshifts ($z\leq 3$), and at very high redshifts ($z\geq 1100$). The
large neutrino mass bound at low redshifts could arise due to the well-known
degeneracy between the neutrino mass sum and the dark energy density (see
below). It could also point at a very late generation of neutrino masses in
the Universe, and would therefore be consistent with the trend seen in Ref.
Lorenz  _et al._ (2019). In the Standard Model of Cosmology, the neutrino
masses are generated during the electroweak or earlier phase transitions.

Figure 1: Reconstructed upper bounds of the neutrino mass as function of
redshift. The panels from top to bottom show incremental addition of data
starting from Planck 2018 CMB temperature and polarization, and then
subsequently adding Planck CMB lensing, BAO from BOSS DR12, 6dF and MGS, eBOSS
DR14 quasars and Lyman-alpha, and SN from Pantheon 18. The left panels show
the results obtained with the binned parameterization and the right panels the
limits from the reconstruction with linear splines and variable knots. In all
cases, we also report with the black solid line the 95% CL constraint from the
same data combination for $\sum m_{\nu}(z)={\rm const.}$ (i.e., the standard
single parameter extension of $\Lambda$CDM for massive neutrinos).

Our plots also allow us to note the contribution to the constraint from
different probes.

_Top row:_ In the top row of Fig. 1, we plot the $z$-dependent constraints on
the neutrino mass sum from CMB anisotropy data alone. Because the binned
reconstruction introduces more free parameters, the constraints in all bins
are consistent but systematically larger than the standard neutrino mass bound
from the latest Planck release Aghanim  _et al._ (2018); Group  _et al._
(2020)333Note that the $\sum m_{\nu}(z)=$ const. line should be compared to
the 95% shaded region.. The CMB anisotropy data put a strong constraint on
$\sum m_{\nu}$ at high redshifts, but only a weak constraint at low redshifts.
We also note that the neutrino mass sum in the lowest redshift bin is
unconstrained. In this bin the CMB data alone provides insufficient
constraining power. As described earlier, at low redshifts, the CMB anisotropy
data would mostly be sensitive to neutrino masses via the late integrated
Sachs-Wolfe effect at low multipoles $\ell$, which is cosmic variance limited
and subject to significant parameter degeneracies. In the lowest redshift bin
($0<z<0.5$), the CMB has only little constraining power, and the prior (with
prior mean $\sum m_{\nu}=2.5$ eV) dominates the result of the reconstruction.
At intermediate redshifts, the data strongly constrain the neutrino mass
through the lensing of the CMB anisotropies which will be picked up here as a
smoothing of the temperature and polarization data, and the impact of
neutrinos on the background evolution (see the discussion in Sec. II). At high
redshifts close to recombination, we observe a small increase in the
uncertainty of $\sum m_{\nu}$. This might be due to the fact that at
recombination the constraint on $\sum m_{\nu}$ comes mostly from the early
integrated Sachs-Wolfe effect, the CMB damping tail, and background effects,
while all other effects from neutrino masses only become relevant at
intermediate to low redshifts.

_Second row:_ In the second row in Fig. 1, we plot the constraints on $\sum
m_{\nu}(z)$ from CMB anisotropies plus the CMB lensing reconstruction. The CMB
lensing kernel peaks around $z\sim 2$ (see, e.g., Refs. Hassani  _et al._
(2016); Manzotti (2018)) and carries constraining power on the bin between
$0.5\leq z<3$. Other constraints are also improved by an indirect effect: CMB
lensing is also sensitive to the total matter density and the dark energy
density Sherwin  _et al._ (2011), and therefore helps to constrain $\sum
m_{\nu}(0<z<0.5)$ by partially breaking the tridimensional degeneracy between
$\Omega_{\Lambda}$, $\Omega_{m}$, and $\sum m_{\nu}$. This can also be seen in
Fig. 2, where we show the results for $\sum m_{\nu}$ in the different bins, as
well as its correlation with $\Omega_{\Lambda}$ and $\Omega_{m}$.

_Third row:_ The BAO data points used in this work are distributed over the
whole redshift range between $z=0.1-2.35$. This mostly improves the
constraints in the first two neutrino mass bins. In addition, the inclusion of
BAO reduces significantly the degeneracy between $\Sigma m_{\nu}$ and
$\Omega_{m}$, and therefore also improves the overall constraints in all
redshift bins.

_Bottom row:_ We include SN data from the Pantheon 18 data set to constrain
further the dark energy component. We observe the smallest uncertainty of
$\sum m_{\nu}$ in the bin between $10\leq z\leq 100$. The sensitivity of the
neutrino mass bound starts to weaken at a similar redshift as the onset of
dark energy domination (see below). The degeneracy between the sum of neutrino
masses and the dark energy density might play an important role here.

Dataset | CMB | CMB + CMBL | CMB + CMBL + BAO | CMB + CMBL + BAO + SN
---|---|---|---|---
$\sum m_{\nu}(z\leq 0.5)$ (eV) | 4.86 | 3.07 | 1.73 | 1.13
$\sum m_{\nu}(0.5\leq z\leq 3)$ (eV) | 3.19 | 2.85 | 0.51 | 0.42
$\sum m_{\nu}(3\leq z\leq 10)$ (eV) | 1.64 | 1.63 | 0.40 | 0.37
$\sum m_{\nu}(10\leq z\leq 100)$ (eV) | 0.47 | 0.50 | 0.20 | 0.19
$\sum m_{\nu}(100\leq z\leq 1100)$ (eV) | 0.39 | 0.37 | 0.33 | 0.32
$\sum m_{\nu}(z\geq 1100)$ (eV) | 0.44 | 0.42 | 0.41 | 0.40
Table 2:  95% upper limits on the amplitudes of the neutrino mass at different
redshifts reconstructed with a binned parametrization and using different data
combinations. Figure 2: 2-dimentional contours showing the 68% and 95% CL and
1-dimensional posteriors for the neutrino mass parameters in the individual
bins, $\Omega_{\Lambda}$ and $\Omega_{m}$. Different colors show different
data combinations, as in Fig. 1.

We note in Fig. 2 that in all cases there is little correlation between the
individual neutrino mass amplitudes. We also studied how the inclusion of
growth rate measurements obtained with redshift-space distortions (RSD) Alam
_et al._ (2017) would affect our results presented here. RSD can be used to
probe the growth rate $f(k,z)$ at different scales and redshifts Percival and
White (2009). The scale-dependent feature in this quantity is a source of
information for neutrino masses Hernández (2017); Villaescusa-Navarro  _et
al._ (2018); Marulli  _et al._ (2011); Boyle and Komatsu (2018). Including the
growth rate measurements, we found a decrease of the bound of $\sum m_{\nu}$
of 38% in the lowest redshift bin and only a marginal improvement in other
bins. The exact modeling of RSD is under active development and therefore we
present our final results without including RSD.

###  IV.2 Reconstruction with splines

Next we reconstruct the neutrino mass sum with linear splines and variable
knots (see Sec. III). We investigate at which redshift the sensitivity on
$\sum m_{\nu}$ starts to weaken by allowing the knots $z_{1}$ and $z_{2}$ to
vary as well. We show our results for this reconstruction on the right-hand
side of Fig. 1. We also present our results in Table 3. For the final data
combination we find

$\displaystyle\sum m_{\nu}(z=0)$ $\displaystyle<1.46\text{ eV}$
$\displaystyle\sum m_{\nu}(z=1100)$ $\displaystyle<0.53\text{ eV}\hskip
30.00005pt\text{(95\% CL)}.$

Dataset | CMB | CMB + CMBL | CMB + CMBL + BAO | CMB + CMBL + BAO + SN
---|---|---|---|---
$\sum m_{\nu}(z=0)$ (eV) | 4.75 | 2.91 | 1.43 | 1.46
$\sum m_{\nu}(z=0.5)$ (eV) | 3.50 | 2.74 | 2.31 | 0.76
$\sum m_{\nu}(z=3)$ (eV) | 2.11 | 1.84 | 0.18 | 0.18
$\sum m_{\nu}(z=10)$ (eV) | 1.07 | 1.11 | 0.16 | 0.15
$\sum m_{\nu}(z=100)$ (eV) | 0.69 | 0.47 | 0.17 | 0.18
$\sum m_{\nu}(z=1100)$ (eV) | 0.52 | 0.50 | 0.53 | 0.53
Table 3: 95% upper limits on the amplitudes of the neutrino mass at different
redshifts reconstructed with linear splines and variable knots and using
different data combinations. Figure 3: The results for the change points
$z_{1}$ and $z_{2}$ from the reconstruction with linear splines and variable
knots. Figure 4: Final result: Neutrino mass limits as function of redshift
obtained from Planck 2018 CMB temperature, polarization, and lensing data, BAO
data from BOSS DR12, 6dF, and MGS, eBOSS DR14 quasars, and Lyman-alpha, and SN
data from Pantheon 18 (zoom of the lower right-hand panel in Fig. 1). The
onset of dark energy domination is denoted with a dashed violet line. The
dashed red line shows $\sum m_{\nu}$ in the case when the neutrino with the
largest neutrino mass becomes nonrelativistic at redshift $z$, assuming a
normal neutrino mass hierarchy.

In general, when comparing the results from the two different reconstruction
methods, we note that the constraints for $\sum m_{\nu}$ at different
redshifts are similar. There are small differences, in particular for $\sum
m_{\nu}(z=0)$, because in the binned reconstruction $\sum m_{\nu}$ is fixed
between the bin margins. We also note that, as expected, we obtain smooth
curves for the reconstruction with linear splines. This is due to the fact
that the positions of the knots $z_{1}$ and $z_{2}$ are free parameters in the
analysis, smoothing the posterior means and credible bands. Compared to the
binned reconstruction, we observe that the credible bands of the reconstructed
curve via splines with variable knots exhibit more features despite the same
number of free parameters.

For constant neutrino masses in time, we would, in principle, expect flat
posteriors for the two knots $z_{1}$ and $z_{2}$, as the position of the knots
does not matter when fitting a constant function. However, even in that case,
the posterior for $z_{1}$ and $z_{2}$ can still vary because of possible
degeneracies between the neutrino mass sum and other cosmological parameters
and because of the different amount of data points at different redshifts,
resulting in local changes of the uncertainty of the neutrino mass sum.

In Fig. 3, we see that there is no clear preference for the position of the
knots $z_{1}$ and $z_{2}$. For the position of the first knot (change point)
we find $z_{1}<0.89$ (68% CL) for the full data combination. This limit is
close to the onset of dark energy domination at $z\sim 0.5$ Velten  _et al._
(2014). The preference for a small value for the position of $z_{1}$ is the
most pronounced for the final data combination.

The position of the second knot, $z_{2}$, is less constrained, but also moves
toward lower redshifts when including more data. Whereas the first knot,
$z_{1}$, mostly is needed to model the increase of the uncertainty of $\sum
m_{\nu}(z)$ at low redshifts, the second knot, $z_{2}$, determines the width
of the regime where the uncertainty of $\sum m_{\nu}(z)$ is the smallest.

###  IV.3 Implications for neutrino mass decay models

In this section, we note that our a priori model-independent reconstruction
can constrain specific models that predict new physics beyond the Standard
Model of Particle Physics. In particular, as discussed in Sec. II.2, it has
been proposed that neutrinos can decay into dark radiation in the late
Universe. To apply our constraints to these models, we need to carefully
distinguish two cases. First, in the scenario in which neutrinos decay while
they are still relativistic, the anisotropic stress changes before the decay
Barenboim  _et al._ (2020). This makes it hard to directly transfer our bounds
to such models, because these models change the neutrino properties both for
$z<z_{\rm nr}$ and for $z>z_{\rm nr}$. If instead neutrinos decay after they
become nonrelativistic, all cosmological neutrino properties for $z>z_{\rm
nr}$ are the same as in the standard case. Thus, any change of neutrino
properties for $z<z_{\rm nr}$ due to the decay will not affect the properties
for $z>z_{\rm nr}$. In particular, any constraints on cosmological parameters
before the decay, such as $\sum m_{\nu}$ or $\omega_{\nu}h^{2}$, will be
unaltered in the model compared to the standard $\Lambda$CDM model with
constant-mass, nondecaying neutrinos. This implies that our constraints for
$z>z_{\rm nr}$ can be directly transferred to nonrelativistic neutrino decay
models.

Our constraint is shown in our final result plot, Fig. 4, which will be
discussed in more detail in the next section. In this plot, the red dashed
line is the neutrino mass sum as a function of the redshift $z_{\rm nr}$ at
which the largest neutrino mass eigenstate $m_{\nu_{3}}$ becomes
nonrelativistic Ichikawa  _et al._ (2005),

$1+z_{\mathrm{nr}}\sim\frac{2\,m_{\nu_{3}}}{\text{meV}},$ (3)

assuming a normal neutrino mass hierarchy. The smaller neutrino mass
eigenstates will become nonrelativistic later and thus can be neglected in
this curve. Taking the intersection of the red dashed line and our result
curve, we find a limit of

$\sum m_{\nu}<0.21\text{ eV}\hskip 30.00005pt\text{(95\% CL)}$ (4)

for $z>z_{\rm nr}$. For the nonrelativistic decay scenario, we bound $\sum
m_{\nu}$ with our result for $\sum m_{\nu}(z=203)$. Note that
$103<z_{\mathrm{nr}}<203$ are the redshifts at which the individual neutrinos
would become nonrelativistic if they had a combined mass of $\sum
m_{\nu}=0.21$ eV (95% CL) (see Fig. 4).

Using Planck 2015 data, it has been previously shown that neutrino decay
models can relax the cosmological neutrino mass bound to $\sum m_{\nu}\lesssim
0.9$ eV for the non-relativistic decay scenario Escudero and Fairbairn (2019);
Chacko  _et al._ (2020a, b); Escudero  _et al._ (2020). Note that the
cosmological impact of the resulting dark radiation is negligible in these
models, as the neutrinos decay during the matter domination era. Compared to
the earlier results in Escudero and Fairbairn (2019); Chacko  _et al._ (2020a,
b); Escudero  _et al._ (2020), our bound in Eq. (4) might be tightened due to
two different effects. First, we use more recent CMB data from the latest 2018
Planck release Aghanim  _et al._ (2018). Second, the redshift of the
nonrelativistic neutrino transition decreases with decreasing neutrino mass,
such that the bound gets tighter for later neutrino decays (see Fig. 4). We
also note that the neutrino mass bound of $\sum m_{\nu}\lesssim 0.9$ eV
Escudero and Fairbairn (2019); Chacko  _et al._ (2020a, b); Escudero  _et al._
(2020) relied on simplified assumptions for the equations of motions, as noted
in Ref. Barenboim  _et al._ (2020). This might have an impact on this model-
specific mass bound but does not alter our results due to the reasons
discussed above.

The limits extracted above do not include all fundamental assumptions of the
nonrelativistic decay models, e.g., a decrease in the total neutrino mass. In
principle, it is possible that the wider limits that we find at low redshifts
could influence the limits at high redshifts. To check the robustness of our
estimates, we extended the extraction to include a theoretical prior on the
neutrino mass sum with $\sum m_{\nu,z_{i}}\leq\sum m_{\nu,z_{j}}$ for
$z_{i}\leq z_{j}$.Using this prior, we found an upper limit of $\sum
m_{\nu}<0.26$ eV (95% CL) for the nonrelativistic decay scenario. This shows
that the neutrino mass bound for the decay models is stable with respect to a
potential late-time change of the neutrino mass sum.

##  V Discussion and Conclusion

In this paper, we reconstructed the cosmological neutrino mass sum as a
function of redshift. This reconstruction was model independent, such that no
specific mass model or interaction with the dark energy or dark radiation
sector was assumed. Our final result is shown in Fig. 4. The figure shows the
reconstructed neutrino mass sum for the full data combination of CMB
temperature, polarization and lensing, combined with BAO and SN data. We
highlight the redshift when dark energy starts to dominate. We also highlight
the redshift $z_{\rm nr}$ at which the largest neutrino mass eigenstate
becomes nonrelativistic, in order to constrain models Escudero and Fairbairn
(2019); Chacko  _et al._ (2020a, b); Escudero  _et al._ (2020) that predict
neutrino decay after this nonrelativistic transition (see Sec. IV and below).

Our result is consistent with neutrino masses that are constant in time, as
predicted by the Standard Model of Cosmology. We observe a small increase of
the bound on $\sum m_{\nu}$ at high redshifts ($z\gtrsim 500$), as well as a
large increase of this bound at low redshifts coinciding with the onset of
dark energy domination. The large mass bound at low redshifts could be
explained in two different ways.

On the one hand, it could simply arise due to the strong degeneracy between
$\Omega_{\Lambda}$ and $\sum m_{\nu}$ at low redshifts, which is shown in Fig.
2. We have explicitly tested this scenario by imposing a strong prior on
$\Omega_{\Lambda}=0.69\pm 0.01$. In that case, the constraint on $\Sigma
m_{\nu}(z=0)$ decreased by 42% to $\sum m_{\nu}(z=0)=0.84$ eV, which shows an
impact from correlations with dark energy but is still larger than standard
neutrino mass bounds. This latter effect might be due to the large number of
additional free parameters needed for the reconstruction.

On the other hand, the large mass bound at low redshifts could be explained by
new physics. It has been a longstanding puzzle why the energy scales of
neutrino masses and dark energy are close
($m_{\nu}\sim\sqrt[4]{\rho_{\Lambda}}\sim{\rm meV}$) but far away from all
other known fundamental energy scales, such as the Higgs or Planck scales. If
cosmological data allow for larger neutrino masses in the late Universe,
$0\leq z\leq 1$ (such as suggested in Refs. Brookfield  _et al._ (2006a, b);
Koksbang and Hannestad (2017); Lorenz  _et al._ (2019); Muir  _et al._ (2020);
Battye and Moss (2014); Wyman  _et al._ (2014); Beutler  _et al._ (2014);
Poulin  _et al._ (2018)), an intriguing theoretical connection to the dark
energy sector could appear Fardon  _et al._ (2004); Bjaelde  _et al._ (2008);
Ayaita  _et al._ (2016); Mandal  _et al._ (2019); D’Amico  _et al._ (2019);
Dvali and Funcke (2016a). In particular, our results agree with the insights
from previous work presented in Ref. Lorenz  _et al._ (2019). In that study, a
specific model was assumed that predicted the late neutrino mass generation in
the Universe arising from a supercooled phase transition. We note that our
results also show the same trend recently obtained in Ref. Muir  _et al._
(2020) by the Dark Energy Survey collaboration, where a higher neutrino mass
and a corresponding low $\Omega_{\Lambda}$ was found.

The analysis presented in this paper could be extended in different ways. On
the theory side, it would be interesting to reconstruct $\sum m_{\nu}(z)$
fully implementing potential interactions with the dark energy or the dark
radiation sector (see below). In particular, a coupling between the dark
energy and the neutrino sectors would affect the dark energy perturbations,
leading to more features in the model and additional ways to decrease the
degeneracy between dark energy and massive neutrinos.

Our approach is conservative in the modeling of the neutrino mass but not
necessarily in the constraints that we obtain: we choose the simplest
parametrization of the mass, without any model-specific assumptions, but the
constraints could be weakened or strengthened when model-specific physics is
introduced. For example, the neutrino mass bounds could become much weaker for
modified gravity models Bellomo  _et al._ (2017); Hagstotz  _et al._ (2019a,
b) or when introducing a nontrivial coupling between the neutrino and dark
energy sectors, for example in the model considered in Ref. Lorenz  _et al._
(2019). On the other hand, the bounds are expected to become stronger in
neutrino decay models when taking into account neutrino perturbations Escudero
and Fairbairn (2019); Chacko  _et al._ (2020a, b); Escudero  _et al._ (2020).
Such features could be captured by a more general reconstruction of the
neutrino mass including nontrivial couplings to different sectors.

On the methodology side, one could model the number of knots or bins with a
hyperparameter Dimatteo and Kass (2001) in order to allow the data to decide
on the number of knots. To speed up computations, this can be fitted with
empirical Bayes, where the hyperparameters are chosen as the maximizers of the
marginalized likelihood Seikel  _et al._ (2012b); Kern  _et al._ (2017). The
exact number of knots or bins can impact the results by under- or overfitting
the reconstructed function.

On the data side, it would be interesting to reconstruct $\sum m_{\nu}$
including more large scale structure data, in particular galaxy clustering,
cosmic shear, and cross-correlations of either the $\ell$ISW effect and galaxy
clustering or of CMB lensing and tracers of the large scale structure. These
additional measurements of the suppression of the matter power spectrum would
also be an additional way to break the degeneracy with dark energy (see, e.g.
Lahav  _et al._ (2010); Mishra-Sharma  _et al._ (2018); Brinckmann  _et al._
(2019); Boyle (2019); Boyle  _et al._ (2020); Pearson and Zahn (2014);
Giusarma  _et al._ (2018); Yu  _et al._ (2018); Lesgourgues  _et al._ (2008);
Xu (2016); Bayer  _et al._ (2021)). To fully exploit these data will however
require careful modeling of the nonlinear scales (see, e.g. Saito  _et al._
(2008); Brandbyge and Hannestad (2010); Lesgourgues  _et al._ (2009); Wong
(2011); Bird  _et al._ (2012); Adamek  _et al._ (2017); Villaescusa-Navarro
_et al._ (2020); Mead  _et al._ (2020); Nascimento and Loverde (2021); Bose
_et al._ (2021)), and of potential systematics, such as the scale-dependent
galaxy bias in the presence of massive neutrinos Castorina  _et al._ (2015);
Raccanelli  _et al._ (2019); Vagnozzi  _et al._ (2018b); Chiang  _et al._
(2019); Xu  _et al._ (2020).

Finally, we comment on the implications of our results for particle physics
models and experiments. In particular, our study strengthens the previously
reported cosmological mass bound $\sum m_{\nu}<0.9$ eV (95% CL) of
nonrelativistic neutrino mass decay scenarios Escudero and Fairbairn (2019);
Chacko  _et al._ (2020a, b); Escudero  _et al._ (2020) to $\sum m_{\nu}<0.21$
eV (95% CL). We note that the constraints are still weaker than for constant
neutrino masses in time, because decay scenarios are insensitive to late-time
cosmological data from BAO and SN. However, our study pushes the neutrino mass
bounds of such models below the sensitivity of the KATRIN experiment. This
implies that a neutrino mass discovery at KATRIN would hint toward other
models, such as models predicting post-recombination neutrino mass generation
and subsequent relic neutrino annihilation, such as proposed in Ref. Dvali and
Funcke (2016a). The cosmological disappearance of the neutrino mass parameter
might also be explained in the context of modified gravity Bellomo  _et al._
(2017); Hagstotz  _et al._ (2019a, b), but this degeneracy between massive
neutrino and modify gravity effects will be broken by future surveys, such as
Euclid Hagstotz  _et al._ (2019b).

## Acknowledgments

We thank Alexandre Réfrégier, Thejs Brinckmann, Sunny Vagnozzi, Pedro Machado,
Roni Harnik, Raphaël Sgier, Jessie Muir, Yvonne Wong and Neal Dalal for
helpful discussions and/or comments on the draft. We also acknowledge helpful
discussions at the Cosmology from Home Conference in 2020. CSL acknowledges
financial support by Alexandre Réfrégier’s Cosmology research group at ETH
Zürich. Research at Perimeter Institute is supported in part by the Government
of Canada through the Department of Innovation, Science and Industry Canada
and by the Province of Ontario through the Ministry of Colleges and
Universities. EC acknowledges support from the STFC Ernest Rutherford
Fellowship ST/M004856/2 and STFC Consolidated Grant ST/S00033X/1, and from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant agreement No. 849169). ML has been
funded in part by ETH Foundations of Data Science (ETH-FDS). Parts of this
research was conducted using the computer cluster Euler at ETH Zürich.

## References

  * Group  _et al._ (2020) Particle Data Group, P. A. Zyla, _et al._ , “Review of Particle Physics,” Prog. Theor. Exp. Phys. 2020 (2020).
  * Aker  _et al._ (2021a) M. Aker _et al._ , “First direct neutrino-mass measurement with sub-eV sensitivity,” (2021a), arXiv:2105.08533 [hep-ex] .
  * Aker  _et al._ (2021b) M. Aker _et al._ , “Analysis methods for the first KATRIN neutrino-mass measurement,” (2021b), arXiv:2101.05253 [hep-ex] .
  * Aghanim  _et al._ (2018) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” (2018), arXiv:1807.06209 [astro-ph.CO] .
  * Beacom  _et al._ (2004) J. F. Beacom, N. F. Bell, and S. Dodelson, “Neutrinoless universe,” Phys. Rev. Lett. 93, 121302 (2004), arXiv:astro-ph/0404585 [astro-ph] .
  * Hannestad (2005) S. Hannestad, “Structure formation with strongly interacting neutrinos - Implications for the cosmological neutrino mass bound,” J. Cosmol. Astropart. Phys. 0502, 11 (2005), arXiv:astro-ph/0411475 [astro-ph] .
  * Lancaster  _et al._ (2017) L. Lancaster, F.-Y. Cyr-Racine, L. Knox, and Z. Pan, “A tale of two modes: Neutrino free-streaming in the early universe,” JCAP 1707, 033 (2017), arXiv:1704.06657 [astro-ph.CO] .
  * Follin  _et al._ (2015) B. Follin, L. Knox, M. Millea, and Z. Pan, “First Detection of the Acoustic Oscillation Phase Shift Expected from the Cosmic Neutrino Background,” Phys. Rev. Lett. 115, 091301 (2015), arXiv:1503.07863 [astro-ph.CO] .
  * Dvali and Funcke (2016a) G. Dvali and L. Funcke, “Small neutrino masses from gravitational $\theta$-term,” Phys. Rev. D 93, 113002 (2016a), arXiv:1602.03191 [hep-ph] .
  * Bellomo  _et al._ (2017) N. Bellomo, E. Bellini, B. Hu, R. Jimenez, C. Pena-Garay, and L. Verde, “Hiding neutrino mass in modified gravity cosmologies,” JCAP 02, 043 (2017), arXiv:1612.02598 [astro-ph.CO] .
  * Hagstotz  _et al._ (2019a) S. Hagstotz, M. Gronke, D. Mota, and M. Baldi, “Breaking cosmic degeneracies: Disentangling neutrinos and modified gravity with kinematic information,” Astron. Astrophys. 629, A46 (2019a), arXiv:1902.01868 [astro-ph.CO] .
  * Hagstotz  _et al._ (2019b) S. Hagstotz, M. Costanzi, M. Baldi, and J. Weller, “Joint halo-mass function for modified gravity and massive neutrinos – I. Simulations and cosmological forecasts,” Monthly Notices of the Royal Astronomical Society 486, 3927–3941 (2019b).
  * Esteban and Salvado (2021) I. Esteban and J. Salvado, “Long range interactions in cosmology: Implications for neutrinos,” (2021), arXiv:2101.05804 [hep-ph] .
  * Joudaki and Kaplinghat (2012) S. Joudaki and M. Kaplinghat, “Dark energy and neutrino masses from future measurements of the expansion history and growth of structure,” Phys. Rev. D 86, 023526 (2012), arXiv:1106.0299 [astro-ph.CO] .
  * Joudaki (2013) S. Joudaki, “Constraints on neutrino mass and light degrees of freedom in extended cosmological parameter spaces,” Phys. Rev. D 87, 083523 (2013).
  * Allison  _et al._ (2015) R. Allison, P. Caucal, E. Calabrese, J. Dunkley, and T. Louis, “Towards a cosmological neutrino mass detection,” Phys. Rev. D 92, 123535 (2015), arXiv:1509.07471 [astro-ph.CO] .
  * Lorenz  _et al._ (2017) C. S. Lorenz, E. Calabrese, and D. Alonso, “Distinguishing between Neutrinos and time-varying Dark Energy through Cosmic Time,” Phys. Rev. D 96, 043510 (2017), arXiv:1706.00730 [astro-ph.CO] .
  * Yang  _et al._ (2017) Weiqiang Yang, Rafael C. Nunes, Supriya Pan, and David F. Mota, “Effects of neutrino mass hierarchies on dynamical dark energy models,” Phys. Rev. D 95, 103522 (2017), arXiv:1703.02556 [astro-ph.CO] .
  * Vagnozzi  _et al._ (2018a) S. Vagnozzi, S. Dhawan, M. Gerbino, K. Freese, A. Goobar, and O. Mena, “Constraints on the sum of the neutrino masses in dynamical dark energy models with $w(z)\geq-1$ are tighter than those obtained in $\Lambda$CDM,” Phys. Rev. D 98, 083501 (2018a), arXiv:1801.08553 [astro-ph.CO] .
  * Roy Choudhury and Choubey (2018) S. Roy Choudhury and S. Choubey, “Updated Bounds on Sum of Neutrino Masses in Various Cosmological Scenarios,” JCAP 09, 017 (2018), arXiv:1806.10832 [astro-ph.CO] .
  * Roy Choudhury and Hannestad (2019) S. Roy Choudhury and S. Hannestad, “Updated results on neutrino mass and mass hierarchy from cosmology with Planck 2018 likelihoods,” (2019), arXiv:1907.12598 [astro-ph.CO] .
  * Oldengott  _et al._ (2019) I. M. Oldengott, G. Barenboim, S. Kahlen, J. Salvado, and D. J. Schwarz, “How to relax the cosmological neutrino mass bound,” JCAP 04, 049 (2019), arXiv:1901.04352 [astro-ph.CO] .
  * Escudero and Fairbairn (2019) M. Escudero and M. Fairbairn, “Cosmological constraints on invisible neutrino decays revisited,” Phys. Rev. D 100, 103531 (2019).
  * Chacko  _et al._ (2020a) Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, “Cosmological Limits on the Neutrino Mass and Lifetime,” JHEP 04, 020 (2020a), arXiv:1909.05275 [hep-ph] .
  * Chacko  _et al._ (2020b) Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, “Determining the Neutrino Lifetime from Cosmology,” (2020b), arXiv:2002.08401 [astro-ph.CO] .
  * Escudero  _et al._ (2020) M. Escudero, J. Lopez-Pavon, N. Rius, and S. Sandner, “Relaxing Cosmological Neutrino Mass Bounds with Unstable Neutrinos,” (2020), arXiv:2007.04994 [hep-ph] .
  * Barenboim  _et al._ (2020) Ga. Barenboim, J. Z. Chen, S. Hannestad, I. M. Oldengott, T. Tram, and Y. Y. Y. Wong, “Invisible neutrino decay in precision cosmology,” (2020), arXiv:2011.01502 [astro-ph.CO] .
  * Fardon  _et al._ (2004) R. Fardon, A. E. Nelson, and N. Weiner, “Dark energy from mass varying neutrinos,” JCAP 10, 005 (2004), arXiv:astro-ph/0309800 .
  * Koksbang and Hannestad (2017) S. M. Koksbang and S. Hannestad, “Constraining dynamical neutrino mass generation with cosmological data,” JCAP 09, 014 (2017), arXiv:1707.02579 [astro-ph.CO] .
  * Lorenz  _et al._ (2019) C. S. Lorenz, L. Funcke, E. Calabrese, and S. Hannestad, “Time-varying neutrino mass from a supercooled phase transition: current cosmological constraints and impact on the $\Omega_{m}$-$\sigma_{8}$ plane,” Phys. Rev. D99, 023501 (2019), arXiv:1811.01991 [astro-ph.CO] .
  * Brookfield  _et al._ (2006a) A. W. Brookfield, C. van de Bruck, D. F. Mota, and D. Tocchini-Valentini, “Cosmology with massive neutrinos coupled to dark energy,” Physical Review Letters 96 (2006a), 10.1103/physrevlett.96.061301.
  * Brookfield  _et al._ (2006b) A. W. Brookfield, C. van de Bruck, D. F. Mota, and D. Tocchini-Valentini, “Cosmology of mass-varying neutrinos driven by quintessence: Theory and observations,” Phys. Rev. D 73, 083515 (2006b).
  * Battye and Moss (2014) R. A. Battye and A. Moss, “Evidence for Massive Neutrinos from Cosmic Microwave Background and Lensing Observations,” Phys. Rev. Lett. 112, 051303 (2014), arXiv:1308.5870 [astro-ph.CO] .
  * Wyman  _et al._ (2014) M. Wyman, D. H. Rudd, R. A. Vanderveld, and W. Hu, “Neutrinos Help Reconcile Planck Measurements with the Local Universe,” Phys. Rev. Lett. 112, 051302 (2014), arXiv:1307.7715 [astro-ph.CO] .
  * Beutler  _et al._ (2014) F. Beutler _et al._ (BOSS), “The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: signs of neutrino mass in current cosmological data sets,” Mon. Not. Roy. Astron. Soc. 444, 3501–3516 (2014), arXiv:1403.4599 [astro-ph.CO] .
  * Poulin  _et al._ (2018) V. Poulin, K. K. Boddy, S. Bird, and M. Kamionkowski, “Implications of an extended dark energy cosmology with massive neutrinos for cosmological tensions,” Phys. Rev. D 97, 123504 (2018), arXiv:1803.02474 [astro-ph.CO] .
  * Bjaelde  _et al._ (2008) O. E. Bjaelde, A. W. Brookfield, C. van de Bruck, S. Hannestad, D. F. Mota, L. Schrempp, and D. Tocchini-Valentini, “Neutrino Dark Energy – Revisiting the Stability Issue,” JCAP 01, 026 (2008), arXiv:0705.2018 [astro-ph] .
  * Ayaita  _et al._ (2016) Y. Ayaita, M. Baldi, F. Führer, E. Puchwein, and C. Wetterich, “Nonlinear growing neutrino cosmology,” Phys. Rev. D 93, 063511 (2016), arXiv:1407.8414 [astro-ph.CO] .
  * Mandal  _et al._ (2019) S. Mandal, G. Y. Chitov, O. Avsajanishvili, B. Singha, and T. Kahniashvili, “Mass varying neutrinos with different quintessence potentials,” (2019), arXiv:1911.06099 [hep-ph] .
  * D’Amico  _et al._ (2019) G. D’Amico, T. Hamill, and N. Kaloper, “Neutrino Masses from Outer Space,” Phys. Lett. B 797, 134846 (2019), arXiv:1804.01542 [hep-ph] .
  * Drexlin  _et al._ (2013) G. Drexlin, V. Hannen, S. Mertens, and C. Weinheimer, “Current direct neutrino mass experiments,” Adv. High Energy Phys. 2013, 293986 (2013), arXiv:1307.0101 [physics.ins-det] .
  * Abazajian  _et al._ (2015) K. N. Abazajian _et al._ (Topical Conveners: K.N. Abazajian, J.E. Carlstrom, A.T. Lee), “Neutrino Physics from the Cosmic Microwave Background and Large Scale Structure,” Astropart. Phys. 63, 66–80 (2015), arXiv:1309.5383 [astro-ph.CO] .
  * Mishra-Sharma  _et al._ (2018) S. Mishra-Sharma, D. Alonso, and J. Dunkley, “Neutrino masses and beyond- $\Lambda$CDM cosmology with LSST and future CMB experiments,” Phys. Rev. D 97, 123544 (2018), arXiv:1803.07561 [astro-ph.CO] .
  * Holsclaw  _et al._ (2010) T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib, and D. Higdon, “Nonparametric Dark Energy Reconstruction from Supernova Data,” Phys. Rev. Letters 105, 241302 (2010), arXiv:1011.3079 [astro-ph.CO] .
  * Holsclaw  _et al._ (2011) T. Holsclaw, U. Alam, B. Sansó, H. Lee, K. Heitmann, S. Habib, and D. Higdon, “Nonparametric reconstruction of the dark energy equation of state from diverse data sets,” Phys. Rev. D 84, 083501 (2011), arXiv:1104.2041 [astro-ph.CO] .
  * Seikel  _et al._ (2012a) M. Seikel, C. Clarkson, and M. Smith, “Reconstruction of dark energy and expansion dynamics using Gaussian processes,” JCAP 2012, 036 (2012a), arXiv:1204.2832 [astro-ph.CO] .
  * Seikel and Clarkson (2013) M. Seikel and C. Clarkson, “Optimising Gaussian processes for reconstructing dark energy dynamics from supernovae,” (2013), arXiv:1311.6678 [astro-ph.CO] .
  * Alberto Vazquez  _et al._ (2012) J. Alberto Vazquez, M. Bridges, M. P. Hobson, and A. N. Lasenby, “Reconstruction of the Dark Energy equation of state,” JCAP 1209, 020 (2012), arXiv:1205.0847 [astro-ph.CO] .
  * Montiel  _et al._ (2014) A. Montiel, R. Lazkoz, I. Sendra, C. Escamilla-Rivera, and V. Salzano, “Nonparametric reconstruction of the cosmic expansion with local regression smoothing and simulation extrapolation,” Phys. Rev. D89, 043007 (2014), arXiv:1401.4188 [astro-ph.CO] .
  * Hee  _et al._ (2016) S. Hee, W. Handley, M. P. Hobson, and A. N. Lasenby, “Bayesian model selection without evidences: application to the dark energy equation-of-state,” Mon. Not. Roy. Astron. Soc. 455, 2461–2473 (2016), arXiv:1506.09024 [astro-ph.CO] .
  * Hee  _et al._ (2017) S. Hee, J. A. Vázquez, W. J. Handley, M. P. Hobson, and A. N. Lasenby, “Constraining the dark energy equation of state using Bayes theorem and the Kullback–Leibler divergence,” Mon. Not. Roy. Astron. Soc. 466, 369–377 (2017), arXiv:1607.00270 [astro-ph.CO] .
  * Gerardi  _et al._ (2019) F. Gerardi, M. Martinelli, and A. Silvestri, “Reconstruction of the Dark Energy equation of state from latest data: the impact of theoretical priors,” JCAP 1907, 042 (2019), arXiv:1902.09423 [astro-ph.CO] .
  * Joudaki  _et al._ (2018) S. Joudaki, M. Kaplinghat, R. Keeley, and D. Kirkby, “Model independent inference of the expansion history and implications for the growth of structure,” Phys. Rev. D 97, 123501 (2018), arXiv:1710.04236 [astro-ph.CO] .
  * Keeley  _et al._ (2019) R. E. Keeley, S. Joudaki, M. Kaplinghat, and D. Kirkby, “Implications of a transition in the dark energy equation of state for the $H_{0}$ and $\sigma_{8}$ tensions,” JCAP 12, 035 (2019), arXiv:1905.10198 [astro-ph.CO] .
  * Park  _et al._ (2021) M. Park, M. Raveri, and B. Jain, “Reconstructing Quintessence,” (2021), arXiv:2101.04666 [astro-ph.CO] .
  * Bonilla  _et al._ (2021) Alexander Bonilla, Suresh Kumar, Rafael C. Nunes, and Supriya Pan, “Reconstruction of the dark sectors’ interaction: A model-independent inference and forecast from GW standard sirens,” (2021), arXiv:2102.06149 [astro-ph.CO] .
  * Shafieloo  _et al._ (2012) A. Shafieloo, A. G. Kim, and E. V. Linder, “Gaussian process cosmography,” Phys. Rev. D 85, 123530 (2012), arXiv:1204.2272 [astro-ph.CO] .
  * Seikel  _et al._ (2012b) M. Seikel, S. Yahya, R. Maartens, and C. Clarkson, “Using H(z) data as a probe of the concordance model,” Phys. Rev. D 86, 083001 (2012b), arXiv:1205.3431 [astro-ph.CO] .
  * Bernal  _et al._ (2016) J. L. Bernal, L. Verde, and A. G. Riess, “The trouble with $H_{0}$,” JCAP 10, 019 (2016), arXiv:1607.05617 [astro-ph.CO] .
  * Gómez-Valent and Amendola (2018) A. Gómez-Valent and L. Amendola, “$H_{0}$ from cosmic chronometers and Type Ia supernovae, with Gaussian Processes and the novel Weighted Polynomial Regression method,” JCAP 04, 051 (2018), arXiv:1802.01505 [astro-ph.CO] .
  * Renzi and Silvestri (2020) F. Renzi and A. Silvestri, “A look at the Hubble speed from first principles,” (2020), arXiv:2011.10559 [astro-ph.CO] .
  * Bernal  _et al._ (2021) José Luis Bernal, Licia Verde, Raul Jimenez, Marc Kamionkowski, David Valcin, and Benjamin D. Wandelt, “The trouble beyond $H_{0}$ and the new cosmic triangles,” (2021), arXiv:2102.05066 [astro-ph.CO] .
  * Bridle  _et al._ (2003) S. L. Bridle, A. M. Lewis, J. Weller, and G. Efstathiou, “Reconstructing the primordial power spectrum,” Mon. Not. Roy. Astron. Soc. 342, L72 (2003), arXiv:astro-ph/0302306 .
  * Guo  _et al._ (2011) Zong-Kuan Guo, Dominik J. Schwarz, and Yuan-Zhong Zhang, “Reconstruction of the primordial power spectrum from CMB data,” JCAP 08, 031 (2011), arXiv:1105.5916 [astro-ph.CO] .
  * Hazra  _et al._ (2014) D. K. Hazra, A. Shafieloo, and T. Souradeep, “Primordial power spectrum from Planck,” JCAP 11, 011 (2014), arXiv:1406.4827 [astro-ph.CO] .
  * Hunt and Sarkar (2014) P. Hunt and S. Sarkar, “Reconstruction of the primordial power spectrum of curvature perturbations using multiple data sets,” JCAP 01, 025 (2014), arXiv:1308.2317 [astro-ph.CO] .
  * Lévy  _et al._ (2020) B. Lévy, R. Mohayaee, and S. von Hausegger, “A fast semi-discrete optimal transport algorithm for a unique reconstruction of the early Universe,” (2020), arXiv:2012.09074 [astro-ph.CO] .
  * Garcia-Quintero  _et al._ (2020) C. Garcia-Quintero, M. Ishak, and O. Ning, “Current constraints on deviations from General Relativity using binning in redshift and scale,” (2020), arXiv:2010.12519 [astro-ph.CO] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A.G. Riess, “Tensions between the Early and the Late Universe,” (2019) arXiv:1907.10625 [astro-ph.CO] .
  * Kunz (2009) M. Kunz, “The dark degeneracy: On the number and nature of dark components,” Phys. Rev. D 80, 123001 (2009), arXiv:astro-ph/0702615 .
  * Amendola  _et al._ (2013) L. Amendola, M. Kunz, M. Motta, I. D. Saltas, and I. Sawicki, “Observables and unobservables in dark energy cosmologies,” Phys. Rev. D 87, 023501 (2013), arXiv:1210.0439 [astro-ph.CO] .
  * Busti and Clarkson (2016) V. C. Busti and C. Clarkson, “Dodging the dark matter degeneracy while determining the dynamics of dark energy,” JCAP 1605, 008 (2016), arXiv:1505.01821 [astro-ph.CO] .
  * Ichikawa  _et al._ (2005) K. Ichikawa, M. Fukugita, and M. Kawasaki, “Constraining neutrino masses by CMB experiments alone,” Phys. Rev. D 71, 043001 (2005), arXiv:astro-ph/0409768 .
  * Archidiacono  _et al._ (2017) M. Archidiacono, T. Brinckmann, J. Lesgourgues, and V. Poulin, “Physical effects involved in the measurements of neutrino masses with future cosmological data,” JCAP 02, 052 (2017), arXiv:1610.09852 [astro-ph.CO] .
  * Lesgourgues  _et al._ (2013) J. Lesgourgues, G. Mangano, G. Miele, and S. Pastor,  _Neutrino Cosmology_ (Cambridge University Press, 2013).
  * Vagnozzi  _et al._ (2017) S. Vagnozzi, E. Giusarma, O. Mena, K. Freese, M. Gerbino, S. Ho, and M. Lattanzi, “Unveiling $\nu$ secrets with cosmological data: neutrino masses and mass hierarchy,” Phys. Rev. D 96, 123503 (2017), arXiv:1701.08172 [astro-ph.CO] .
  * Howlett  _et al._ (2012) C. Howlett, A. Lewis, A. Hall, and A. Challinor, “CMB power spectrum parameter degeneracies in the era of precision cosmology,” JCAP 2012, 027 (2012), arXiv:1201.3654 [astro-ph.CO] .
  * Sutherland (2018) W. Sutherland, “The CMB neutrino mass/vacuum energy degeneracy: a simple derivation of the degeneracy slopes,” Mon. Not. Roy. Astron. Soc. 477, 1913–1920 (2018), arXiv:1803.02298 [astro-ph.CO] .
  * Moresco  _et al._ (2012) M. Moresco, L. Verde, L. Pozzetti, R. Jimenez, and A. Cimatti, “New constraints on cosmological parameters and neutrino properties using the expansion rate of the Universe to z~1.75,” JCAP 07, 053 (2012), arXiv:1201.6658 [astro-ph.CO] .
  * Hall and Challinor (2012) A. C. Hall and A. Challinor, “Probing the neutrino mass hierarchy with cosmic microwave background weak lensing,” MNRAS 425, 1170–1184 (2012), arXiv:1205.6172 [astro-ph.CO] .
  * Cabass  _et al._ (2015) G. Cabass, M. Gerbino, E. Giusarma, A. Melchiorri, L. Pagano, and L. Salvati, “Constraints on the early and late integrated Sachs-Wolfe effects from the Planck 2015 cosmic microwave background anisotropies in the angular power spectra,” Phys. Rev. D 92, 063534 (2015), arXiv:1507.07586 [astro-ph.CO] .
  * Hou  _et al._ (2014) Z. Hou, C. L. Reichardt, K. T. Story, B. Follin, Keisler, _et al._ , “Constraints on Cosmology from the Cosmic Microwave Background Power Spectrum of the 2500 deg2 SPT-SZ Survey,” Astrophysical Journal 782, 74 (2014), arXiv:1212.6267 [astro-ph.CO] .
  * Lesgourgues  _et al._ (2008) J. Lesgourgues, W. Valkenburg, and E. Gaztanaga, “Constraining neutrino masses with the ISW-galaxy correlation function,” Phys. Rev. D 77, 063505 (2008), arXiv:0710.5525 [astro-ph] .
  * Lewis and Challinor (2006) A. Lewis and A. Challinor, “Weak gravitational lensing of the CMB,” Phys. Rept. 429, 1–65 (2006), arXiv:astro-ph/0601594 .
  * Lesgourgues and Pastor (2006) J. Lesgourgues and S. Pastor, “Massive neutrinos and cosmology,” Phys. Rept. 429, 307–379 (2006), arXiv:astro-ph/0603494 .
  * Hu  _et al._ (1998) W. Hu, D. J. Eisenstein, and M. Tegmark, “Weighing neutrinos with galaxy surveys,” Phys. Rev. Lett. 80, 5255–5258 (1998), arXiv:astro-ph/9712057 .
  * Henning  _et al._ (2018) J.W. Henning _et al._ (SPT), “Measurements of the Temperature and E-Mode Polarization of the CMB from 500 Square Degrees of SPTpol Data,” Astrophys. J. 852, 97 (2018), arXiv:1707.09353 [astro-ph.CO] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), “Planck 2018 results. V. CMB power spectra and likelihoods,” Astron. Astrophys. 641, A5 (2020a), arXiv:1907.12875 [astro-ph.CO] .
  * Choi  _et al._ (2020) S. K. Choi _et al._ (ACT), “The Atacama Cosmology Telescope: A Measurement of the Cosmic Microwave Background Power Spectra at 98 and 150 GHz,” (2020), arXiv:2007.07289 [astro-ph.CO] .
  * Bianchini  _et al._ (2020) F. Bianchini _et al._ (SPT), “Constraints on Cosmological Parameters from the 500 deg2 SPTpol Lensing Power Spectrum,” Astrophys. J. 888, 119 (2020), arXiv:1910.07157 [astro-ph.CO] .
  * Darwish  _et al._ (2020) O. Darwish _et al._ , “The Atacama Cosmology Telescope: A CMB lensing mass map over 2100 square degrees of sky and its cross-correlation with BOSS-CMASS galaxies,”  (2020), 10.1093/mnras/staa3438, arXiv:2004.01139 [astro-ph.CO] .
  * Calabrese  _et al._ (2017) E. Calabrese, D. Alonso, and J. Dunkley, “Complementing the ground-based CMB-S4 experiment on large scales with the PIXIE satellite,” Phys. Rev. D 95, 063504 (2017), arXiv:1611.10269 [astro-ph.CO] .
  * Eisenstein  _et al._ (2005) D. J. Eisenstein _et al._ (SDSS), “Detection of the Baryon Acoustic Peak in the Large-Scale Correlation Function of SDSS Luminous Red Galaxies,” Astrophys. J. 633, 560–574 (2005), arXiv:astro-ph/0501171 .
  * Cole  _et al._ (2005) S. Cole _et al._ (2dFGRS), “The 2dF Galaxy Redshift Survey: Power-spectrum analysis of the final dataset and cosmological implications,” Mon. Not. Roy. Astron. Soc. 362, 505–534 (2005), arXiv:astro-ph/0501174 .
  * Weinberg  _et al._ (2013) D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess, and E. Rozo, “Observational probes of cosmic acceleration,” Physics Reports 530, 87–255 (2013), arXiv:1201.2434 [astro-ph.CO] .
  * Bassett and Hlozek (2009) B. A. Bassett and R. Hlozek, “Baryon Acoustic Oscillations,” (2009), arXiv:0910.5224 [astro-ph.CO] .
  * Alam  _et al._ (2020) S. Alam _et al._ (eBOSS), “The Completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Cosmological Implications from two Decades of Spectroscopic Surveys at the Apache Point observatory,” (2020), arXiv:2007.08991 [astro-ph.CO] .
  * de Carvalho  _et al._ (2018) E. de Carvalho, A. Bernui, G.C. Carvalho, C.P. Novaes, and H.S. Xavier, “Angular Baryon Acoustic Oscillation measure at $z=2.225$ from the SDSS quasar survey,” JCAP 04, 064 (2018), arXiv:1709.00113 [astro-ph.CO] .
  * Ata  _et al._ (2018) M. Ata _et al._ , “The clustering of the SDSS-IV extended Baryon Oscillation Spectroscopic Survey DR14 quasar sample: first measurement of baryon acoustic oscillations between redshift 0.8 and 2.2,” Mon. Not. Roy. Astron. Soc. 473, 4773–4794 (2018), arXiv:1705.06373 [astro-ph.CO] .
  * Hou  _et al._ (2021) J. Hou, A. G. Sánchez, A. J. Ross, A. Smith, Neveux, _et al._ , “The completed SDSS-IV extended Baryon Oscillation Spectroscopic Survey: BAO and RSD measurements from anisotropic clustering analysis of the quasar sample in configuration space between redshift 0.8 and 2.2,” Monthly Notices of the Royal Astronomical Society 500, 1201–1221 (2021), arXiv:2007.08998 [astro-ph.CO] .
  * McDonald and Eisenstein (2007) P. McDonald and D. Eisenstein, “Dark energy and curvature from a future baryonic acoustic oscillation survey using the Lyman-alpha forest,” Phys. Rev. D 76, 063009 (2007), arXiv:astro-ph/0607122 .
  * Bautista  _et al._ (2017) J. E. Bautista _et al._ , “Measurement of baryon acoustic oscillation correlations at $z=2.3$ with SDSS DR12 Ly$\alpha$-Forests,” Astron. Astrophys. 603, A12 (2017), arXiv:1702.00176 [astro-ph.CO] .
  * Seljak  _et al._ (2005) U. Seljak _et al._ (SDSS), “Cosmological parameter analysis including SDSS Ly-alpha forest and galaxy bias: Constraints on the primordial spectrum of fluctuations, neutrino mass, and dark energy,” Phys. Rev. D 71, 103515 (2005), arXiv:astro-ph/0407372 .
  * Seljak  _et al._ (2006) U. Seljak, A. Slosar, and P. McDonald, “Cosmological parameters from combining the Lyman-alpha forest with CMB, galaxy clustering and SN constraints,” JCAP 10, 014 (2006), arXiv:astro-ph/0604335 .
  * Palanque-Delabrouille  _et al._ (2015a) N. Palanque-Delabrouille _et al._ , “Constraint on neutrino masses from SDSS-III/BOSS Ly$\alpha$ forest and other cosmological probes,” JCAP 02, 045 (2015a), arXiv:1410.7244 [astro-ph.CO] .
  * Palanque-Delabrouille  _et al._ (2015b) N. Palanque-Delabrouille _et al._ , “Neutrino masses and cosmology with Lyman-alpha forest power spectrum,” JCAP 11, 011 (2015b), arXiv:1506.05976 [astro-ph.CO] .
  * Palanque-Delabrouille  _et al._ (2020) N. Palanque-Delabrouille, C. Yèche, N. Schöneberg, J. Lesgourgues, M. Walther, S. Chabanier, and E. Armengaud, “Hints, neutrino bounds and WDM constraints from SDSS DR14 Lyman-$\alpha$ and Planck full-survey data,” JCAP 04, 038 (2020), arXiv:1911.09073 [astro-ph.CO] .
  * Croft  _et al._ (1998) R. A. C. Croft, D. H. Weinberg, N. Katz, and L. Hernquist, “Recovery of the power spectrum of mass fluctuations from observations of the ly$\alpha$ forest,” The Astrophysical Journal 495, 44–62 (1998).
  * Zhan  _et al._ (2005) H. Zhan, R. Dave, D. Eisenstein, and N. Katz, “Ly-alpha flux power spectrum and its covariance,” Mon. Not. Roy. Astron. Soc. 363, 1145–1154 (2005), arXiv:astro-ph/0504419 .
  * Pedersen  _et al._ (2020) C. Pedersen, A. Font-Ribera, T. D. Kitching, P. McDonald, S. Bird, A. Slosar, K. K. Rogers, and A. Pontzen, “Massive neutrinos and degeneracies in Lyman-alpha forest simulations,” JCAP 04, 025 (2020), arXiv:1911.09596 [astro-ph.CO] .
  * Hui and Gnedin (1997) L. Hui and N. Y. Gnedin, “Equation of state of the photoionized intergalactic medium,” Mon. Not. Roy. Astron. Soc. 292, 27 (1997), arXiv:astro-ph/9612232 .
  * Walther  _et al._ (2020) M. Walther, E. Armengaud, C. Ravoux, N. Palanque-Delabrouille, C. Yèche, and Z. Lukić, “Simulating intergalactic gas for DESI-like small scale Lyman\alpha forest observations,” (2020), arXiv:2012.04008 [astro-ph.CO] .
  * Goobar and Leibundgut (2011) A. Goobar and B. Leibundgut, “Supernova Cosmology: Legacy and Future,” Annual Review of Nuclear and Particle Science 61, 251–279 (2011), arXiv:1102.1431 [astro-ph.CO] .
  * Boyle (2019) A. Boyle, “Understanding the neutrino mass constraints achievable by combining CMB lensing and spectroscopic galaxy surveys,” JCAP 04, 038 (2019), [Erratum: JCAP 10, E01 (2019)], arXiv:1811.07636 [astro-ph.CO] .
  * Minkowski (1977) P. Minkowski, “$\mu\rightarrow$e$\gamma$ at a rate of one out of 109 muon decays?” Phys. Lett. B67, 421 – 428 (1977).
  * M. Gell-Mann, P. Ramond, and R. Slansky (edited by P. Van Nieuwenhuizen and D. Freedman, Elsevier, Amsterdam, 1979) M. Gell-Mann, P. Ramond, and R. Slansky,  _in Proceedings of the Supergravity Stony Brook Workshop_ (edited by P. Van Nieuwenhuizen and D. Freedman, Elsevier, Amsterdam, 1979).
  * Yanagida (1980) T. Yanagida, “Horizontal symmetry and masses of neutrinos,” Prog. Theor. Phys. 64, 1103–1105 (1980).
  * Mohapatra and Senjanović (1980) R. N. Mohapatra and G. Senjanović, “Neutrino mass and spontaneous parity nonconservation,” Phys. Rev. Lett. 44, 912–915 (1980).
  * Schechter and Valle (1980) J. Schechter and J. W. F. Valle, “Neutrino masses in su(2) $\bigotimes$ u(1) theories,” Phys. Rev. D 22, 2227–2235 (1980).
  * Schechter and Valle (1982) J. Schechter and J. W. F. Valle, “Neutrino decay and spontaneous violation of lepton number,” Phys. Rev. D 25, 774–783 (1982).
  * Afshordi  _et al._ (2005) N. Afshordi, M. Zaldarriaga, and K. Kohri, “On the stability of dark energy with mass-varying neutrinos,” Phys. Rev. D 72, 065024 (2005), arXiv:astro-ph/0506663 .
  * Dvali and Funcke (2016b) G. Dvali and L. Funcke, “Domestic Axion,” (2016b), arXiv:1608.08969 [hep-ph] .
  * Funcke  _et al._ (2020) L. Funcke, G. Raffelt, and E. Vitagliano, “Distinguishing Dirac and Majorana neutrinos by their decays via Nambu-Goldstone bosons in the gravitational-anomaly model of neutrino masses,” Phys. Rev. D 101, 015025 (2020), arXiv:1905.01264 [hep-ph] .
  * Betti  _et al._ (2019) M. G. Betti, M. Biasotti, A. Boscá, F. Calle, N. Canci, G. Cavoto, C. Chang, A. G. Cocco, A. P. Colijn, J. Conrad, and et al., “Neutrino physics with the ptolemy project: active neutrino properties and the light sterile case,” Journal of Cosmology and Astroparticle Physics 2019, 047–047 (2019).
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), “Planck 2015 results. XIII. Cosmological parameters,” Astron. Astrophys. 594, A13 (2016), arXiv:1502.01589 [astro-ph.CO] .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), “Planck 2018 results. VIII. Gravitational lensing,” Astron. Astrophys. 641, A8 (2020b), arXiv:1807.06210 [astro-ph.CO] .
  * Beutler  _et al._ (2011) F. Beutler, C. Blake, M. Colless, D. H. Jones, L. Staveley-Smith, L. Campbell, Q. Parker, W. Saunders, and F. Watson, “The 6dF Galaxy Survey: baryon acoustic oscillations and the local Hubble constant,” Monthly Notices of the Royal Astronomical Society 416, 3017–3032 (2011), arXiv:1106.3366 [astro-ph.CO] .
  * Ross  _et al._ (2015) A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden, and M. Manera, “The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at $z=0.15$,” Mon. Not. Roy. Astron. Soc. 449, 835–847 (2015), arXiv:1409.3242 [astro-ph.CO] .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), “The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample,” Mon. Not. Roy. Astron. Soc. 470, 2617–2652 (2017), arXiv:1607.03155 [astro-ph.CO] .
  * Blomqvist  _et al._ (2019) M. Blomqvist _et al._ , “Baryon acoustic oscillations from the cross-correlation of Ly$\alpha$ absorption and quasars in eBOSS DR14,” Astron. Astrophys. 629, A86 (2019), arXiv:1904.03430 [astro-ph.CO] .
  * Scolnic  _et al._ (2018) D. M. Scolnic, D. O. Jones, A. Rest, Y. C. Pan, R. Chornock, R. J. Foley, M. E. Huber, R. Kessler, G. Narayan, A. G. Riess, and et al., “The complete light-curve sample of spectroscopically confirmed sne ia from pan-starrs1 and cosmological constraints from the combined pantheon sample,” The Astrophysical Journal 859, 101 (2018).
  * Schöneberg  _et al._ (2019) N. Schöneberg, J. Lesgourgues, and D. C. Hooper, “The BAO+BBN take on the Hubble tension,” JCAP 10, 029 (2019), arXiv:1907.11594 [astro-ph.CO] .
  * Lewis  _et al._ (2000) A. Lewis, A. Challinor, and A. Lasenby, “Efficient computation of CMB anisotropies in closed FRW models,” Astrophys. J. 538, 473–476 (2000), arXiv:astro-ph/9911177 .
  * Lewis and Bridle (2002) A. Lewis and S. Bridle, “Cosmological parameters from CMB and other data: A Monte Carlo approach,” Phys. Rev. D 66, 103511 (2002), arXiv:astro-ph/0205436 .
  * Wasserman (2006) L. Wasserman,  _All of Nonparametric Statistics_ (Springer-Verlag New York, 2006).
  * Carbone  _et al._ (2008) Carmelita Carbone, Volker Springel, Carlo Baccigalupi, Matthias Bartelmann, and Sabino Matarrese, “Full-sky maps for gravitational lensing of the CMB,” Mon. Not. Roy. Astron. Soc. 388, 1618 (2008), arXiv:0711.2655 [astro-ph] .
  * Muccino  _et al._ (2020) M. Muccino, L. Izzo, O. Luongo, K. Boshkayev, L. Amati, M. Della Valle, G. B. Pisani, and E. Zaninoni, “Tracing dark energy history with gamma ray bursts,” arXiv e-prints , arXiv:2012.03392 (2020), arXiv:2012.03392 [astro-ph.CO] .
  * Dimatteo and Kass (2001) Genovese C.R. Dimatteo, I. and R.E. Kass, “Bayesian curve‐fitting with free‐knot splines,” Biometrika 88, 1055–1071 (2001).
  * Vázquez  _et al._ (2012) J. Alberto Vázquez, M. Bridges, M. P. Hobson, and A. N. Lasenby, “Model selection applied to reconstruction of the Primordial Power Spectrum,” JCAP 2012, 006 (2012), arXiv:1203.1252 [astro-ph.CO] .
  * Millea and Bouchet (2018) Marius Millea and François Bouchet, “Cosmic Microwave Background Constraints in Light of Priors Over Reionization Histories,” Astron. Astrophys. 617, A96 (2018), arXiv:1804.08476 [astro-ph.CO] .
  * Lewis (2019) A. Lewis, “GetDist: a Python package for analysing Monte Carlo samples,” (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Gariazzo  _et al._ (2018) S. Gariazzo, M. Archidiacono, P. F. de Salas, O. Mena, C. A. Ternes, and M. Tórtola, “Neutrino masses and their ordering: Global Data, Priors and Models,” JCAP 03, 011 (2018), arXiv:1801.04946 [hep-ph] .
  * Jimenez  _et al._ (2010) . Jimenez, T. Kitching, C. Pena-Garay, and L. Verde, “Can we measure the neutrino mass hierarchy in the sky?” JCAP 1005, 035 (2010), arXiv:1003.5918 [astro-ph.CO] .
  * Gerbino  _et al._ (2017) M. Gerbino, M. Lattanzi, O. Mena, and K. Freese, “A novel approach to quantifying the sensitivity of current and future cosmological datasets to the neutrino mass ordering through Bayesian hierarchical modeling,” Phys. Lett. B 775, 239–250 (2017), arXiv:1611.07847 [astro-ph.CO] .
  * Lattanzi and Gerbino (2018) M. Lattanzi and M. Gerbino, “Status of neutrino properties and future prospects - Cosmological and astrophysical constraints,” Front. in Phys. 5, 70 (2018), arXiv:1712.07109 [astro-ph.CO] .
  * De Salas  _et al._ (2018) P.F. De Salas, S. Gariazzo, O. Mena, C.A. Ternes, and M. Tórtola, “Neutrino Mass Ordering from Oscillations and Beyond: 2018 Status and Future Prospects,” Front. Astron. Space Sci. 5, 36 (2018), arXiv:1806.11051 [hep-ph] .
  * Mahony  _et al._ (2020) C. Mahony, B. Leistedt, H. V. Peiris, J. Braden, B. Joachimi, A. Korn, L. Cremonesi, and R. Nichol, “Target Neutrino Mass Precision for Determining the Neutrino Hierarchy,” Phys. Rev. D 101, 083513 (2020), arXiv:1907.04331 [astro-ph.CO] .
  * Archidiacono  _et al._ (2020) M. Archidiacono, S. Hannestad, and J. Lesgourgues, “What will it take to measure individual neutrino mass states using cosmology?” JCAP 09, 021 (2020), arXiv:2003.03354 [astro-ph.CO] .
  * Xu  _et al._ (2020) W. L. Xu, N. DePorzio, J. B. Muñoz, and C. Dvorkin, “Accurately Weighing Neutrinos with Cosmological Surveys,” (2020), arXiv:2006.09395 [astro-ph.CO] .
  * Stöcker  _et al._ (2020) P. Stöcker _et al._ (GAMBIT Cosmology Workgroup), “Strengthening the bound on the mass of the lightest neutrino with terrestrial and cosmological experiments,” (2020), arXiv:2009.03287 [astro-ph.CO] .
  * Hergt  _et al._ (2021) Lukas T. Hergt, Will J. Handley, Michael P. Hobson, and Anthony N. Lasenby, “Bayesian evidence for the tensor-to-scalar ratio $r$ and neutrino masses $m_{\nu}$: Effects of uniform vs logarithmic priors,” (2021), arXiv:2102.11511 [astro-ph.CO] .
  * Hassani  _et al._ (2016) F. Hassani, S. Baghram, and H. Firouzjahi, “Lensing as a Probe of Early Universe: from CMB to Galaxies,” JCAP 05, 044 (2016), arXiv:1511.05534 [astro-ph.CO] .
  * Manzotti (2018) A. Manzotti, “Future cosmic microwave background delensing with galaxy surveys,” Phys. Rev. D 97, 043527 (2018), arXiv:1710.11038 [astro-ph.CO] .
  * Sherwin  _et al._ (2011) B. D. Sherwin, J. Dunkley, S. Das, _et al._ , “Evidence for Dark Energy from the Cosmic Microwave Background Alone Using the Atacama Cosmology Telescope Lensing Measurements,” Phys. Rev. Lett. 107, 021302 (2011), arXiv:1105.0419 [astro-ph.CO] .
  * Percival and White (2009) W. J. Percival and M. White, “Testing cosmological structure formation using redshift-space distortions,” Mon. Not. Roy. Astron. Soc. 393, 297 (2009), arXiv:0808.0003 [astro-ph] .
  * Hernández (2017) O. F. Hernández, “Neutrino Masses, Scale-Dependent Growth, and Redshift-Space Distortions,” JCAP 06, 018 (2017), arXiv:1608.08298 [astro-ph.CO] .
  * Villaescusa-Navarro  _et al._ (2018) F. Villaescusa-Navarro, A. Banerjee, N. Dalal, E. Castorina, R. Scoccimarro, R. Angulo, and D. N. Spergel, “The imprint of neutrinos on clustering in redshift-space,” Astrophys. J. 861, 53 (2018), arXiv:1708.01154 [astro-ph.CO] .
  * Marulli  _et al._ (2011) F. Marulli, C. Carbone, M. Viel, L. Moscardini, and A. Cimatti, “Effects of massive neutrinos on the large-scale structure of the Universe,” Monthly Notices of the Royal Astronomical Society 418, 346–356 (2011), arXiv:1103.0278 [astro-ph.CO] .
  * Boyle and Komatsu (2018) A. Boyle and E. Komatsu, “Deconstructing the neutrino mass constraint from galaxy redshift surveys,” JCAP 03, 035 (2018), arXiv:1712.01857 [astro-ph.CO] .
  * Velten  _et al._ (2014) H. E. S. Velten, R. F. vom Marttens, and W. Zimdahl, “Aspects of the cosmological “coincidence problem”,” Eur. Phys. J. C 74, 3160 (2014), arXiv:1410.2509 [astro-ph.CO] .
  * Muir  _et al._ (2020) J. Muir _et al._ (DES), “DES Y1 results: Splitting growth and geometry to test $\Lambda$CDM,” (2020), arXiv:2010.05924 [astro-ph.CO] .
  * Kern  _et al._ (2017) Nicholas S. Kern, Adrian Liu, Aaron R. Parsons, Andrei Mesinger, and Bradley Greig, “Emulating simulations of cosmic dawn for 21 cm power spectrum constraints on cosmology, reionization, and x-ray heating,” The Astrophysical Journal 848, 23 (2017).
  * Lahav  _et al._ (2010) O. Lahav, A. Kiakotou, F. B. Abdalla, and C. Blake, “Forecasting neutrino masses from galaxy clustering in the Dark Energy Survey combined with the Planck measurements,” Monthly Notices of the Royal Astronomical Society 405, 168–176 (2010), arXiv:0910.4714 [astro-ph.CO] .
  * Brinckmann  _et al._ (2019) T. Brinckmann, D. C. Hooper, M. Archidiacono, J. Lesgourgues, and T. Sprenger, “The promising future of a robust cosmological neutrino mass measurement,” JCAP 01, 059 (2019), arXiv:1808.05955 [astro-ph.CO] .
  * Boyle  _et al._ (2020) A. Boyle, C. Uhlemann, O. Friedrich, A. Barthelemy, S. Codis, F. Bernardeau, C. Giocoli, and M. Baldi, “Nuw CDM cosmology from the weak lensing convergence PDF,” (2020), arXiv:2012.07771 [astro-ph.CO] .
  * Pearson and Zahn (2014) R. Pearson and O. Zahn, “Cosmology from cross correlation of CMB lensing and galaxy surveys,” Phys. Rev. D 89, 043516 (2014), arXiv:1311.0905 [astro-ph.CO] .
  * Giusarma  _et al._ (2018) E. Giusarma, S. Vagnozzi, S. Ho, S. Ferraro, K. Freese, R. Kamen-Rubio, and K.-B. Luk, “Scale-dependent galaxy bias, CMB lensing-galaxy cross-correlation, and neutrino masses,” Phys. Rev. D98, 123526 (2018), arXiv:1802.08694 [astro-ph.CO] .
  * Yu  _et al._ (2018) B. Yu, R. Z. Knight, B. D. Sherwin, S. Ferraro, L. Knox, and M. Schmittfull, “Towards Neutrino Mass from Cosmology without Optical Depth Information,” (2018), arXiv:1809.02120 [astro-ph.CO] .
  * Xu (2016) L. Xu, “Probing the neutrino mass through the cross correlation between the Rees-Sciama effect and weak lensing,” JCAP 08, 059 (2016), arXiv:1605.02403 [astro-ph.CO] .
  * Bayer  _et al._ (2021) Adrian E. Bayer, Francisco Villaescusa-Navarro, Elena Massara, Jia Liu, David N. Spergel, Licia Verde, Benjamin Wandelt, Matteo Viel, and Shirley Ho, “Detecting neutrino mass by combining matter clustering, halos, and voids,” (2021), arXiv:2102.05049 [astro-ph.CO] .
  * Saito  _et al._ (2008) Shun Saito, Masahiro Takada, and Atsushi Taruya, “Impact of massive neutrinos on the nonlinear matter power spectrum,” Phys. Rev. Lett. 100, 191301 (2008).
  * Brandbyge and Hannestad (2010) Jacob Brandbyge and Steen Hannestad, “Resolving cosmic neutrino structure: a hybrid neutrinoN-body scheme,”  2010, 021–021 (2010).
  * Lesgourgues  _et al._ (2009) Julien Lesgourgues, Sabino Matarrese, Massimo Pietroni, and Antonio Riotto, “Non-linear Power Spectrum including Massive Neutrinos: the Time-RG Flow Approach,” JCAP 06, 017 (2009), arXiv:0901.4550 [astro-ph.CO] .
  * Wong (2011) Y. Y. Y. Wong, “Neutrino mass in cosmology: status and prospects,” Ann. Rev. Nucl. Part. Sci. 61, 69–98 (2011), arXiv:1111.1436 [astro-ph.CO] .
  * Bird  _et al._ (2012) Simeon Bird, Matteo Viel, and Martin G. Haehnelt, “Massive neutrinos and the non-linear matter power spectrum,” MNRAS 420, 2551–2561 (2012), arXiv:1109.4416 [astro-ph.CO] .
  * Adamek  _et al._ (2017) J. Adamek, R. Durrer, and M. Kunz, “Relativistic N-body simulations with massive neutrinos,” JCAP 11, 004 (2017), arXiv:1707.06938 [astro-ph.CO] .
  * Villaescusa-Navarro  _et al._ (2020) F. Villaescusa-Navarro _et al._ , “The Quijote simulations,” Astrophys. J. Suppl. 250, 2 (2020), arXiv:1909.05273 [astro-ph.CO] .
  * Mead  _et al._ (2020) A. Mead, S. Brieden, T. Tröster, and C. Heymans, “HMcode-2020: Improved modelling of non-linear cosmological power spectra with baryonic feedback,” (2020), arXiv:2009.01858 [astro-ph.CO] .
  * Nascimento and Loverde (2021) Caio Bastos de Senna Nascimento and Marilena Loverde, “Neutrinos in N-body simulations,” (2021), arXiv:2102.05690 [astro-ph.CO] .
  * Bose  _et al._ (2021) Benjamin Bose, Bill S. Wright, Matteo Cataneo, Alkistis Pourtsidou, Carlo Giocoli, Lucas Lombriser, Ian G. McCarthy, Marco Baldi, Simon Pfeifer, and Qianli Xia, “On the road to percent accuracy V: the non-linear power spectrum beyond $\Lambda$CDM with massive neutrinos and baryonic feedback,” (2021), arXiv:2105.12114 [astro-ph.CO] .
  * Castorina  _et al._ (2015) E. Castorina, C. Carbone, J. Bel, E. Sefusatti, and K. Dolag, “DEMNUni: The clustering of large-scale structures in the presence of massive neutrinos,” JCAP 07, 043 (2015), arXiv:1505.07148 [astro-ph.CO] .
  * Raccanelli  _et al._ (2019) A. Raccanelli, L. Verde, and F. Villaescusa-Navarro, “Biases from neutrino bias: to worry or not to worry?” Mon. Not. Roy. Astron. Soc. 483, 734–743 (2019), arXiv:1704.07837 [astro-ph.CO] .
  * Vagnozzi  _et al._ (2018b) S. Vagnozzi, T. Brinckmann, M. Archidiacono, K. Freese, M. Gerbino, J. Lesgourgues, and T. Sprenger, “Bias due to neutrinos must not uncorrect’d go,” JCAP 09, 001 (2018b), arXiv:1807.04672 [astro-ph.CO] .
  * Chiang  _et al._ (2019) Chi-Ting Chiang, Marilena LoVerde, and Francisco Villaescusa-Navarro, “First detection of scale-dependent linear halo bias in $N$-body simulations with massive neutrinos,” Phys. Rev. Lett. 122, 041302 (2019), arXiv:1811.12412 [astro-ph.CO] .