# Strengthening the bound on the mass of the lightest neutrino
with terrestrial and cosmological experiments

The GAMBIT Cosmology Workgroup: Patrick Stöcker<EMAIL_ADDRESS>aachen.de Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH
Aachen University, D-52056 Aachen, Germany  Csaba Balázs  School of Physics
and Astronomy, Monash University, Melbourne, VIC 3800, Australia  Sanjay Bloor
School of Mathematics and Physics, The University of Queensland, St. Lucia,
Brisbane, QLD 4072, Australia  Department of Physics, Imperial College London,
Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK  Torsten
Bringmann  Department of Physics, University of Oslo, Box 1048, Blindern,
N-0316 Oslo, Norway  Tomás E. Gonzalo  School of Physics and Astronomy, Monash
University, Melbourne, VIC 3800, Australia  Will Handley  Cavendish
Laboratory, University of Cambridge, JJ Thomson Avenue, Cambridge, CB3 0HE, UK
Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge, CB3 0HA, UK  Gonville & Caius College, Trinity Street, Cambridge,
CB2 1TA, UK  Selim Hotinli  Department of Physics, Imperial College London,
Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK  Cullan Howlett
<EMAIL_ADDRESS>School of Mathematics and Physics, The University of
Queensland, St. Lucia, Brisbane, QLD 4072, Australia  Felix Kahlhoefer
Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen
University, D-52056 Aachen, Germany  Janina J. Renk<EMAIL_ADDRESS>School of Mathematics and Physics, The University of Queensland, St. Lucia,
Brisbane, QLD 4072, Australia  Department of Physics, Imperial College London,
Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK  Oskar Klein
Centre for Cosmoparticle Physics, AlbaNova University Centre, SE-10691
Stockholm, Sweden  Pat Scott<EMAIL_ADDRESS>School of Mathematics and
Physics, The University of Queensland, St. Lucia, Brisbane, QLD 4072,
Australia  Department of Physics, Imperial College London, Blackett
Laboratory, Prince Consort Road, London SW7 2AZ, UK  Aaron C. Vincent
Department of Physics, Engineering Physics and Astronomy, Queen’s University,
Kingston ON K7L 3N6, Canada  Arthur B. McDonald Canadian Astroparticle Physics
Research Institute, Kingston ON K7L 3N6, Canada  Perimeter Institute for
Theoretical Physics, Waterloo ON N2L 2Y5, Canada  Martin White  ARC Centre for
Dark Matter Particle Physics, Department of Physics, University of Adelaide,
Adelaide, SA 5005, Australia

###### Abstract

We determine the upper limit on the mass of the lightest neutrino from the
most robust recent cosmological and terrestrial data. Marginalizing over
possible effective relativistic degrees of freedom at early times
($N_{\mathrm{eff}}$) and assuming normal mass ordering, the mass of the
lightest neutrino is less than 0.037 eV at 95% confidence; with inverted
ordering, the bound is 0.042 eV. These results improve upon the strength and
robustness of other recent limits and constrain the mass of the lightest
neutrino to be barely larger than the largest mass splitting. We show the
impacts of realistic mass models, and different sources of $N_{\mathrm{eff}}$.

††preprint: TTK-20-28, gambit-physics-2020

##  I Introduction

Neutrino masses are arguably the most concrete evidence to date of physics
beyond the Standard Model (SM). Measurements of their flavor oscillations at
reactor Gando  _et al._ (2013); An  _et al._ (2017); Adey  _et al._ (2018);
Serra ; Bak  _et al._ (2018), accelerator Adamson  _et al._ (2013a, b);
Izmaylov ; Friend ; Sanchez (2018); Acero  _et al._ (2019), solar Cleveland
_et al._ (1998); Kaether  _et al._ (2010); Abdurashitov  _et al._ (2009);
Aharmim  _et al._ (2013); Hosaka  _et al._ (2006); Cravens  _et al._ (2008);
Abe  _et al._ (2011); Ikeda (2018); Bellini  _et al._ (2011, 2010, 2014) and
atmospheric Aartsen  _et al._ (2015); Abe  _et al._ (2018) experiments show
that at least two of the three SM neutrinos must be massive. While oscillation
experiments probe mass-squared differences between eigenstates, the expansion
history of the Universe, growth of cosmic structure and the cosmic microwave
background (CMB) are sensitive to absolute masses, which determine when a
neutrino becomes nonrelativistic; for most cosmological applications, this is
expressed in terms of the sum of masses $\sum m_{\nu}$. Robust and precise
inference on the mass of the lightest state can therefore only be obtained by
combining the latest results of all these probes self-consistently, including
associated uncertainties from each, as well as constraints on other relevant
parameters from e.g. big bang nucleosynthesis (BBN) and late-time cosmological
observables Wong (2011); Lesgourgues and Pastor (2014); Vagnozzi  _et al._
(2017); Aghanim  _et al._ (2020a); Loureiro  _et al._ (2019); Ivanov  _et al._
(2020); Archidiacono  _et al._ (2020).

As no probe has yet directly measured the mass of a single neutrino, the most
convenient three-flavor parametrisation is in terms of the mass $m_{\nu_{0}}$
of the lightest neutrino and two squared mass splittings, $\Delta
m_{21}^{2}\equiv m_{2}^{2}-m_{1}^{2}$ and $\Delta m_{3l}^{2}\equiv
m_{3}^{2}-m_{l}^{2}$. Here 1, 2 and 3 label the mass eigenstates with the
largest component of $\nu_{e}$, $\nu_{\mu}$ and $\nu_{\tau}$, respectively.
Two mass orderings are presently permitted by the data: the normal
($m_{1}<m_{2}\ll m_{3}$; NH) and inverted hierarchies ($m_{3}\ll m_{1}<m_{2}$;
IH). $\Delta m_{3l}^{2}$ refers to the splitting between the lightest and
heaviest states, i.e. $l=1$ for the NH and $l=2$ for the IH. In terms of the
splitting parameters, the physical masses are

$\displaystyle\mathrm{NH:}\ $ $\displaystyle(m_{1}^{2},m_{2}^{2},m_{3}^{2})=$
$\displaystyle m_{\nu_{0}}^{2}+(0,\Delta m_{21}^{2},\Delta m_{3l}^{2})$ (1)
$\displaystyle\mathrm{IH:}\ $ $\displaystyle(m_{3}^{2},m_{1}^{2},m_{2}^{2})=$
$\displaystyle m_{\nu_{0}}^{2}+(0,|\Delta m_{3l}^{2}|-\Delta
m_{21}^{2},|\Delta m_{3l}^{2}|).$

In this article, we make use of the new cosmology module CosmoBit Renk  _et
al._ (2021) within the beyond-the-SM global fitting package GAMBIT Athron  _et
al._ (2017) in order to perform the most precise and robust combination to
date of cosmological and experimental constraints on the mass of the lightest
neutrino. We include the most recent CMB likelihoods from _Planck_ Aghanim
_et al._ (2020b), recent three-flavor neutrino global fit results from NuFit
Esteban  _et al._ (2019), and correlated measurements of the baryon acoustic
oscillation (BAO) scale by 6dF Beutler  _et al._ (2011), SDSS-MGS Ross  _et
al._ (2015), BOSS DR12 Alam  _et al._ (2017), eBOSS DR14 Ata  _et al._ (2018);
Bautista  _et al._ (2018); Blomqvist  _et al._ (2019) and DES Abbott  _et al._
(2019). We also compute and propagate the primordial helium abundance and
number of effective relativistic degrees of freedom $N_{\mathrm{eff}}$ through
all our calculations and likelihoods self-consistently and account for the
uncertainty on the lifetime of the neutron. When computing bounds on neutrino
masses, we illustrate the impact of different physical assumptions about
$N_{\mathrm{eff}}$ and show the impact of the choice of neutrino mass model on
the derived value of the Hubble parameter $H_{0}$, of particular interest
given the present tension between expansion measurements at early and late
times Riess  _et al._ (2019); Wong  _et al._ (2020); Freedman  _et al._
(2019).

##  II Methodology

Our likelihoods are based on the latest and most constraining data implemented
in CosmoBit Renk  _et al._ (2021), NeutrinoBit Chrzaszcz  _et al._ (2020) and
PrecisionBit Athron  _et al._ (2018a):

  * (i)

Neutrino oscillations: Two-dimensional NH and IH $\Delta\chi^{2}$ tables for
$\Delta m_{3l}^{2}$ and $\Delta m_{21}^{2}$ from NuFit 4.1 Esteban  _et al._
(2019). These come from fits to the data from solar (Homestake chlorine
Cleveland  _et al._ (1998), Gallex/GNO Kaether  _et al._ (2010), SAGE
Abdurashitov  _et al._ (2009), SNO Aharmim  _et al._ (2013), four phases of
Super-Kamiokande Hosaka  _et al._ (2006); Cravens  _et al._ (2008); Abe  _et
al._ (2011); Ikeda (2018), two phases of Borexino Bellini  _et al._ (2011,
2010, 2014)), atmospheric (IceCube/DeepCore Aartsen  _et al._ (2015), Super-
Kamiokande Abe  _et al._ (2018)), reactor (KamLAND Gando  _et al._ (2013),
Double Chooz Serra , Daya Bay An  _et al._ (2017); Adey  _et al._ (2018), Reno
Bak  _et al._ (2018)), and accelerator experiments (MINOS Adamson  _et al._
(2013a, b), T2K Izmaylov ; Friend , NO$\nu$A Sanchez (2018); Acero  _et al._
(2019)). The other oscillation parameters (mixing angles $\theta_{ij}$ and CP
phase $\delta_{CP}$) have no bearing on neutrino mass studies and do not enter
our analysis. Updates contained in NuFit 5.0 Esteban  _et al._ (2020) released
when this article was in the final stages of preparation include only small
improvements to the likelihoods for $\Delta m_{3l}^{2}$ and $\Delta
m_{21}^{2}$, so they have minimal impact on the results we show here.

  * (ii)

BBN: Primordial abundances of 4He, $Y_{\mathrm{p}}=0.245\pm 0.003$ Tanabashi
_et al._ (date) and deuterium, ${\rm D/H}=(2.527\pm 0.030)\times 10^{-5}$
Cooke  _et al._ (2018).

  * (iii)

CMB: _Planck_ 2018 baseline likelihoods consisting of high-$\ell$ and
low-$\ell$ temperature and polarization data, plus CMB lensing Aghanim  _et
al._ (2020b).

  * (iv)

Supernovae Type Ia (SN Ia): $1048$ SN Ia included in the Pantheon compilation
Scolnic  _et al._ (2018).

  * (v)

BAO scale: Measurements of the transverse comoving distance $D_{M}$ and the
Hubble parameter $H(z)$ from the BOSS DR12 anisotropic consensus Alam  _et
al._ (2017), $D_{M}$ from DES Y1 Abbott  _et al._ (2019), and the volume-
averaged distance $D_{V}$ from the combined 6dF and MGS galaxy surveys Beutler
_et al._ (2011); Ross  _et al._ (2015); Carter  _et al._ (2018) and the eBOSS
DR14 lumious red galaxy (LRG) and quasi-stellar objects (QSO) samples Ata  _et
al._ (2018); Bautista  _et al._ (2018). All measurements are relative to
$r_{s}$, the radius of the sound horizon at the baryon-drag epoch.

We carefully consider correlations between overlapping samples in the BAO
scale measurements. The 6dF+MGS result can be considered independent of the
others, as these samples do not overlap in redshift with the others.
Similarly, we treat the DES results as independent of all others, as less than
$10\%$ of the DES footprint overlaps with BOSS DR12 or eBOSS DR14. The DES
sample also consists of very different targets to BOSS and eBOSS, and uses a
different methodology (photometric rather than spectroscopic redshifts).
However, the BOSS DR12 and eBOSS results are correlated: The eBOSS LRG sample
actually contains some of the same galaxies as BOSS DR12, while the eBOSS QSO
sample overlaps substantially with the LRGs both on the sky and in redshift.
Overall, there are nonzero correlations that should be accounted for between
the measurements $(D_{M}/r_{s})^{\mathrm{BOSS}}$, $(Hr_{s})^{\mathrm{BOSS}}$,
$(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$, and
$(D_{V}/r_{s})^{\mathrm{eBOSS,QSO}}$.

To do this in a way that accounts for variation with cosmological parameters,
we implement a novel method to compute the cross-correlation coefficients
using Fisher matrices, following BAO forecasting techniques Seo and Eisenstein
(2007); Font-Ribera  _et al._ (2014a). We sum the Fisher information that each
sample contributes to the four overlapping measurements listed above,
accounting for redshift and sky overlap. Inverting the full Fisher matrix then
gives the correlation coefficients. We do this separately for every
combination of cosmological parameters in the fit, using the number density of
objects, matter power spectrum and growth rate of structure to model the
BOSS/eBOSS galaxy power spectra and their covariance matrices as a function of
redshift. We split the models into smooth and oscillatory components to obtain
the derivatives of the BAO feature in the power spectrum with respect to the
distance measurements, and fix the galaxy bias and nonlinear damping of the
power spectra to their best-fit values as reported in the original works. The
Fisher matrix calculation then uses these all as inputs, integrating over
angles and scales in the clustering measurements consistent with the range
used in the original measurements.

We find a value for the cross-correlation between BOSS DR12 and the eBOSS LRGs
comparable to the one reported in Ref. Bautista  _et al._ (2018). The benefits
of our technique are that it includes information from all scales included in
the measurement and does so self-consistently for each combination of
cosmological parameters. More details on the calculation and the computed
correlation coefficients can be found in Appendix A.

We compute the BAO scale and SN Ia likelihoods via an interface to MontePython
3.3.0 Audren  _et al._ (2013); Brinckmann and Lesgourgues (2019); our novel
BAO scale correlation treatment will appear in a future release. For computing
observable predictions, we use routines in CosmoBit and associated interfaces
to AlterBBN 2.2 (Arbey (2012); Arbey  _et al._ (2020); for BBN yields) and
CLASS 2.9.3 (Blas  _et al._ (2011); for solving the background cosmology and
Boltzmann equations).

Sector | Parameter | Range
---|---|---
$\nu$ masses | $m_{\nu_{0}}$ | [0, 1.1] eV
| $\Delta m^{2}_{21}$ |  $[6,\ 9]\times 10^{-5}$ eV2
(NH) | $\Delta m^{2}_{3l}$ |  $[2.2,\ 2.8]\times 10^{-3}$ eV2
(IH) | $\Delta m^{2}_{3l}$ |  $[-2.8,\ -2.2]\times 10^{-3}$ eV2
$\Lambda$CDM | $H_{0}$ | [50, 80] km s-1 Mpc-1
| $\Omega_{\mathrm{b}}h^{2}$ | [0.020, 0.024]
| $\Omega_{\mathrm{cdm}}h^{2}$ | [0.10, 0.15]
| $\tau_{\mathrm{reionization}}$ | [0.004, 0.20]
| $\mathrm{ln}\left(10^{10}\,A_{s}\right)$ | [2.5, 3.5]
| $n_{s}$ | [0.90, 1.10]
$N_{\mathrm{eff}}$ | $r_{\nu}$ | [0.75, 1.15]
Nuisance | SN Ia abs. magnitude $M$ | [$-$20, $-$18]
| Neutron lifetime $\tau_{\mathrm{n}}$ | [876, 883] s
|  Planck likelihood | 21 parameters varied
Table 1: Parameters and ranges varied in the main analysis of this article.
All parameters are sampled with linear priors. For the nuisance parameters
associated with the Planck likelihood, we adopt the same prior ranges as done
for the Planck baseline analysis Aghanim  _et al._ (2020b) and apply the
recommended Gaussian priors as additional likelihood contributions.

We perform separate fits of the NH and IH, varying $m_{\nu_{0}}$, $\Delta
m_{21}^{2}$ and $\Delta m_{3l}^{2}$, the six free parameters of the standard
Lambda cold dark matter cosmology ($\Lambda$CDM; see Renk  _et al._ (2021) for
detailed definitions), and $\Delta N_{\mathrm{eff}}\equiv
N_{\mathrm{eff}}-N_{\mathrm{SM}}$ with $N_{\mathrm{SM}}=3.045$ de Salas and
Pastor (2016); Akita and Yamaguchi (2020); Froustey  _et al._ (2020) (Table
1). For our main analysis, we conservatively adopt a linear prior on
$m_{\nu_{0}}$ between 0 and 1.1 eV. In Appendix B, we show how the limit on
$m_{\nu_{0}}$ strengthens if we change to a logarithmic prior above
$m_{\nu_{0}}=0.0003$ eV. We adopt linear priors on the $\Lambda$CDM
parameters, as these are sufficiently well constrained by data that their
priors are inconsequential.

We assume that $N_{\mathrm{eff}}$ has the same value during BBN and
recombination, varying $\Delta N_{\mathrm{eff}}$ by scanning linearly over the
effective neutrino temperature ratio $r_{\nu,\mathrm{eff}}\equiv
T_{\nu}/T_{\nu,\text{SM}}=(\Delta
N_{\mathrm{eff}}/N_{\mathrm{SM}}+1)^{\frac{1}{4}}$. This approach enables us
to explore the full range of the number of effective relativistic degrees of
freedom in the early Universe, corresponding to both positive and negative
values of $\Delta N_{\mathrm{eff}}$. We later explore the impact of
restricting our analysis to $\Delta N_{\mathrm{eff}}\geq 0$, in which case,
$\Delta N_{\mathrm{eff}}$ can be interpreted as the contribution of additional
ultrarelativistic (radiation) species and to the pure SM case ($\Delta
N_{\mathrm{eff}}=0$).

Finally, we vary a total of 23 nuisance parameters representing uncertainties
in the SN Ia absolute magnitude, the Planck analysis, and the neutron
lifetime. These nuisance parameters are constrained, respectively, with the
MontePython likelihood for the magnitude of Pantheon supernovae, a likelihood
implementation of the Planck nuisance priors Aghanim  _et al._ (2020b), and
the combination of all “bottle” measurements of the lifetime of the neutron
$\tau_{n\,,\text{bottle}}=879.4\pm 0.6$ s Zyla  _et al._ (2020). We employ the
nested sampler PolyChord 1.17.1 Handley  _et al._ (2015) in fast-slow mode via
ScannerBit Martinez  _et al._ (2017) in order to oversample Planck and SN Ia
nuisances. Our fits use 500 live points, 5000 initial samples from the prior,
a stopping tolerance of 0.01, $n_{\mathrm{repeats}}=2n_{\mathrm{slow}}=22$, a
1:3 fast-slow timing split (leading to approximately 340 likelihood
evaluations with different nuisance parameters per combination of the
remaining parameters), and default values for all other settings.

##  III Results

Assuming $\Lambda$CDM cosmology plus a free $N_{\mathrm{eff}}$ and normal mass
ordering, we find a global 95% confidence upper bound on the lightest neutrino
mass of $m_{\nu_{0}}<0.037$ eV; for inverted ordering, this increases slightly
to $m_{\nu_{0}}<0.042$ eV. In terms of the sum of neutrino masses, this
corresponds to $0.058<\sum m_{\nu}/\mathrm{eV}<0.139$ (NH) and $0.098<\sum
m_{\nu}/\mathrm{eV}<0.174$ (IH).

Figure 1: 1D and joint 2D posteriors (bottom left) and profile likelihoods
(upper right) on the lightest neutrino mass, the sum of neutrino masses, the
number of effective neutrino species at the time of CMB formation, and the
Hubble parameter, based on the most robust and complete combination to date of
CMB, BAO scale, SN Ia, BBN and neutrino oscillation data. Posteriors are shown
for normal and inverted neutrino hierarchies and are compared with the often-
seen (but unphysical) scenario of three degenerate-mass neutrinos. Profile
likelihoods are for the normal hierarchy only and assume the best-fit values
from the corresponding posterior scan for the 21 Planck nuisance parameters.
Shading indicates 68% and 95% credible/confidence regions.

The lower-left triangle of Figure 1 compares the one- and two-dimensional
marginalized posterior distributions for the normal and inverted hierarchies.
Here we show both the lightest neutrino mass and the sum of neutrino masses,
as well as their correlations with $N_{\mathrm{eff}}$ and $H_{0}$. We also
show the result for the canonical scenario considered in most previous
analyses, such as those by both Planck Aghanim  _et al._ (2020a) and eBOSS
Alam  _et al._ (2021), where a single parameter specifies a common degenerate
mass for all three neutrinos. In all cases, the maximum posterior probability
density is achieved for a massless lightest neutrino, reflecting the fact that
there remains no positive cosmological hint for neutrino mass to date. The
degenerate-mass assumption would lead one to erroneously infer that $\sum
m_{\nu}/\mathrm{eV}<0.115$ at 95% confidence. This result is plainly biased
toward lower values due to the fact that the majority of the probability
distribution lies within the unphysical region excluded by neutrino
oscillation experiments. Using a physically realistic mass model shifts the
95% interval for $H_{0}$ from $67.7\pm 1.7$ km s-1 Mpc-1 to $67.5\pm 1.8$ km
s-1 Mpc-1 (NH) or $67.4\pm 1.7$ km s-1 Mpc-1 (IH), and the 95% interval for
$N_{\mathrm{eff}}$ from $3.04\pm 0.24$ to $3.06\pm 0.24$ (NH) or $3.08\pm
0.24$ (IH).

For comparison, in the upper-right triangle of Figure 1, we also show prior-
independent profile likelihoods for the NH obtained from the differential
evolution sampler Diver Martinez  _et al._ (2017) with a population of
$10^{4}$, a convergence threshold of $10^{-4}$, and the Planck nuisance
parameters fixed to their best-fit values from the NH PolyChord fit (as Diver
has no fast-slow feature). The results match the posteriors reasonably
closely, but give slightly stronger implied constraints at 95% confidence:
$m_{\nu_{0}}<0.033$ eV and $0.058<\sum m_{\nu}/\mathrm{eV}<0.127$.111The
slightly higher profile likelihood than posterior for much of the allowed
range of $\sum m_{\nu}$ should be understood in the context of frequentist
confidence levels deriving from isolikelihood contours, rather than integrated
posterior probabilities as in the case of Bayesian credible regions. These
findings confirm the robustness of our main (Bayesian) results. We have also
calculated Bayesian evidences and find that the IH is disfavored relative to
the NH with Bayes factors between $\log B=5.6$ and $\log B=7.0$ (depending on
the treatment of $N_{\text{eff}}$), a result driven mostly by the neutrino
oscillation likelihoods de Salas  _et al._ (2021).

Figure 2: Comparison of posterior probabilities for the mass of the lightest
neutrino under three different assumptions: $\Delta N_{\mathrm{eff}}=0$,
$\Delta N_{\mathrm{eff}}\geq 0$ corresponding to dark radiation, or $\Delta
N_{\mathrm{eff}}$ free to take on positive or negative values, corresponding
to a modified effective neutrino temperature.

In Figure 2, we examine the impacts of different physical sources of
$N_{\mathrm{eff}}$: (a) changes in the neutrino temperature, where $\Delta
N_{\mathrm{eff}}$ is allowed to be positive or negative, as in our benchmark
analyses, (b) dark radiation, where $\Delta N_{\mathrm{eff}}\geq 0$, and (c)
the pure SM case, where $\Delta N_{\mathrm{eff}}=0$. The resulting posteriors
only change very slightly, corresponding to shifts of the order of 0.002–0.003
eV in the 95% limit on $m_{\nu_{0}}$. Our final results can therefore be
considered rather robust to assumptions about $\Delta N_{\mathrm{eff}}$.
Figure 2 suggests that at the 68% confidence level, allowing positive $\Delta
N_{\mathrm{eff}}$ may weaken the limits slightly compared to $\Delta
N_{\mathrm{eff}}=0$, with the effect offset to some extent by also allowing
$\Delta N_{\mathrm{eff}}<0$. Any such effect is however small enough that it
is difficult to distinguish from sampling noise. Notably, bounds on
$m_{\nu_{0}}$ and $\sum m_{\nu}$ can be weakened substantially in cosmologies
featuring both dark radiation and a modified neutrino temperature Renk  _et
al._ (2021), neutrino self-interactions, or exotic dark energy, even to the
level where a direct measurement of the neutrino mass may be within reach of
the KATRIN experiment Aker  _et al._ (2019).

##  IV Discussion

Our analysis provides a more precise and robust limit on the sum of neutrino
masses than either those of Planck Aghanim  _et al._ (2020a) or eBOSS Alam
_et al._ (2021), mainly due to our use of physical neutrino mass models rather
than the assumption of degenerate masses. Comparing results for the unphysical
degenerate-mass model however provides an indication of the constraining power
of the cosmological data used in each case. Our limit ($\sum
m_{\nu}/\mathrm{eV}<0.115$) is very close to the most similar combination in
Ref. Aghanim  _et al._ (2020a) ($\sum m_{\nu}/\mathrm{eV}<0.12$, Planck
TT,TE,EE+lowE+lensing+BAO, with varying $N_{\rm{eff}}$), indicating that
within $\Lambda$CDM+$N_{\text{eff}}$ the eBOSS DR14 and DES BAO and Pantheon
SN Ia measurements only add limited additional information.

The most similar eBOSS limit ($\sum m_{\nu}/\mathrm{eV}<0.099$) is stronger.
Reference Alam  _et al._ (2021) includes slightly more up-to-date eBOSS BAO
scale measurements compared to our analysis, but the improved sensitivity is
mainly driven by their inclusion of redshift-space distortions (RSD). We do
not include RSD measurements, as to date they have been based on templates for
the matter power spectrum that assume a particular neutrino mass, and the fits
then neglect the scale dependence in the growth rate of structure. As such we
believe they cannot robustly be used to constrain neutrino masses. Reference
Alam  _et al._ (2021) also includes new (DR16) Ly$\alpha$ constraints on the
BAO scale. We checked that including DR14 Ly$\alpha$ measurements has no
impact on our limit. Given this result, the fact that the redshifts probed by
Ly$\alpha$ data are intermediate between those of the CMB and other BAO
measurements, that Ly$\alpha$ BAO are slightly discrepant with other BAO, and
that the Ly$\alpha$ BAO results require more precise control over
observational and astrophysical systematics than galaxy BAO Blomqvist  _et
al._ (2019), we argue that excluding the DR14 Ly$\alpha$ result gives a more
robust limit on neutrino masses at no decrease in the statistical constraining
power.

The limits that we present here on the mass of the lightest neutrino are
almost 60% stronger than those derived from a combined fit to both mass
hierarchies in a recent similar analysis ($m_{\nu_{0}}<0.086$ eV Loureiro  _et
al._ (2019)), and 14% stronger than those appearing in a similar
contemporaneous analysis de Salas  _et al._ (2021). At the limiting value of
$m_{\nu_{0}}$, this brings the absolute scale of neutrino masses down to a
level comparable to the larger of the two mass splittings ($m_{3}-m_{1}=0.025$
eV in the NH, $m_{2}-m_{3}=0.023$ eV in the IH). We use the same BBN and SN Ia
data as Ref. Loureiro  _et al._ (2019) but improved neutrino and CMB data:
results from NuFit 4.1 rather than 2.1 for neutrino experiments (leading to
roughly 20% stronger bounds), and a CMB likelihood based on 2018 rather than
2015 Planck data (leading to roughly 30% stronger bounds). Compared to de
Salas  _et al._ (2021), we add BBN constraints, DES, and eBOSS BAO scale data
and marginalize over $N_{\mathrm{eff}}$. We also propagate the primordial
helium abundance fully and incorporate the uncertainty on the lifetime of the
neutron (see also discussion in Ref. Renk  _et al._ (2021)). From galaxy
surveys, we and Ref. de Salas  _et al._ (2021) rely exclusively on scale
measurements, as a correct statistical combination of scale data can provide a
strong limit on neutrino masses that is also very robust. Reference Loureiro
_et al._ (2019) instead used an angular clustering reanalysis of the BOSS DR12
data. Their approach used the full shape of the clustering more self-
consistently than template-based RSD measurements but still required several
uncertainties such as the galaxy bias, redshift error dispersions, and
spectroscopic redshift errors to be modeled, and data cuts to be made, adding
28 more nuisance parameters. In principle, nonlinear scales should provide the
most constraining power on neutrino mass parameters but difficulty in the
modeling limits the accessibility of this information Font-Ribera  _et al._
(2014b). With current analysis techniques, BAO scale measurements still add
more constraining power for neutrino masses than data encoding the full shape
of the galaxy power spectrum—although this is not expected to remain true for
much longer Ivanov  _et al._ (2020).

##  V Summary

We have presented a comprehensive combined analysis of recent neutrino
oscillation, CMB, SN Ia, BBN, and BAO scale data deriving the most accurate
and precise limits to date on the mass of the lightest neutrino and the sum of
neutrino masses. Assuming normal mass ordering and standard cosmology plus
$\Delta N_{\mathrm{eff}}\neq 0$, we have found $m_{\nu_{0}}<0.037$ eV and
$0.058<\sum m_{\nu}/\mathrm{eV}<0.139$. With inverted ordering,
$m_{\nu_{0}}<0.042$ eV and $0.098<\sum m_{\nu}/\mathrm{eV}<0.174$. These
results should serve as a benchmark in the coming years, as neutrino cosmology
continues its inexorable progress toward a measurement of the absolute
neutrino mass scale.

All input files and parameter samples produced for this article can be found
on Zenodo Cos .

## Acknowledgements

We thank the GAMBIT Community, Alex Arbey, Thejs Brinckmann, Tamara Davis,
Martina Gerbino, Deanna Hooper, Julien Lesgourgues, Vivian Poulin, Nils
Schöneberg, Jesús Torrado and Sunny Vagnozzi for helpful discussions, PRACE
for access to Joliot-Curie at CEA, RWTH Aachen University for access to JARA
under Project No. jara0184 and the University of Cambridge for access to CSD3
resources. P. St. and F. K. acknowledge funding from Deutsche
Forschungsgemeinschaft Grant No. KA 4662/1-1, C. B., T. E. G., M. J. W., P.
Sc. and C. H. from Australian Research Cuncil Grants No. DP180102209, No.
FL180100168, No. CE200100008, and No. FT190100814, J. J. R. from Swedish
Research Council Contract No. 638-2013-8993, W. H. from a Gonville & Caius
Research Fellowship and the George Southgate visiting fellowship, and A. C. V.
from the Arthur B. McDonald Canadian Astroparticle Physics Research Institute,
Canada Foundation for Innovation and Ontario Ministry of Economic Development,
Job Creation and Trade (MEDJCT). Research at Perimeter Institute is supported
by the Government of Canada through the Department of Innovation, Science, and
Economic Development, and by the Province of Ontario through MEDJCT. This
article made use of matplotlib Hunter (2007), GetDist Lewis (2019), and pippi
Scott (2012).

##  Appendix A Details of BAO scale correlation coefficients

Here we provide more detail on our novel method to account for the
correlations between overlapping BAO experiments in a cosmology-dependent way.
We acknowledge that the use of cosmology-dependent Gaussian covariance
matrices can lead to an overestimation of the Fisher information (Carron,
2013); however, in our method we are only computing the cross-correlation
coefficients, which are then scaled by the original measurement errors. We do
not modify these original errors and so are not underestimating the
uncertainty. More so, we argue that including cosmology dependence in the
correlation coefficients is the correct thing to do, as the overlap between
BAO measurements is mainly related to the (effective) cosmological volume
shared by the experiments, which is clearly a function of the cosmological
model.

We begin by expressing the joint covariance matrix of BAO measurements from
the overlapping BOSS DR12 and eBOSS results $\bm{\mathsf{C}}$ in terms of the
inverse of the Fisher matrix $\bm{\mathsf{F}}$,

$\bm{\mathsf{F}}^{-1}=\bm{\mathsf{C}}=\bm{\mathsf{E}}\begin{pmatrix}1&\rho^{\mathrm{BOSS}}_{D_{M}/r_{s},Hr_{s}}&c_{0}&0\\\
\rho^{\mathrm{BOSS}}_{D_{M}/r_{s},Hr_{s}}&1&c_{1}&0\\\ c_{0}&c_{1}&1&c_{2}\\\
0&0&c_{2}&1\\\ \end{pmatrix}\bm{\mathsf{E}}.$ (2)

Here, $c_{0}$ is the correlation coefficient between measurements
$(D_{M}/r_{s})^{\mathrm{BOSS}}$ and $(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$,
$c_{1}$ is the correlation between $(Hr_{s})^{\mathrm{BOSS}}$ and
$(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$, and $c_{2}$ is the correlation between
$(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$ and $(D_{V}/r_{s})^{\mathrm{eBOSS,QSO}}$.
The cross-correlation $\rho^{BOSS}_{D_{M}/r_{s},Hr_{s}}$ between the BOSS DR12
measurements of $D_{M}/r_{s}$ and $Hr_{s}$ is provided as part of the BOSS
DR12 results in MontePython and is not replaced in this analysis, or made
dependent upon cosmological parameters. The matrix
$\bm{\mathsf{E}}=\mathrm{diag}(\sigma^{\mathrm{BOSS}}_{D_{M}/r_{s}},\sigma^{\mathrm{BOSS}}_{Hr_{s}},\sigma^{\mathrm{eBOSS,LRG}}_{D_{V}/r_{s}},\sigma^{\mathrm{eBOSS,QSO}}_{D_{V}/r_{s}})$
contains the experimental uncertainties from the BAO measurements, which are
also included in MontePython. Hence the only unknowns are $c_{i}$, which we
calculate for each set of cosmological parameters by inverting the full Fisher
matrix. In practice, we compute the Fisher matrix for the BAO scale parameters

$\alpha=\frac{D_{V}r^{\mathrm{fid}}_{s}}{D^{\mathrm{fid}}_{V}r_{s}},\quad\alpha_{\perp}=\frac{D_{M}r^{\mathrm{fid}}_{s}}{D^{\mathrm{fid}}_{M}r_{s}},\quad\alpha_{||}=\frac{H^{\mathrm{fid}}r^{\mathrm{fid}}_{s}}{Hr_{s}},$
(3)

where “fid” corresponds to the fiducial cosmology used to make the original
clustering measurements. It is trivial to convert the correlation coefficients
to those for the distance scales that we actually fit in MontePython.

We do this based on Howlett  _et al._ (2017), writing the Fisher matrix
elements for each parameter of interest as a sum over the information from the
different redshift bins and overlapping/nonoverlapping sky areas for each
survey. For example, consider the redshift range $0.60<z<0.75$. For each
redshift bin in this range, we have information contributing to three
parameters: $\alpha^{\mathrm{BOSS}}_{\perp}$, $\alpha_{||}^{\mathrm{BOSS}}$
and $\alpha^{\mathrm{eBOSS,LRG}}$. We first compute the $2\times 2$ Fisher
“submatrix” for $\alpha^{\mathrm{BOSS}}_{\perp}$,
$\alpha_{||}^{\mathrm{BOSS}}$ from the nonoverlapping sky area
$\Omega^{\mathrm{BOSS}}-\Omega^{\mathrm{eBOSS,LRG}}$ and add this to the full
matrix, then add on the $3\times 3$ Fisher submatrix for all three parameters
from the sky area $\Omega^{\mathrm{eBOSS,LRG}}$.

The Fisher matrix element for parameters $\lambda_{i}$ and $\lambda_{j}$
measured from survey $A$ within a redshift bin is Tegmark (1997)

$F^{A}_{ij}=\int_{0}^{k_{\mathrm{max}}}k^{2}dk\int_{0}^{1}d\mu\frac{\partial
P^{A}_{g}(k,\mu)}{\partial\lambda_{i}}C^{-1}_{P^{A}_{g}}(k,\mu)\frac{\partial
P^{A}_{g}(k,\mu)}{\partial\lambda_{j}},$ (4)

where we fix $k_{\mathrm{max}}=0.3h\,\mathrm{Mpc^{-1}}$.
$C_{P^{A}_{g}}(k,\mu)$ is the covariance matrix of the galaxy power spectrum
$P^{A}_{g}(k,\mu)$,

$C_{P^{A}_{g}}(k,\mu)=\frac{4\pi^{2}}{V^{A}}\biggl{[}P^{A}_{g}(k,\mu)+\frac{V^{A}}{N^{A}}\biggl{]}^{2},$
(5)

where $N^{A}$ is the number of galaxies in the redshift bin, and $V^{A}$ is
the cosmological volume contained within the redshift bin and sky area
$\Omega^{A}$. The first term in the covariance matrix models cosmic variance;
the second is shot noise from the finite number of galaxies. We model the
galaxy power spectrum as a function of the matter power spectrum $P_{m}(k)$
and potentially scale-dependent growth rate of structure $f(k)$, splitting it
into two components using the smoothed matter power spectrum $P_{sm}(k)$,

$P^{A}_{g}(k,\mu)=(b^{A}+f(k)\mu^{2})^{2}P_{sm}(k)\\\
\times\biggl{[}1+\biggl{(}\frac{P_{m}(k)}{P_{sm}(k)}-1\biggl{)}e^{-\frac{1}{2}k^{2}[\mu^{2}\Sigma^{A}_{nl,||}+(1-\mu^{2})\Sigma^{A}_{nl,\perp}]}\biggl{]}.$
(6)

Here $b^{A}$ represents the linear galaxy bias for the survey, and
$\Sigma^{A}_{nl,||}$ and $\Sigma^{A}_{nl,\perp}$ account for nonlinear damping
of the BAO feature. The values for these are fixed to the best-fit values from
the BOSS and eBOSS analyses. We compute $P_{sm}(k)$ using the method of Ref.
Hinton  _et al._ (2017).

Figure 3: 1D and joint 2D posteriors for the correlation coefficients $c_{0}$
[$(D_{M}/r_{s})^{\mathrm{BOSS}}$ and $(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$],
$c_{1}$ [$(Hr_{s})^{\mathrm{BOSS}}$ and $(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$],
and $c_{2}$ [$(D_{V}/r_{s})^{\mathrm{eBOSS,LRG}}$ and
$(D_{V}/r_{s})^{\mathrm{eBOSS,QSO}}$]. Also shown are the variations of the
correlation coefficients with neutrino mass parameters $N_{\mathrm{eff}}$ and
the Hubble parameter $H_{0}$.

The split into smooth and nonsmooth components of the matter power spectrum
ensures that we are only including information from the BAO scale, and not the
broadband shape of the power spectrum or redshift-space distortions. As such,
the derivatives with respect to $\lambda_{i,j}$ in Eq. (4) are computed only
on the $P_{m}(k)/P_{sm}(k)$ component of $P^{A}_{g}(k,\mu)$. We do this by
finite differencing $P_{m}(k^{\prime})/P_{sm}(k^{\prime})$ evaluated at
$k^{\prime}=k/\alpha$ or

$k^{\prime}=\frac{k}{\alpha_{\perp}}\biggl{[}1+\mu^{2}\biggl{(}\frac{\alpha^{2}_{\perp}}{\alpha^{2}_{||}}-1\biggl{)}\biggl{]}^{1/2}.$
(7)

Overall, the calculation of the Fisher matrix for a particular redshift bin
and survey matches that commonly used in the literature Seo and Eisenstein
(2007); Font-Ribera  _et al._ (2014a) and our model power spectrum is
representative of how BAO constraints are actually extracted from the data
Alam  _et al._ (2017); Ata  _et al._ (2018).

In Figure 3, we show the distribution of the correlation coefficients in our
BAO scale joint likelihood from our main NH fit. The variation of the
correlation coefficients with cosmological parameters is small but
perceptible, with $c_{1}$ and $c_{2}$ both increasing along with $H_{0}$ and
$N_{\mathrm{eff}}$, and $c_{0}$ decreasing with larger $H_{0}$ and
$N_{\mathrm{eff}}$. The trends are weaker but in the opposite direction for
the neutrino mass parameters. In general $c_{1}$ and $c_{2}$ are strongly
correlated with each other and anticorrelated with $c_{0}$. Our method can be
easily extended to include other datasets, and we expect larger variation in
the values of the coefficients for models including nonzero curvature or non-
cosmological-constant models of dark energy.

##  Appendix B Impacts of priors on the mass of the lightest neutrino

Figure 4: 1D and joint 2D posteriors on the neutrino mass and relevant
cosmological parameters, assuming normal ordering. Here we reproduce the
results from Figure 1 based on a linear prior for the mass of the lightest
neutrino and compare them to the results when employing a hybrid prior, linear
below $m_{\nu_{0}}=0.0003$ eV and logarithmic above.

In the main body of this article, we presented only results based on a linear
prior for the lightest neutrino mass $m_{\nu_{0}}$. This is a conservative
choice, as it is an essentially uninformative prior (but see Ref. Heavens and
Sellentin (2018) for an objective Bayesian construction of the most
uninformative prior for this problem). Although more informative, a
logarithmic prior is arguably more physically justified: Given that we have no
information suggesting any preferred scale for $m_{\nu_{0}}$ (or, in fact any
preference for $m_{\nu_{0}}>0$ at all), one might argue that any scale for the
mass is as likely as another. In Figure 4 we compare our result for the NH
with the result if we instead adopt a logarithmic prior for $m_{\nu_{0}}$
above $0.0003$ eV (and retain the linear prior below this value, enforcing
continuity of the prior across the transition). Masses below 0.0003 eV are
indistinguishable from the massless case in the outputs of CLASS, so there is
little point in oversampling this region. As can be seen from Figure 4, this
hybrid linear-logarithmic prior indeed results in a much stronger preference
for very small neutrino masses, giving $m_{\nu_{0}}<0.020$ eV and $0.058<\sum
m_{\nu}/\mathrm{eV}<0.100$ at 95% confidence. The effect is similar in the IH,
producing $m_{\nu_{0}}<0.024$ eV and $0.098<\sum m_{\nu}/\mathrm{eV}<0.136$ at
95% confidence.

## References

  * Gando  _et al._ (2013) A. Gando _et al._ (KamLAND), Phys. Rev. D 88, 033001 (2013), arXiv:1303.4667 .
  * An  _et al._ (2017) F. P. An _et al._ (Daya Bay), Chin. Phys. C41, 013002 (2017), arXiv:1607.05378 .
  * Adey  _et al._ (2018) D. Adey _et al._ (Daya Bay), Phys. Rev. Lett. 121, 241805 (2018), arXiv:1809.02261 .
  * (4) A. C. Serra, Double Chooz Improved Multi-Detector Measurements. CERN EP colloquium, CERN, Switzerland, September 20, 2016.
  * Bak  _et al._ (2018) G. Bak _et al._ (RENO), Phys. Rev. Lett. 121, 201801 (2018), arXiv:1806.00248 .
  * Adamson  _et al._ (2013a) P. Adamson _et al._ (MINOS), Phys. Rev. Lett. 110, 251801 (2013a), arXiv:1304.6335 .
  * Adamson  _et al._ (2013b) P. Adamson _et al._ (MINOS), Phys. Rev. Lett. 110, 171801 (2013b), arXiv:1301.4581 .
  * (8) A. Izmaylov, T2K Neutrino Experiment: Recent Results and Plans. Flavour Physics Conference, Quy Nhon, Vietnam, August 13-19, 2017.
  * (9) M. Friend, Updated Results from the T2K Experiment with $3.13\times 10^{21}$ Protons on Target. KEK seminar, January 10, 2019.
  * Sanchez (2018) M. Sanchez, Nova results and prospects (2018), Neutrino 2018, Heidelberg, June 4, 2018.
  * Acero  _et al._ (2019) M. Acero _et al._ (NOvA), Phys. Rev. Lett. 123, 151803 (2019), arXiv:1906.04907 .
  * Cleveland  _et al._ (1998) B. T. Cleveland, _et al._ , Astrophys. J. 496, 505 (1998).
  * Kaether  _et al._ (2010) F. Kaether, W. Hampel, G. Heusser, J. Kiko, and T. Kirsten, Phys. Lett. B 685, 47 (2010), arXiv:1001.2731 .
  * Abdurashitov  _et al._ (2009) J. N. Abdurashitov _et al._ (SAGE), Phys. Rev. C80, 015807 (2009), arXiv:0901.2200 .
  * Aharmim  _et al._ (2013) B. Aharmim _et al._ (SNO), Phys. Rev. C88, 025501 (2013), arXiv:1109.0763 .
  * Hosaka  _et al._ (2006) J. Hosaka _et al._ (Super-Kamiokande), Phys. Rev. D 73, 112001 (2006), arXiv:hep-ex/0508053 .
  * Cravens  _et al._ (2008) J. P. Cravens _et al._ (Super-Kamiokande), Phys. Rev. D 78, 032002 (2008), arXiv:0803.4312 .
  * Abe  _et al._ (2011) K. Abe _et al._ (Super-Kamiokande), Phys. Rev. D 83, 052010 (2011), arXiv:1010.0118 .
  * Ikeda (2018) M. Ikeda, Solar neutrino measurements with super-kamiokande (2018), Neutrino 2018, Heidelberg, June 5, 2018\.
  * Bellini  _et al._ (2011) G. Bellini _et al._ , Phys. Rev. Lett. 107, 141302 (2011), arXiv:1104.1816 .
  * Bellini  _et al._ (2010) G. Bellini _et al._ (Borexino), Phys. Rev. D 82, 033006 (2010), arXiv:0808.2868 .
  * Bellini  _et al._ (2014) G. Bellini _et al._ (Borexino), Nature 512, 383 (2014).
  * Aartsen  _et al._ (2015) M. G. Aartsen _et al._ (IceCube), Phys. Rev. D 91, 072004 (2015), arXiv:1410.7227 .
  * Abe  _et al._ (2018) K. Abe _et al._ (Super-Kamiokande), Phys. Rev. D 97, 072001 (2018), arXiv:1710.09126 .
  * Wong (2011) Y. Y. Wong, Ann. Rev. Nucl. & Part. Sci. 61, 69 (2011), arXiv:1111.1436 .
  * Lesgourgues and Pastor (2014) J. Lesgourgues and S. Pastor, New J. Phys. 16, 065002 (2014), arXiv:1404.1740 .
  * Vagnozzi  _et al._ (2017) S. Vagnozzi, _et al._ , Phys. Rev. D 96, 123503 (2017), arXiv:1701.08172 .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020a), arXiv:1807.06209 .
  * Loureiro  _et al._ (2019) A. Loureiro _et al._ , Phys. Rev. Lett. 123, 081301 (2019), arXiv:1811.02578 .
  * Ivanov  _et al._ (2020) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, Phys. Rev. D 101, 083504 (2020), arXiv:1912.08208 .
  * Archidiacono  _et al._ (2020) M. Archidiacono, S. Hannestad, and J. Lesgourgues, JCAP 09, 021 (2020), arXiv:2003.03354 .
  * Renk  _et al._ (2021) J. J. Renk, _et al._ (GAMBIT Cosmology Workgroup), JCAP 02, 022 (2021), arXiv:2009.03286 .
  * Athron  _et al._ (2017) P. Athron, _et al._ (GAMBIT Collaboration), Eur. Phys. J. C 77, 784 (2017), addendum in Athron  _et al._ (2018b), arXiv:1705.07908 .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A5 (2020b), arXiv:1907.12875 .
  * Esteban  _et al._ (2019) I. Esteban, M. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz, JHEP 01, 106 (2019), v4.1 www.nu-fit.org, arXiv:1811.05487 .
  * Beutler  _et al._ (2011) F. Beutler, _et al._ , MNRAS 416, 3017 (2011), arXiv:1106.3366 .
  * Ross  _et al._ (2015) A. J. Ross, _et al._ , MNRAS 449, 835 (2015), arXiv:1409.3242 .
  * Alam  _et al._ (2017) S. Alam _et al._ (BOSS), MNRAS 470, 2617 (2017), arXiv:1607.03155 .
  * Ata  _et al._ (2018) M. Ata _et al._ , Mon. Not. Roy. Astron. Soc. 473, 4773 (2018), arXiv:1705.06373 .
  * Bautista  _et al._ (2018) J. E. Bautista, _et al._ , ApJ 863, 110 (2018), arXiv:1712.08064 .
  * Blomqvist  _et al._ (2019) M. Blomqvist, _et al._ , A&A 629, A86 (2019), arXiv:1904.03430 .
  * Abbott  _et al._ (2019) T. M. C. Abbott, _et al._ , MNRAS 483, 4866 (2019), arXiv:1712.06209 .
  * Riess  _et al._ (2019) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Astrophys. J. 876, 85 (2019), arXiv:1903.07603 .
  * Wong  _et al._ (2020) K. C. Wong, _et al._ , MNRAS 498, 1420 (2020), arXiv:1907.04869 .
  * Freedman  _et al._ (2019) W. L. Freedman, _et al._ , ApJ 882, 34 (2019), arXiv:1907.05922 .
  * Chrzaszcz  _et al._ (2020) M. Chrzaszcz, _et al._ , Eur. Phys. J. C 80, 569 (2020), arXiv:1908.02302 .
  * Athron  _et al._ (2018a) P. Athron, _et al._ (GAMBIT Models Workgroup), Eur. Phys. J. C 78, 22 (2018a), arXiv:1705.07936 .
  * Esteban  _et al._ (2020) I. Esteban, M. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178 (2020), arXiv:2007.14792 .
  * Tanabashi  _et al._ (date) M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018 and 2019 update).
  * Cooke  _et al._ (2018) R. J. Cooke, M. Pettini, and C. C. Steidel, Astrophys. J. 855, 102 (2018), arXiv:1710.11129 .
  * Scolnic  _et al._ (2018) D. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 .
  * Carter  _et al._ (2018) P. Carter, _et al._ , Monthly Notices of the Royal Astronomical Society 481, 2371–2383 (2018).
  * Seo and Eisenstein (2007) H. Seo and D. J. Eisenstein, The Astrophysical Journal 665, 14 (2007).
  * Font-Ribera  _et al._ (2014a) A. Font-Ribera, _et al._ , Journal of Cosmology and Astroparticle Physics 2014 (05), 023–023.
  * Audren  _et al._ (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001 (2013), arXiv:1210.7183 .
  * Brinckmann and Lesgourgues (2019) T. Brinckmann and J. Lesgourgues, Phys. Dark Univ. 24, 100260 (2019), arXiv:1804.07261 .
  * Arbey (2012) A. Arbey, Comput. Phys. Commun. 183, 1822 (2012), arXiv:1106.1363 .
  * Arbey  _et al._ (2020) A. Arbey, J. Auffinger, K. P. Hickerson, and E. S. Jenssen, Comput. Phys. Commun. 248, 106982 (2020), arXiv:1806.11095 .
  * Blas  _et al._ (2011) D. Blas, J. Lesgourgues, and T. Tram, JCAP 1107, 034 (2011), arXiv:1104.2933 .
  * de Salas and Pastor (2016) P. F. de Salas and S. Pastor, JCAP 1607 (07), 051, arXiv:1606.06986 .
  * Akita and Yamaguchi (2020) K. Akita and M. Yamaguchi, JCAP 08, 012, arXiv:2005.07047 .
  * Froustey  _et al._ (2020) J. Froustey, C. Pitrou, and M. C. Volpe, JCAP 12, 015 (2020), arXiv:2008.01074 .
  * Zyla  _et al._ (2020) P. A. Zyla _et al._ (Particle Data Group), Prog. Theor. Exp. Phys. 083, C01 (2020).
  * Handley  _et al._ (2015) W. J. Handley, M. P. Hobson, and A. N. Lasenby, MNRAS 453, 4384 (2015), arXiv:1506.00171 .
  * Martinez  _et al._ (2017) G. D. Martinez, _et al._ (GAMBIT Scanner Workgroup), Eur. Phys. J. C 77, 761 (2017), arXiv:1705.07959 .
  * Alam  _et al._ (2021) S. Alam _et al._ (eBOSS), Phys. Rev. D 103, 083533 (2021), arXiv:2007.08991 .
  * de Salas  _et al._ (2021) P. de Salas, _et al._ , JHEP 02, 071 (2021), arXiv:2006.11237 .
  * Aker  _et al._ (2019) M. Aker _et al._ (KATRIN), Phys. Rev. Lett. 123, 221802 (2019), arXiv:1909.06048 .
  * Font-Ribera  _et al._ (2014b) A. Font-Ribera _et al._ (BOSS), JCAP 1405, 027 (2014b), arXiv:1311.1767 .
  * (70) GAMBIT Cosmology Workgroup, Supplementary Data: Strengthening the bound on the mass of the lightest neutrino with terrestrial and cosmological experiments, (2020), https://doi.org/10.5281/zenodo.4005381.
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * Lewis (2019) A. Lewis, (2019), arXiv:1910.13970 .
  * Scott (2012) P. Scott, Eur. Phys. J. Plus 127, 138 (2012), arXiv:1206.2245 .
  * Carron (2013) J. Carron, Astronomy & Astrophysics 551, A88 (2013).
  * Howlett  _et al._ (2017) C. Howlett, L. Staveley-Smith, and C. Blake, Monthly Notices of the Royal Astronomical Society 464, 2517–2544 (2017).
  * Tegmark (1997) M. Tegmark, Physical Review Letters 79, 3806–3809 (1997).
  * Hinton  _et al._ (2017) S. R. Hinton, _et al._ , Monthly Notices of the Royal Astronomical Society 464, 4807–4822 (2017).
  * Heavens and Sellentin (2018) A. F. Heavens and E. Sellentin, JCAP 04, 047, arXiv:1802.09450 .
  * Athron  _et al._ (2018b) P. Athron, _et al._ (GAMBIT Collaboration), Eur. Phys. J. C 78, 98 (2018b), addendum to Athron  _et al._ (2017), arXiv:1705.07908 .