# Light Sterile Neutrinos in the Early Universe: Effects of Altered Dispersion
Relations and a coupling to Axion-Like Dark Matter

Dominik Hellmann1, Heinrich Päs1
1 Fakultät für Physik, Technische Universität Dortmund, Germany

###### Abstract

We investigate the cosmological consequences of light sterile neutrinos with
altered dispersion relations (ADRs) and couplings to an ultra-light, axion-
like scalar field. In particular we study the impact on the number of
additional, light, fermionic degrees of freedom and primordial
nucleosynthesis. While the ADR leads to a new potential term in the
Hamiltonian, the coupling to the scalar field results in a time dependent,
effective mass contribution. We solve the quantum kinetic equations (QKEs) for
the neutrino density matrix and find that in certain parameter regions both
new physics effects can individually yield a suppressed population of sterile
neutrino species and the correct observed amount of helium in nucleosynthesis.
Combining both effects opens up new patches of parameter space excluded by
experimental bounds applying to models featuring only one of the effects.

##  1 Introduction

There exist several anomalies in short baseline (SBL) neutrino experiments
that point towards one or more sterile neutrinos in the
$\mathcal{O}(1\,\mathrm{eV})$ range [1, 2, 3, 4]. Even after the recently
published MicroBooNE data [5, 6, 7, 8] had disfavored the excess observed in
the MiniBooNE experiment, sterile neutrinos in this mass range remain an
interesting option for Standard Model (SM) extensions [9, 10].

As has been shown many years ago, the simple solution of just adding a few
sterile neutrino generations in the correct mass range to the SM is not
viable, since it solves the problem only for SBL experiments while long
baseline (LBL) and atmospheric neutrino oscillation experiments do not observe
any anomalies [11, 12, 13, 14, 15, 16] at the same baseline-energy ratio
$L/E$.

Therefore, if the SBL anomalies are indeed explained by light sterile
neutrinos, the oscillation probability cannot depend on $L/E$ in the same way
as it would be the case in the standard description. There exist several
proposals to resolve this tension (for a recent review see [10]). One such
possibility is to include an additional potential term for sterile neutrinos
in the Hamiltonian that parametrizes possible new physics effects. As a
consequence, an energy dependent mixing pattern arises between active and
sterile neutrinos. In this work, we mainly concentrate on a model proposed and
investigated in refs. [17, 18, 19, 20, 21, 22, 23] where the active-sterile
mixing becomes maximal at some resonance energy $E_{\mathrm{res}}$ and nearly
vanishes at higher energies. Such altered dispersion relations (ADRs) for
example arise if sterile neutrinos can take shortcuts through asymmetrically
warped extra dimensions that are unaccessible to particles charged under the
SM gauge group [24]. In order to solve the SBL anomalies, one has to require
that $E_{\mathrm{res}}$ is located within the energy ranges probed by these
experiments such that oscillations incorporating the new heavier mass
eigenstate associated with the sterile neutrino can reproduce the observed
pattern. At higher energies, the sterile neutrinos decouple from the active
ones due to mixing suppression and the standard oscillation pattern is
restored in the energy ranges typically probed by LBL experiments.

While this modification may help to resolve the tensions encountered between
SBL and LBL experiments, there remains a severe tension with cosmological
observations. For example, the effective number of neutrino generations
$N_{\mathrm{eff}}:=3+\Delta N_{\mathrm{eff}}$ has been shown to be very close
to 3 at the time of recombination according to recent Planck data [25], while
it could be altered dramatically if new, light, sterile degrees of freedom are
equilibrated in the early universe. Furthermore, the helium fraction produced
in big bang nucleosynthesis (BBN), $Y_{\mathrm{{}^{4}He}}$, is sensitive to
two effects. The first and most significant one is the faster expansion of the
universe caused by the presence of additional light degrees of freedom leading
to an earlier freeze out of neutron-proton reactions. The second effect is the
possible conversion of electron neutrinos to sterile neutrinos also resulting
in a potentially earlier freeze out since electron neutrinos are a substantial
part of almost all reactions keeping neutrons and protons in equilibrium with
the plasma. Sterile neutrinos with ADRs may also ameliorate this cosmological
tension by suppressing the population of the sterile flavor at high energies
in the early universe, as has been suggested in [20].

Another idea that has been proposed to reconcile light sterile neutrinos with
cosmological data employs an axion-like particle (ALP) coupling to the sterile
neutrinos via a Yukawa interaction term [26, 27, 28, 29, 30, 31] (for non-
cosmological applications of this particular class of models see also [32, 33,
34, 35, 36]). The additional, ultra-light scalar field $\phi$ is assumed to be
an approximately homogeneous condensate behaving like a classical field. Its
time evolution is governed by the Klein-Gordon equation in an expanding
universe and influences the effective mass of the sterile neutrino and hence
the mixing of the sterile and active neutrino.

Note that similar ideas pursued in the literature [38, 39, 40, 41, 42, 43, 37]
also involve the coupling of sterile neutrinos to dark scalars or massive
gauge bosons inducing an effective potential in the oscillation Hamiltonian
suppressing the thermalization of sterile degrees of freedom in the early
universe. In the models considered in this work, the $\nu_{s}$-$\phi$ coupling
is small and consequently, the induced matter potential and non-linear self
interactions of $\nu_{s}$ are strongly suppressed. The additional sterile
neutrino mass still remains relevant at early times due to the large initial
amplitude of the classical part of $\phi$. Furthermore, the ADR potential
included in the oscillation Hamiltonian is not assumed to be induced by
interactions as described above and can be sufficiently weakly temperature
dependent such that this dependence can be neglected.

In this paper we discuss both, the individual cosmological effects of and the
interplay between altered dispersion relations and a time varying sterile mass
induced by a light scalar field. In order to estimate the influence of the
full model on the chosen set of cosmological quantities, i.e. the effective
number of additional neutrino generations $\Delta N_{\mathrm{eff}}$ and the
helium fraction $Y_{\mathrm{{}^{4}He}}$, we solve the quantum kinetic
equations (QKEs) for the neutrino density matrix $\varrho$. This defines the
main task of this work, since from $\varrho$ we can derive the neutrino phase
space distributions and energy densities. In order to simplify our framework
for this study, we only consider a two flavor system of one active and one
sterile flavor. Moreover, we take the active neutrino to be the electron
neutrino to be able to draw realistic conclusions about nucleosynthesis.

This paper is organized as follows: in section 2, we describe the model in
detail and define the parameter space and crucial parameters. Section 3 is
concerned with setting up and numerically solving the Boltzmann equations.
Furthermore, we discuss the proper definition of the density matrix. In
section 4, we present results for the final helium abundances and $\Delta
N_{\mathrm{eff}}$ for various benchmark points in the model parameter space.
Finally in section 5, we draw our conclusions.

##  2 Sterile Neutrinos with Altered Dispersion Relations coupling to Axion-
Like Dark Matter

###  2.1 Sterile Neutrinos with Altered Dispersion Relations

Assuming that the relativistic dispersion relations of neutrinos are altered
by some unspecified new physics effect such as the presence of asymmetrically
warped extra dimensions111The asymmetric warping leads to the effect that
sterile neutrinos propagating through the extra dimension between two points
$x_{1}$ and $x_{2}$ on the SM brane need less time to complete their journey
than SM neutrinos traveling between those points on geodesics bound to the SM
brane [17, 24]. An observer located on the brane will hence come to the
conclusion that the usual energy momentum relation does not hold for sterile
neutrinos and needs to incorporate an effective potential into
$E^{2}=p^{2}+m^{2}$ in order to describe their behavior using brane bound
quantities. leads to new terms in the propagation Hamiltonian. Employing the
usual ultra-relativistic expansion of the neutrino dispersion relation yields

$\displaystyle\mathcal{H}(p)$
$\displaystyle=\frac{1}{2p}M^{\dagger}M+V_{s}(p)\,,$ (1) $\displaystyle
V_{s}(p)$ $\displaystyle=-\frac{bp}{2}\mathbb{P}_{s}\,,$ (2)

where $p$ is the average neutrino momentum, $b$ is the so called ADR parameter
controlling the strength of the ADR effect, $\mathbb{P}_{s}$ is the sterile
neutrino projector and $M$ is the neutrino mass matrix.

In the flavor basis, the last two quantities read

$\displaystyle M$ $\displaystyle=\begin{pmatrix}m_{ee}&m_{es}\\\
m_{es}^{\ast}&m_{ss}\\\ \end{pmatrix}\,,$ (3) $\displaystyle\mathbb{P}_{s}$
$\displaystyle=\begin{pmatrix}0&0\\\ 0&1\\\ \end{pmatrix}\,,$ (4)

for the $2\times 2$ neutrino system under consideration. Furthermore, we
choose all mass parameters to be real valued meaning that there is no CP
violation in the active-sterile mixing. To be specific, we use

$\displaystyle m_{ee}\approx 0\,\mathrm{eV}\,,\quad m_{es}\approx
0.12\,\mathrm{eV}\,,\quad m_{ss}\approx 1.1\,\mathrm{eV}\,,$ (5)

in accordance with ref. [29] and fits to SBL data.

Since we consider the system in the early universe, we also have to include a
potential for the electron neutrino induced by elastic scattering processes
modifying the effective masses of the neutrino matter eigenstates. Therefore,
the full Hamiltonian reads [44]

$\displaystyle\mathcal{H}(p)$
$\displaystyle=\frac{1}{2p}M^{\dagger}M+V_{e}(p)+V_{s}(p)$ (6) $\displaystyle
V_{e}(p)$
$\displaystyle=-\frac{8\sqrt{2}G_{F}p}{3}\left(\frac{\rho_{e^{\pm}}}{m_{W}^{2}}+\frac{\rho_{\nu_{e}}}{m_{Z}^{2}}\right)\mathbb{P}_{e}\,,$
(7)

with the Fermi coupling $G_{F}$, the electron and neutrino energy densities
$\rho_{\alpha}$ and the $W,Z$-Boson masses $m_{W,Z}$. Here, we neglected the
usual MSW potential proportional to the particle-antiparticle asymmetries
because we assume them to be small (of the order of the baryon asymmetry).
Thus, $V_{e}$ only contains the more significant higher order elastic
scattering contributions from the inverse mass expansion of the $W,Z$-Boson
propagators.

In order to analyze the resonance structure of this two flavor system, we
diagonalize the Hamiltonian with the general ansatz

$\displaystyle U(\theta(p)):=\begin{pmatrix}\cos(\theta(p))&\sin(\theta(p))\\\
-\sin(\theta(p))&\cos(\theta(p))\\\ \end{pmatrix}$ (8)

and find the following relation for the mixing angle

$\displaystyle\tan(2\theta(p))$
$\displaystyle=\frac{2m_{es}(m_{ee}+m_{ss})}{(m_{ss}^{2}-m_{ee}^{2})+2p(V_{s}(p)-V_{e}(p))}\,.$
(9)

The mixing becomes maximal as soon as $\theta=\pi/4$ because then electron and
sterile neutrinos equally constitute both mass eigenstates and the mass gap
between these two eigenstates is minimal. This in turn leads to higher
transition rates between active and sterile neutrinos in the energy regime
close to the resonance.

From eq. (9) we can infer the resonance condition

$\displaystyle\theta\rightarrow\frac{\pi}{4}$
$\displaystyle\Rightarrow\tan(2\theta(p))\rightarrow\pm\infty$
$\displaystyle\Rightarrow(m_{ss}^{2}-m_{ee}^{2})+2p(V_{s}(p)-V_{e}(p))\rightarrow
0^{\pm}\,.$ (10)

The momentum $p_{\mathrm{res}}$ fulfilling the condition (10) is called the
resonance momentum. As follows from eq. (10), the interplay of the two
appearing potentials and the neutrino masses determines the resonance
structure of the system.

This will become important again as soon as we study the effects of resonant
conversion between sterile and active neutrinos in combination with collisions
in the early universe plasma.

###  2.2 Sterile Neutrinos coupling to Axion-Like Dark Matter

In addition to the ADR effects, we also introduce an ultra-light, real, scalar
field. The corresponding scalar particles are assumed to be produced non-
thermally and to form a coherent condensate. This condensate behaves like a
homogeneous, classical field in an expanding universe and its time evolution
is given by [29]

$\displaystyle\phi(t)$ $\displaystyle=\phi_{0}\eta(t)$ (11)
$\displaystyle\eta(t)$
$\displaystyle=1.08\frac{J_{\frac{1}{4}}(m_{\phi}t)}{\sqrt[4]{m_{\phi}t}}\quad\mathrm{with}\quad\lim\limits_{t\rightarrow
0}\eta(t)=1\,.$ (12)

Here, $J_{\frac{1}{4}}$ is the regular Bessel $J$ function of fractional order
$1/4$, $m_{\phi}$ is the mass of the scalar and $t$ is the cosmic time.
Moreover, we allow for a feeble222The coupling has to be feeble in order for
the scalar field not to thermalize if the sterile neutrinos are thermalized
via oscillations. coupling to sterile neutrinos via a Yukawa interaction term
in its Lagrangian,

$\displaystyle\mathcal{L}(\phi,\partial_{\mu}\phi)$
$\displaystyle=\frac{1}{2}(\partial_{\mu}\phi)(\partial^{\mu}\phi)-\frac{m_{\phi}^{2}}{2}\phi^{2}-\frac{\lambda}{2}\phi\bar{\nu}_{s}\nu_{s}\,.$
(13)

Since we assume the coupling constant $\lambda$ to be very small, we neglect
possible interactions between $\nu_{s}$ and $\phi$ quanta and interactions of
$\nu_{s}$ with itself mediated by $\phi$. Nevertheless, we account for an
additional time dependent sterile neutrino mass term

$\displaystyle m_{\mathrm{eff}}(t)$
$\displaystyle=m_{ss}+\lambda\phi_{0}\eta(t)\,,$ (14)

modifying the mixing between electron and sterile neutrinos. Although this
additional mass term is proportional to $\lambda$, it is indeed significant if
the amplitude $\phi_{0}$ of the classical component of $\phi$ is sufficiently
large. Assuming $\phi$ to contribute a substantial amount to the dark matter
density, one can deduce [29]

$\displaystyle\phi_{0}\sim
10^{15}\,\mathrm{GeV}\left(\frac{10^{-15}\,\mathrm{eV}}{m_{\phi}}\right)^{\frac{1}{4}}\,.$
(15)

Thus, for scalar masses on the order of $m_{\phi}\sim 10^{-10}\,\mathrm{eV}$
one can still obtain $\lambda\cdot\phi_{0}\sim 100\,\mathrm{eV}$ for
$\lambda\sim 10^{-21}$.

At early times, the mass contribution remains approximately constant leading
to a constant sterile neutrino mass $m_{\mathrm{eff}}\approx
m_{ss}+\lambda\phi_{0}$. If $\phi_{0}$ is sufficiently large, the effective
mixing of electron and sterile neutrinos is negligible since $\tan(2\theta)$
is suppressed by the large mass gap $m_{\mathrm{eff}}^{2}-m_{ee}^{2}$ in the
denominator of eq. (9) with $m_{s}s$ replaced by $m_{\mathrm{eff}}$. Hence,
sterile neutrinos are not populated via oscillations at these early times.

From the point in time where the Bessel function starts its damped oscillating
behavior ($t\sim 1/m_{\phi}$) the effective sterile neutrino mass approaches
$m_{ss}$, allowing for significant oscillations again. Therefore, the neutrino
oscillation behavior depends on the mass parameter $m_{\phi}$, i.e. the
smaller $m_{\phi}$ the longer $\phi(t)$ will remain approximately constant and
active-sterile mixing will be suppressed by $m_{s0}$.

###  2.3 Central Quantities and Parameter Space

In the following analysis, we focus on two major quantities:

  1. 1.

The effective number of additional neutrino generations

$\displaystyle\Delta N_{\mathrm{eff}}(t):$
$\displaystyle=\frac{8}{7}\left(\frac{11}{4}\right)^{\frac{4}{3}}\sum_{k\in\mathcal{V}}\frac{\rho_{k}(t)}{\rho_{\gamma}(t)}-3\,,\quad\mathrm{with}\quad\mathcal{V}=\\{e,\mu,\tau,s\\}\,.$
(16)

Here, $\rho_{k}(t)$ are the neutrino energy densities at time $t$ and
$\rho_{\gamma}(t)$ is the corresponding photon density333By using the energy
density in the definition of $\Delta N_{\mathrm{eff}}$ we slightly differ from
the methods employed in ref. [29] where the number density is used instead..
Moreover, we need to incorporate factors of $8/7$ and $(11/4)^{4/3}$ to
directly compare fermionic and bosonic energy densities with a temperature
deviation of $(11/4)^{1/3}$ to each other.

  2. 2.

The helium mass fraction

$\displaystyle Y_{{}^{4}\mathrm{He}}$
$\displaystyle:=\frac{4n_{\mathrm{He}}}{n_{\mathrm{B}}}\,,$ (17)

with the helium and baryon number densities $n_{\mathrm{He}}$,
$n_{\mathrm{B}}$.

In order to compare these quantities with observations, these observables have
to be computed at different times in the cosmic evolution. The value of the
number of additional neutrino generations $\Delta N_{\mathrm{eff}}(t)$ is
inferred from the Hubble rate measurement from the CMB [25] and thus needs to
be known at the time of the last photon scattering, while the value of the
helium fraction $Y_{{}^{4}\mathrm{He}}$ has to be computed shortly after BBN.
It is, however, sufficient to evaluate $\Delta N_{\mathrm{eff}}$ directly
after $e^{\pm}$ annihilation since from there on it remains constant. This is
because after $e^{\pm}$ annihilation the total neutrino energy density only
changes due to the expansion of the universe444We prove that the total
neutrino energy density only changes due to the expansion of the universe
after neutrino decoupling in Appendix A using the neutrino density matrix
formalism introduced in the next section. and by normalizing it to the photon
energy density we cancel the dependence on the scale factor.

The helium fraction remains constant right after BBN, hence it is appropriate
to evaluate $Y_{4^{\mathrm{He}}}$ as soon as deuterium dissociation has ceased
to be efficient. Hence, our analysis solely focuses on the era of radiation
domination.

Finally, we want to discuss the three dimensional parameter space of the model
under consideration. It is parameterized by

  1. 1.

The scalar field mass $m_{\phi}$ determining when the scalar field starts to
oscillate.

  2. 2.

The amplitude $m_{s0}:=\lambda\phi_{0}$ of the additional mass contribution
for the sterile neutrino.

  3. 3.

The ADR parameter $b$.

In the following, we assume a range for
$m_{s0}\in[10\,\mathrm{eV},250\,\mathrm{eV}]$ in accordance with ref. [29].
Furthermore, for $m_{\phi}$ we choose the interval of possible values to be
$[10^{-22}\,\mathrm{eV},10^{-14}\,\mathrm{eV}]$ since for $m_{\phi}\leq
10^{-22}\,\mathrm{eV}$ the scalar field starts oscillating so late that the
effective sterile neutrino mass remains constant during the considered
temperatures. Thus, for $m_{\phi}\lesssim 10^{-22}\,\mathrm{eV}$ our scenario
becomes independent of $m_{\phi}$ and yields the same constraints on $m_{s0}$
as for $m_{\phi}\sim 10^{-22}\,\mathrm{eV}$.

At values larger than $m_{\phi}=10^{-14}\,\mathrm{eV}$ the addition of the
scalar field to the model becomes meaningless since the sterile neutrino mass
already gets close to its bare value at the relevant temperatures and the
sterile species equilibrates.

For the ADR parameter, we choose benchmark values in $I_{b}=[0,\infty)$, i.e.
we consider anything between zero ADR potential and an arbitrarily large ADR
effect. A special point in this parameter range is $b\sim 10^{-17}$, since
this is the order of magnitude needed in order to explain SBL anomalies [21].
The reason why we allow for arbitrarily high (low) ADR parameters in the early
universe is that there might be some mechanism in the extradimensional
realization of ADRs causing the curvature of the extra dimension to change
from early times until today and correspondingly implying the ADR parameter to
decrease (increase), accordingly.

##  3 Neutrino Quantum Kinetic Equations and Numerical Strategy

The density matrix describing the oscillations of neutrinos in the early
universe is defined as the thermal average of creation and annihilation
operators, $a_{j}^{(\dagger)}(\vec{p})$ of neutrino mass eigenstates, i.e.
[44, 45]

$\displaystyle(2\pi)^{3}\delta^{3}(\vec{p}-\vec{q})\varrho_{jk}(\vec{p}):=\langle
a_{k}(\vec{p})^{\dagger}a_{j}(\vec{q})\rangle$ (18)

where the thermal average of an operator $\hat{O}$ is defined using the
density operator555Despite their very similar names the density operator and
density matrix are by no means equivalent quantities. $\hat{\Pi}$ of the
thermal system as usual
$\langle\hat{O}\rangle:=\mathrm{Tr}(\hat{\Pi}\hat{O})$.

Thus, if we consider the diagonal elements ($j=k$) of the density matrix it
simplifies to the thermal average of the occupation number operator of
$\nu_{k}$ which in turn results in its phase space distribution function. This
already gives some intuition of its physical meaning: The diagonal of the
density matrix contains the information how many $\nu_{k}$ momentum states on
average are occupied in the system. Moreover, by inspecting the off-diagonal
elements of $\varrho$, we obtain information about the average correlation
between $\nu_{j}(p)$ and $\nu_{k}(p)$. Such correlations arise for example due
to neutrino oscillations. This makes the density matrix the quantity of choice
if we want to consider incoherent particle collisions and oscillations in a
thermal environment [46] at the same time.

In the following, we will mainly work in the flavor basis instead of the mass
basis because the collision terms and Hamiltonian potentials are easier to
calculate in the flavor basis. Therefore, we need to transform $\varrho_{jk}$
into this basis using the neutrino mixing matrix $U$ from eq. (8) via

$\displaystyle\varrho_{jk}^{f}=\sum_{l,m=1}^{2}U_{jl}\varrho_{lm}(U^{\dagger})_{mk}\,.$
(19)

From now on, we will work with $\varrho^{f}$ only and drop the superscript
$f$.

###  3.1 Quantum Kinetic Equations for the Density Matrix

The time evolution of $\varrho$ in an expanding, homogeneous and isotropic
universe is governed by a Boltzmann-like, quantum kinetic equation (QKE) [44,
45, 47, 48, 49]

$\displaystyle(\partial_{t}-pH\partial_{p})\varrho(t,p)$
$\displaystyle=-i[\mathcal{H}(t,p),\varrho(t,p)]+\mathcal{C}[t,p,\varrho]\,,$
(20)

where $p$ is the modulus of the neutrino momentum, $t$ is the cosmic time, $H$
is the Hubble rate, $\mathcal{H}$ is the Hamiltonian from eq. (6) and
$\mathcal{C}$ is the collison operator.

The convectional derivative operator, $\partial_{t}-pH\partial_{p}$, on the
left hand side includes the effect of the expansion of the universe
redshifting the neutrino momentum as $p\propto a^{-1}$, where $a$ is the scale
factor of the Robertson Walker metric. Furthermore, the right hand side
contains the commutator of the neutrino Hamiltonian with the density matrix
and the collision operator $\mathcal{C}$. While the commutator part governs
the evolution of $\rho$ due to neutrino oscillations, the collision part
determines how many neutrinos are annihilated, created or scattered to other
momentum modes in interactions with the background plasma. Since there are
also neutrinos in the background plasma, this last term is non-linear.

The collision operator is the sum of the individual collision operators
corresponding to a scattering process involving neutrinos

$\displaystyle\mathcal{C}[t,p,\varrho]$
$\displaystyle=\sum_{k\in\mathcal{P}}\mathcal{C}_{k}[t,p,\varrho]\,,$ (21)

where the set of all processes $\mathcal{P}$ contains the interactions given
in table 1. Here we neglect neutrino-nucleon scattering processes due to
strong Boltzmann suppression of the nucleon distribution functions at
temperatures of $\mathcal{O}(100\,\mathrm{MeV})$ and smaller.

Table 1: All relevant processes considered in the neutrino collision terms. k | process
---|---
1 | $\nu e^{-}\leftrightarrow\nu e^{-}$
2 | $\nu e^{+}\leftrightarrow\nu e^{+}$
3 | $\nu\bar{\nu}\leftrightarrow e^{-}e^{+}$
4 | $\nu\nu\leftrightarrow\nu\nu$
5 | $\nu\bar{\nu}\leftrightarrow\nu\bar{\nu}$

For example, the collision term for the process $\nu e^{-}\leftrightarrow\nu
e^{-}$ is given by

$\displaystyle\mathcal{C}_{1}[t,p,\varrho]=$
$\displaystyle\frac{8G_{F}^{2}}{p}\int\mathrm{d}^{3}\vec{\pi}_{1}\,\mathrm{d}^{3}\vec{\pi}_{2}\,\mathrm{d}^{3}\vec{\pi}_{3}\,(2\pi)^{4}\delta^{4}(p^{\mu}+p_{1}^{\mu}-p_{2}^{\mu}-p_{3}^{\mu})$
$\displaystyle\times\left\\{4g_{l}^{2}(p_{\alpha}p_{3}^{\alpha})(p_{1\beta}p_{2}^{\beta})+4g_{R}^{2}(p_{\alpha}p_{1}^{\alpha})(p_{2\beta}p_{3}^{\beta})-4g_{L}g_{R}(p_{\alpha}p_{2}^{\alpha})m_{e}^{2}\right\\}$
$\displaystyle\times\phi_{\nu e^{-}\nu
e^{-}}(t,p^{\mu},p_{1}^{\mu},p_{2}^{\mu},p_{3}^{\mu})\,,$ (22)

where
$\mathrm{d}^{3}\vec{\pi}_{j}:=(2E_{j}(2\pi)^{3})^{-1}\mathrm{d}^{3}\vec{p}_{j}$
denotes the Lorentz invariant phase space measure,
$g_{L}=1/2+\sin^{2}(\theta_{W})$, $g_{R}=\sin^{2}(\theta_{W})$, $\theta_{W}$
is the Weinberg angle and the statistical factor $\phi_{\nu e^{-}\nu e^{-}}$
is given by

$\displaystyle\phi_{\nu e^{-}\nu
e^{-}}(t,p^{\mu},p_{1}^{\mu},p_{2}^{\mu},p_{3}^{\mu}):=\phi_{\nu e^{-}\nu
e^{-}}^{+}(t,p^{\mu},p_{1}^{\mu},p_{2}^{\mu},p_{3}^{\mu})-\phi_{\nu e^{-}\nu
e^{-}}^{-}(t,p^{\mu},p_{1}^{\mu},p_{2}^{\mu},p_{3}^{\mu})$ (23)

with

$\displaystyle\phi_{\nu e^{-}\nu
e^{-}}^{+}(t,p^{\mu},p_{1}^{\mu},p_{2}^{\mu},p_{3}^{\mu})$
$\displaystyle:=[1-f_{e^{-}}(t,p_{1})]f_{e^{-}}(t,p_{3})\\{\mathbb{P}_{e}\varrho(t,p_{2})\mathbb{P}_{e},\mathbb{I}-\varrho(t,p)\\}\,,$
$\displaystyle\phi_{\nu e^{-}\nu
e^{-}}^{-}(t,p^{\mu},p_{1}^{\mu},p_{2}^{\mu},p_{3}^{\mu})$
$\displaystyle:=f_{e^{-}}(t,p_{1})[1-f_{e^{-}}(t,p_{3})]\\{\mathbb{P}_{e}[\mathbb{I}-\varrho(t,p_{2})]\mathbb{P}_{e},\varrho(t,p)\\}\,.$

Furthermore, $f_{e^{-}}$ denotes the electron phase space distribution. The
collision terms are calculated applying the methods described in [48] to the
current scenario.

In addition to the density matrix for neutrino states, in principle there is
also an analogous one for antineutrinos which has to be solved at the same
time. In the following, we assume that the lepton-antilepton asymmetry is of
the order of the baryon asymmetry and hence negligible compared to the total
phase space densities. This implies that the antineutrino density matrix
behaves the same as the neutrino density matrix and therefore we just have to
consider the QKE for neutrinos.

In order to keep track of the temperature $T_{\gamma}$ of the electromagnetic
plasma, we need to solve the continuity equation of the universe

$\displaystyle\dot{\rho}$ $\displaystyle=-3H(\rho+P)\,,$ (24)

where $\rho$ and $P$ are the total energy density and total pressure of all
radiation species, respectively. By substituting in the equilibrium
expressions for electrons and photons and assuming these particles to be in
thermal equilibrium666This assumption is valid due to the rapid
electromagnetic interactions between photons and electrons roughly until the
time of last scattering., this equation can be reformulated into a
differential equation for $T_{\gamma}$.

###  3.2 Numerical Solution of the Quantum Kinetic Equations

In order to prepare the numerical solution of the previously introduced QKE
(20), we define a new set of dimensionless variables

$\displaystyle x(t,p)$ $\displaystyle:=m_{0}a(t)\,,$ (25) $\displaystyle
y(t,p)$ $\displaystyle:=a(t)p\,,$ (26)

where we choose $m_{0}=1\,\mathrm{MeV}$. Therefore, $x$ represents the
dimensionless scale factor and $y$ is a dimensionless momentum variable not
being redshifted over time, since $p\propto a^{-1}$. Moreover, $x$ takes the
role of the reciprocal of the neutrino temperature which is equal to
$T_{\gamma}$ at early times but deviates from it after neutrino decoupling and
electron positron annihilation. Transformed to these new variables eq. (20)
assumes the form

$\displaystyle xH\partial_{x}\tilde{\varrho}(x,y)$
$\displaystyle=-i[\mathcal{H}(x,y),\tilde{\varrho}(x,y)]+\mathcal{C}[x,y,\tilde{\varrho}(x,y)]\,,$
(27)

with $\tilde{\varrho}$ being the density matrix expressed in the new set of
variables. From now on, we will only refer to this quantity and hence drop the
tilde, i.e. $\tilde{\varrho}\to\varrho$.

In order to integrate eq. (27), in principle we had to start at $x_{0}=0$ and
set $\varrho_{ik}(x_{0},y)\equiv 0$. But since the QKE described in the last
section are only valid after the strong phase transition, we have to find a
finite starting point $x_{0}$ matching all criteria of validity of our
equations of motion which are

  1. 1.

Active neutrinos are in thermal equilibrium with the electromagnetic plasma.

  2. 2.

Quarks and gluons are bound into hadrons.

  3. 3.

Contributions from processes involving muons are negligble.

Furthermore, we assume the sterile neutrino density and correlations between
active and sterile neutrinos to be negligible at $x_{0}$ such that our initial
condition for $\varrho$ is given by

$\displaystyle\varrho(x_{0},y)\approx\begin{pmatrix}(\exp(y)+1)^{-1}&0\\\
0&0\\\ \end{pmatrix}\,.$ (28)

We found that $x_{0}=0.01$, i.e. $T_{\gamma,0}=100\,\mathrm{MeV}$, fulfills
these criteria. For more discussions see appendix B.

We terminate the integration at $x_{1}$ which we require to fulfill the
following criteria:

  1. 1.

Neutrino interactions are completely frozen out

  2. 2.

All free neutrons are bound into light nuclei

  3. 3.

The neutrino distribution functions have reached their asymptotic values

  4. 4.

The relativistic approximation for the oscillation Hamiltonian and collisions
is valid

We found $x_{1}=50$ to be a suitable final point fulfilling these criteria
while still being located in radiation domination.

To solve eq. (27) in the interval $X:=[x_{0},x_{1}]$, we discretize the
momentum space $\Omega_{y}=[0,\infty)$ and integrate the resulting set of
ordinary differential equations. To do so, we choose $N_{y}$ equidistant
points between the minimal and maximal momentum values
$y_{\mathrm{min}},y_{\mathrm{max}}$ at which we cut off the distribution
function. Choosing a minimal value is necessary because the ultra-relativistic
approximation employed within the oscillation Hamiltonian is not valid for all
momentum values, especially not for $y=0$. Hence, the minimal $y$ value is
chosen to be $y_{\mathrm{min}}=10^{-4}$ in order to still yield a reasonably
good approximation for the neutrino energy density. Furthermore, the maximal
momentum value is chosen to be $y=20$ since the neutrino distribution at this
$y$ value fulfills

$\displaystyle f_{\nu}(x,y=20)$ $\displaystyle\leq
f_{\mathrm{eq}}(x,y=20)=(\exp(y)+1)^{-1}|_{y=20}\approx 2\cdot 10^{-9}\,,$
(29)

which is sufficiently close to zero. Therefore, the total, relative error
induced within the neutrino energy density needed to calculate our central
quantities is of the order

$\displaystyle\epsilon_{\mathrm{rel}}:=\frac{|\rho_{\nu}^{\mathrm{approx}}-\rho_{\nu}|}{\rho_{\nu}}\sim
10^{-6}\,.$

The discretized version of $\Omega_{y}$ then reads

$\displaystyle\tilde{\Omega}_{y}(N_{y}):=\left\\{y_{k}=y_{\mathrm{min}}+k\cdot\Delta
y\,\Big{|}\,k\in\\{0,\ldots,N_{y}-1\\},\Delta
y=\frac{y_{\mathrm{max}}-y_{\mathrm{min}}}{N_{y}-1}\right\\}\,.$

Hence, we arrive at $N_{y}$ coupled, ordinary, differential equations for the
density matrix values at the chosen momentum nodes plus one equation for the
photon temperature. Furthermore, decomposing the Hermitian density matrix into
its 4 independent, real components yields a total of $4N_{y}+1$ coupled
differential equations which need to be solved777For the numerical details see
appendix C.

###  3.3 Calculating the helium Abundance

In this section, we present how the ${}^{4}\mathrm{He}$ mass fraction
$Y_{{}^{4}\mathrm{He}}$ is estimated from the neutrino distribution functions.
Our explanations and notation closely follow the book [50] by Bernstein on
kinetic theory in an expanding universe.

At first, we define the neutron fraction, i.e.

$\displaystyle X_{n}(t):=\frac{n_{n}(t)}{n_{n}(t)+n_{p}(t)}\,,$ (30)

where $n_{n}$ and $n_{p}$ are the neutron and proton number densities,
respectively. This choice greatly simplifies the Boltzmann equation for
$n_{n}$ since the $a^{-3}$ dependence of the number densities cancel.
Furthermore, we find

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\\{a^{3}(t)(n_{n}(t)+n_{p}(t))\\}\equiv\mathrm{const.}\,,$
(31)

because the baryon number in a comoving volume, $N_{B}\approx
a^{3}(n_{n}+n_{p})$, is conserved within all relevant processes shortly before
neutron freeze out. These processes are

$\displaystyle n+e^{+}$ $\displaystyle\leftrightarrow p+\bar{\nu}_{e}$ (32)
$\displaystyle n+\nu_{e}$ $\displaystyle\leftrightarrow p+e^{-}$ (33)
$\displaystyle n$ $\displaystyle\leftrightarrow p+\bar{\nu}_{e}+e^{-}\,.$ (34)

The differential equation governing the evolution of $X_{n}$ reads [50]

$\displaystyle\frac{\mathrm{d}X_{n}(x)}{\mathrm{d}x}=\frac{\lambda_{pn}(x)}{xH}(1-X_{n}(x))-\frac{\lambda_{np}(x)}{xH}X_{n}(x)\,.$
(35)

in terms of $x=m_{0}a(t)$. Here $\lambda_{np}$ is the thermal interaction rate
of all processes converting neutrons to protons, while $\lambda_{pn}$ is that
of all processes converting protons to neutrons. They are given and discussed
in appendix D.

We solve eq. (35) from $x(T=5\,\mathrm{MeV})$ where neutrons and protons are
still in thermal equilibrium up until $x(T=0.07\,\mathrm{MeV})$ where
deuterium dissociation ceases to be efficient. During the solution of this
differential equation, we interpolate the neutrino distribution functions and
temperatures obtained on the grid of $(x,y)$ values using the methods
described in the previous subsections.

Finally in order to estimate the produced helium abundance, we convert the
neutron fraction into the helium mass fraction, i.e.

$\displaystyle Y_{{}^{4}\mathrm{He}}$
$\displaystyle=\frac{4n_{{}^{4}\mathrm{He}}}{n_{\mathrm{B}}}=\frac{2n_{n}}{n_{\mathrm{B}}}=2X_{n}\,,$
(36)

where we used the assumption that approximately all free neutrons are bound
into helium-4 nuclei at the end of BBN.

Of course this method is subject to several approximations especially since we
neglect the nuclear reaction rates, hence our estimate cannot be compared
directly to observations from [51]. Nevertheless, we can inspect if the
different models lead to a relative deviation from the expected value which is
on the order of experimental uncertainty or if it exceeds this uncertainty
significantly.

##  4 Predicted Effective Degrees of Freedom and helium Abundance

Now, we present the results for different benchmark points within the 3
dimensional parameter space. In the following, we first discuss our results
for the pure ADR scenario, the pure scalar field scenario and afterwards for
the combination of both effects. For each chosen benchmark point, we calculate
the resulting effective, additional number of degrees of freedom $\Delta
N_{\mathrm{eff}}$ and the estimated helium-4 abundance
$Y_{{}^{4}\mathrm{He}}$. These simulated values for $\Delta N_{\mathrm{eff}}$
are compared to bounds obtained by the Planck collaboration, i.e.

  * •

TT + lowE ($95\%$ CL): $N_{\mathrm{eff}}=3.00_{-0.53}^{+0.57}$ $\Rightarrow$
$\mathbf{\Delta N_{\mathrm{eff}}\leq 0.57}$,

  * •

TT, TE, EE + lowE ($95\%$ CL): $N_{\mathrm{eff}}=2.92_{-0.37}^{+0.36}$
$\Rightarrow$ $\mathbf{\Delta N_{\mathrm{eff}}\leq 0.28}$,

  * •

TT + lowE + lensing + BAO ($95\%$ CL): $N_{\mathrm{eff}}=3.11_{-0.43}^{+0.44}$
$\Rightarrow$ $\mathbf{\Delta N_{\mathrm{eff}}\leq 0.55}$,

  * •

TT, TE, EE + lowE + lensing + BAO ($95\%$ CL):
$N_{\mathrm{eff}}=2.99_{-0.33}^{+0.34}$ $\Rightarrow$ $\mathbf{\Delta
N_{\mathrm{eff}}\leq 0.33}$.

Here the abbreviations TT, TE, EE, lowE, lensing, BAO refer to different
measurement techniques / features of the CMB data (i.e. TT $\hat{=}$ intensity
(temperature) only, TE $\hat{=}$ temperature + curl free polarization data, EE
$\hat{=}$ curl free polarization data only, lowE $\hat{=}$ curl free
polarization data only at low multipole moments, lensing $\hat{=}$ grav.
lensing measurement, BAO $\hat{=}$ baryon acoustic oscillations). Afterwards,
we turn towards the helium abundance and compare its deviation for different
benchmark points from the expected SM value to the experimental uncertainties
on the helium mass fraction from ref. [51], i.e.
$\sigma_{{}^{4}\mathrm{He}}=0.004$.

###  4.1 $\Delta N_{\mathrm{eff}}$ and $Y_{{}^{4}\mathrm{He}}$ in the pure ADR
Scenario

The values of $\Delta N_{\mathrm{eff}}$ obtained after the full integration of
QKEs for different ADR parameters are shown in table 2.

Table 2: Estimated additional light degrees of freedom $\Delta N_{\mathrm{eff}}$ at $x=50$ for different ADR parameters $b$. $b$ | $0$ | $10^{-17}$ | $10^{-15}$ | $10^{-12}$ | $10^{-6}$ | $10^{-4}$ | $10^{-2}$
---|---|---|---|---|---|---|---
$\Delta N_{\mathrm{eff}}$ | $1.36$ | $1.36$ | $1.36$ | $1.38$ | $0.04$ | $0.04$ | $0.04$

Here, we see that the resulting $\Delta N_{\mathrm{eff}}$ values for the
smallest ADR parameters exceed all Planck bounds by far while bigger $b$
values lead to excellent agreement with all four bounds $\Delta
N_{\mathrm{eff}}<\\{0.28,0.33,0.55,0.57\\}$. Moreover, we can infer that for
small $b$ values the no-ADR scenario is resembled and sterile neutrinos are
equilibrated via oscillations. Turning $b$ up to much larger values ($b\gtrsim
10^{-6}$) leads to a decrease in $\Delta N_{\mathrm{eff}}$ and sterile
neutrinos are not close to equilibrium anymore. In figure 1, we show the final
sterile neutrino distributions compared to the equilibrium distribution for
four $b$ values differing by many orders of magnitude to emphasize this
statement.

Figure 1: Final ($x=50$) $\nu_{s}$ phase space densities (solid) as a function
of the comoving momentum $y$ for various ADR parameters
$b\in\\{0,10^{-17},10^{-12},10^{-6}\\}$ compared to the equilibrium density
(dashed).

The behavior described above can be explained by considering the resonance
structure of each parameter configuration. In figure 2, we show the resonance
momentum $y_{\mathrm{res}}$ for multiple ADR parameters in the temperature
range $T_{\nu}\in\mathcal{T}_{\nu}:=[3,100]\,\mathrm{MeV}$.

Figure 2: Resonance momentum $y_{\mathrm{res}}$ plotted against the neutrino
temperature $T_{\nu}\approx T_{\gamma}$ for various ADR parameters
$b\in\\{10^{-17},10^{-15},10^{-12},10^{-6},10^{-4},10^{-2}\\}$. The relevant
momentum region $10^{-4}\leq y\leq 20$ is shown as the green shaded region.
All modes above the corresponding resonance curve are subject to strong mixing
suppression.

For larger ADR parameters the resonance curve passes through the relevant $y$
region from above leading to resonantly enhanced $\nu_{e}$-$\nu_{s}$
conversion. Momentum modes located underneath the respective resonance curve
neither experience mixing enhancement nor mixing suppression and approximately
behave as in the no-ADR scenario. On the other hand, momentum modes well above
the resonance curve, i.e. $y\gg y_{\mathrm{res}}(T_{\nu})$ $\forall
T_{\nu}\in\mathcal{T}_{\nu}$, are subject to effective mixing suppression
since $\nu_{1}\approx\nu_{e}$ and $\nu_{2}\approx\nu_{s}$. Thus, $\nu_{s}$
remains unpopulated in this regime.

Note that the strength of this suppression/enhancement effect is momentum
dependent since the active as well as the sterile potentials are proportional
to $y$. Hence, the mixing of neutrinos with small momenta is closer to the
vacuum case leading to a faster population of these modes even if the
resonance momentum is much smaller. But since $\Delta N_{\mathrm{eff}}$
depends more strongly on the high momentum region, this does not affect its
estimate significantly.

Now one could ask what happens as soon as multiple (high) momentum modes pass
through one or two resonances, as is the case for e.g. $b=10^{-15}$. A sterile
$y$ mode can be strongly populated by passing through a resonance, but it is
even more important how big the long term mixing around this resonance is. If
the mixing is sufficiently smaller than vacuum mixing, the respective sterile
mode experiences some enhancement by the resonance, but it will not reach its
equilibrium value. On the other hand, if the mixing is very large already
before or after passing the resonance, the sterile momentum mode will be
equilibrated irregardless of the resonance. Thus what matters is only the fact
whether the relevant momentum modes are subject to mixing suppression for
sufficiently long or if they are closer to the (large) vacuum mixing. This can
be seen in figure 3 which shows the temperature evolution of the modulus
$|\varrho_{es}|:=\sqrt{\mathrm{Re}(\varrho_{es})^{2}+\mathrm{Im}(\varrho_{es})^{2}}$
for $y=5$.

Figure 3: Time evolution of the off-diagonal density matrix element
$|\varrho_{es}(y=5)|$ for several ADR parameters. The $y=5$ mode undergoes a
resonance at $T_{\nu}\approx 7\,\mathrm{MeV}$ for $b=10^{-15}$ leading to a
delayed increase of $|\varrho_{es}|$ afterwards. The temperature at which this
resonance occurs is marked by the orange dashed vertical line.

The off-diagonal element $\varrho_{es}$ is important since it contains
information about the energy transfer from $\nu_{e}$ to $\nu_{s}$. Therefore,
the plot shows that shortly after $T=\mathcal{O}(100\,\mathrm{MeV})$ for
models with vanishing or small ADR parameters a dip occurs in $\varrho_{es}$
leading to a significant enhancement of $\rho_{ss}$, afterwards. Considering,
the curve for $b=10^{-15}$, we see that the resonance around
$\mathcal{O}(7\,\mathrm{MeV})$ leads to a significant impact in
$|\varrho_{es}|$ but doesn’t lead to a significant increase of the sterile
neutrino density, since it already reached thermal equilibrium. Moreover, we
see that for $b\geq 10^{-6}$ the off-diagonal matrix elements stay much closer
to zero due to sizeable mixing suppression.

Despite the fact that the excess of light degrees of freedom, $\Delta
N_{\mathrm{eff}}$, is a good estimator for the degree of population of sterile
neutrinos, it is not sufficient to rely on this number alone. After neutrino
decoupling around $T_{\gamma}=\mathcal{O}(3\,\mathrm{MeV})$, the active-
sterile oscillations could lead to a depletion of the density of active
neutrinos which has strong impact on the neutron-proton equilibrium, since
fewer $\nu_{e}$ lead to an early freeze out of $n$-$p$ reactions and hence to
a larger neutron abundance. This leads to an excess of helium that contradicts
the very good agreement of the predictions from standard cosmology and
cosmological observations. Hence, we need to carefully estimate how much
helium is produced for our chosen parameter configurations. To do so, we now
estimate the impact of the $\nu_{e}$ depletion on nucleosynthesis by
proceeding as described in section 3.3 and solve the Boltzmann equation for
the neutron fraction $X_{n}=n_{n}/(n_{n}+n_{p})$. The neutron fraction
$X_{n}=n_{n}/(n_{n}+n_{p})$ can then be translated into the helium mass
fraction $Y_{{}^{4}\mathrm{He}}\simeq 2X_{n}$ at $T_{\gamma}\approx
0.07\,\mathrm{MeV}$.

The final helium abundances for different ADR parameters are shown and
compared to the expected standard value of
$Y_{{}^{4}\mathrm{He}}^{\mathrm{std}}\approx 0.227$ in table 3.

Table 3: Estimated helium abundances for different ADR parameters $b$ compared to the standard value $Y_{{}^{4}\mathrm{He}}^{\mathrm{SM}}=0.227\pm 0.004$. $b$ | $0$ | $10^{-17}$ | $10^{-15}$ | $10^{-12}$ | $10^{-6}$ | $10^{-4}$ | $10^{-2}$
---|---|---|---|---|---|---|---
$Y_{{}^{4}\mathrm{He}}$ | $0.235$ | $0.235$ | $0.235$ | $0.235$ | $0.227$ | $0.227$ | $0.227$

Here, we observe the same consistent picture as for our $\Delta
N_{\mathrm{eff}}$ observable. Small ADR parameters lead to a deviation $\Delta
Y_{{}^{4}\mathrm{He}}=\mathcal{O}(0.01)$ from the SM expectation much larger
than the experimental uncertainty $\sigma_{{}^{4}\mathrm{He}}\sim 0.004$ of
the observable $Y_{{}^{4}\mathrm{He}}$. On the other hand, very large ADR
parameters, i.e. $b\gtrsim 10^{-6}$ lead to discrepancies much smaller than
$\sigma_{{}^{4}\mathrm{He}}$ and hence would be in agreement with experiment.

The argument here is exactly the same as before since the depletion of
$\nu_{e}$ solely follows from the mixing behavior which in turn is dominantly
influenced by the resonance structure. In figure 4, we compare the temperature
evolution of $X_{n}$ for two different scenarios, i.e. for $b=10^{-17}$ and
for $b=10^{-4}$.

Figure 4: Evolution of the neutron abundance $X_{n}$ with respect to the
photon temperature $T_{\gamma}$ for two different ADR parameters, i.e.
$b\in\\{10^{-17},10^{-4}\\}$, after neutron freeze out. The dashed black curve
represents the SM expectation while the two solid curves correspond to the
benchmark scenarios, respectively.

Around the temperature of $\mathcal{O}(1\,\mathrm{MeV})$ $X_{n}$ departs from
equilibrium, as expected, in each scenario. However, $X_{n}(b=10^{-17})$
leaves equilibrium a little earlier and adopts a higher value compared to the
SM curve after neutron freeze out. At the temperature when all neutrons are
bound into helium nuclei, this leads to a higher helium abundance.

In the large ADR parameter scenario, i.e. $b=10^{-4}$, $X_{n}$ essentially
stays in agreement with the SM curve for the relevant temperatures. Therefore,
the corresponding helium abundance would also be in agreement with the SM
expectation within the experimental margin of error. This is due to the
negligible presence of $\nu_{s}$ and much higher interaction rates
$\lambda_{np},\lambda_{pn}$ compared to the previous case caused by the non-
dilution of the $\nu_{e}$ density.

###  4.2 $\Delta N_{\mathrm{eff}}$ and $Y_{{}^{4}\mathrm{He}}$ in the ALP only
Scenario

Next, we look at the behavior of ALP only models, where we turn off the ADR
potential and turn on the coupling of the ALP field to the sterile neutrino.
This results in a time dependent, additional mass for the sterile neutrino
mass matrix element, $m_{ss}\rightarrow m_{ss}+m_{s0}\eta(t)$, where $\eta$ is
given by eq. (12).

We integrate the QKEs for different parameter values shown in table 4. The
obtained results for $\Delta N_{\mathrm{eff}}$ show the clear pattern that
higher $m_{s0}$ values and lower $m_{\phi}$ values are favored by experimental
observation.

Table 4: Estimated additional light degrees of freedom $\Delta N_{\mathrm{eff}}$ at $x=50$ for different scalar field parameters. $m_{s0}$ / $\mathrm{eV}$ | $50$ | $100$ | $250$
---|---|---|---
$m_{\phi}$ / $\mathrm{eV}$ | $10^{-20}$ | $10^{-16}$ | $10^{-12}$ | $10^{-20}$ | $10^{-16}$ | $10^{-12}$ | $10^{-20}$ | $10^{-16}$ | $10^{-12}$
$\Delta N_{\mathrm{eff}}$ | $1.23$ | $1.22$ | $1.30$ | $0.80$ | $0.80$ | $0.93$ | $0.26$ | $0.26$ | $0.30$

We can explain this by looking at the behavior of the time dependent part
$\propto\eta$ of the sterile mass matrix element which is shown in figure 5
and the $\tan(2\theta)$ of the effective mixing angle shown in figure 6. The
first plot shows that if the scalar field is too heavy, it starts to oscillate
earlier leading to a decrease of the additional mass contribution of the
sterile neutrino. According to the latter figure, smaller sterile masses lead
to larger effective mixing angles which in turn lead to faster population of
the sterile species. Therefore, we need a large $m_{s0}$ and a small
$m_{\phi}$ to reconcile the existence of the sterile species with experimental
observations.

Figure 5: Behavior of the normalized scalar vev $\eta$ for different ALP mass
parameters $m_{\phi}$ at temperatures in the integration range of the QKEs.
Figure 6: Effective mixing angle between $\nu_{e}$ and $\nu_{s}$ for different
values of $m_{s0}$ at temperatures in the integration range of the QKEs. Here
we fixed $m_{\phi}=10^{-22}\,\mathrm{eV}$ so that $\eta$ remains constant for
all temperatures of interest.

These conjectures are supported by figures 7(a) and 7(b) showing the off-
diagonal element $|\varrho_{es}(y=5)|$. There we see that after the scalar
field has started to oscillate the correlations between active and sterile
neutrinos increase, whereas for an overall smaller $m_{s0}$ parameter the
correlations are also overall bigger.

(a) $m_{\phi}=10^{-20}\,\mathrm{eV}$.

(b) $m_{s0}=250\,\mathrm{eV}$.

Figure 7: As figure 3 for the ALP only scenario, for fixed $m_{\phi}$ (left)
and for fixed $m_{s0}$ (right). Shown are three different parameter
configurations per panel.

Now, we consider the obtaind helium abundances for the different benchmark
points. In table 5 the resulting values for $Y_{{}^{4}\mathrm{He}}$ are shown.
Here the same pattern arises as for $\Delta N_{\mathrm{eff}}$ with the
difference that our BBN observable is more sensitive to the ALP mass. For
masses $m_{\phi}\gtrsim 10^{-18}\,\mathrm{eV}$ the condensate oscillates
during or before nucleosynthesis has started resulting in a depleted electron
neutrino density. The earlier this oscillation occurs, the more $f_{\nu_{e}}$
is depleted and the bigger the deviation of the final helium abundance from
the SM value becomes as we can see from table 5. Only the benchmark points
with $m_{\phi}\gtrsim 10^{-16}\,\mathrm{eV}$ and $m_{s0}\gtrsim
100\,\mathrm{eV}$ are within the uncertainty around the standard value.

Table 5: Estimated helium abundance at $x=50$ for different scalar field parameters compared to the standard value $Y_{{}^{4}\mathrm{He}}^{\mathrm{SM}}=0.227\pm 0.004$. $m_{s0}$ / $\mathrm{eV}$ | $50$ | $100$ | $250$
---|---|---|---
$m_{\phi}$ / $\mathrm{eV}$ | $10^{-20}$ | $10^{-16}$ | $10^{-12}$ | $10^{-20}$ | $10^{-16}$ | $10^{-12}$ | $10^{-20}$ | $10^{-16}$ | $10^{-12}$
$Y_{{}^{4}\mathrm{He}}$ | $0.234$ | $0.234$ | $0.235$ | $0.232$ | $0.232$ | $0.233$ | $0.229$ | $0.229$ | $0.229$

To underline this statement we also show the evolution of the neutron fraction
for the benchmark points which are least and most compatible with observations
in figure 8. Note that the red curve deviates more significantly from the SM
expectation than the green one describing the most compatible parameter
configuration under consideration. The departure of the model curve for the
least compatible configuration due to the depletion of electron neutrinos is
more prominent than the almost non-existing one in the most compatible case.

Figure 8: As figure 4, for ALP only scenarios after neutron freeze out.

As we have seen, a bigger sterile mass matrix element due to the coupling of a
scalar field suppresses the mixing of $\nu_{e}$ and $\nu_{s}$ such that the
resulting helium fraction and effective degrees of freedom become compatible
with experimental bounds. Next, we consider the combined ADR and scalar field
scenario.

###  4.3 $\Delta N_{\mathrm{eff}}$ and $Y_{{}^{4}\mathrm{He}}$ in the combined
ADR and ALP Scenario

In the ADR only case, we have concluded that for sufficiently large ADR
parameters $b$ the equilibration of $\nu_{s}$ is suppressed. Choosing smaller
ADR parameters leads to a strong population of $\nu_{s}$ and hence large
corrections to $N_{\mathrm{eff}}$ that exceed experimental bounds. We expect
that in the combined scenario even small ADR parameters can be brought into
agreement with experiment by invoking the mixing suppression by the scalar
field $\phi$ from section 2.2. This expectation is further substantiated from
inspecting figure 9, which shows a (on average) smaller $|\varrho_{es}|$ than
for small $b$ values in the ADR only case.

Figure 9: Comparison of the evolution of $|\varrho_{es}(y=5)|$ for the ADR
only scenarios $b\in\\{10^{-17},10^{-6}\\}$ versus the combined scenario with
$m_{\phi}=10^{-22}\,\mathrm{eV}$, $m_{s0}=250\,\mathrm{eV}$ and $b=10^{-17}$.
The ADR only plots are shown as dashed lines, while the line corresponding to
the combined scenario is solid. Compare Figs. 3 and 7 for the individual ADR
and scalar field cases.

We confirm this expectation by choosing $m_{\phi}\sim 10^{-22}\,\mathrm{eV}$
and $m_{s0}\sim 250\,\mathrm{eV}$. For this value of $m_{\phi}$ the scalar
field starts oscillating long after nucleosynthesis has ceased and leads to a
constant addition to the sterile neutrino mass during the time of integration.
In table 6, we display the values of $\Delta N_{\mathrm{eff}}$ again for
$b\in\\{0,10^{-17},10^{-15},10^{-12}\\}$ with the addition of the scalar
field.

Table 6: Estimated additional light degrees of freedom $\Delta N_{\mathrm{eff}}$ at $x=50$ for different ADR parameters $b$ and $m_{\phi}\sim 10^{-22}\,\mathrm{eV}$ and $m_{s0}\sim 250\,\mathrm{eV}$. $b$ | $0$ | $10^{-17}$ | $10^{-15}$ | $10^{-12}$
---|---|---|---|---
$\Delta N_{\mathrm{eff}}$ | $0.25$ | $0.25$ | $0.25$ | $0.26$

We see that the values of $\Delta N_{\mathrm{eff}}$ decrease significantly
compared to the pure ADR scenario which is due to mixing suppression because
of the additional sterile mass as discussed in the last subsection. Now the
corrections to $N_{\mathrm{eff}}$ are in agreement with all bounds from
Planck, i.e. $\Delta N_{\mathrm{eff}}<\\{0.28,0.33,0.55,0.57\\}$. Thus, even
cases where the ADR scenario alone does not explain cosmological observation
an additional mass contribution for the sterile neutrino can reconcile them
with experiment. For comparison with the pure ADR case, we show the final
neutrino distributions for $b=10^{-17}$ and $b=10^{-12}$ in figure 10.

Figure 10: As figure 1 for the combined ADR and ALP scenario, for two
different ADR parameters. The shown $b$ values are excluded in the ADR only
scenario but are reconciled with experiment due to the addition of the scalar
field.

Of course this suppression effect disappears if one allows for higher scalar
masses such that the asymptotic sterile mass at $\mathcal{O}(1\,\mathrm{MeV})$
is reached before BBN or neutrino freeze out or if $m_{s0}$ is too small in
the first place. Consider for example the benchmark points

$\displaystyle p_{1}$
$\displaystyle:=(m_{\phi},m_{s0},b)=(10^{-22}\,\mathrm{eV},\mathbf{50}\,\mathrm{eV},10^{-17})\,,$
(37) $\displaystyle p_{2}$
$\displaystyle:=(m_{\phi},m_{s0},b)=(\mathbf{10^{-12}}\,\mathrm{eV},250\,\mathrm{eV},10^{-17})\,,$
(38)

where we adopted a smaller value for $m_{s0}$ for $p_{1}$ and a bigger value
for $m_{\phi}$ for $p_{2}$, respectively. Integrating the QKEs for the first
configuration yields

$\displaystyle\Delta N_{\mathrm{eff}}(p_{1})$ $\displaystyle\approx 1.23\,,$
(39)

which is slightly smaller than the original result $\Delta
N_{\mathrm{eff}}(b=10^{-17})\approx 1.36$ in the pure ADR case, but is still
much larger than for $m_{s0}=250\,\mathrm{eV}$. Hence, as expected, a smaller
effective sterile mass contribution leads to a reduced mixing suppression.
Increasing the scalar field mass as specified in eq. (38) for the second
benchmark point, we get

$\displaystyle\Delta N_{\mathrm{eff}}(p_{2})$ $\displaystyle\approx 0.30\,,$
(40)

which is slightly higher than the value we obtained for
$m_{\phi}=10^{-22}\,\mathrm{eV}$ but still in agreement with three out of four
Planck bounds. This is because after neutrino freeze out $\Delta
N_{\mathrm{eff}}$ remains constant and $\phi(t)$ only starts oscillating
shortly before. Thus, we obtain a similar result for $\Delta
N_{\mathrm{eff}}(p_{2})$ as for $m_{\phi}=10^{-22}\,\mathrm{eV}$. For
$m_{\phi}\gg 10^{-12}\,\mathrm{eV}$ the scalar field oscillates long enough
before neutrino freeze out to imply $\Delta
N_{\mathrm{eff}}(p_{2})\rightarrow\Delta N_{\mathrm{eff}}(b)$.

After having discussed the combined scenario for $\Delta N_{\mathrm{eff}}$, we
now turn towards our estimate of the helium abundance. Here, we expect the
same mechanism to apply as in section 4.2: assuming a scalar field mass of
$m_{\phi}\lesssim 10^{-18}\,\mathrm{eV}$ the oscillatory behavior of the
scalar field starts after $t_{\mathrm{BBN}}\sim 300\,\mathrm{s}$. Hence,
active-sterile oscillations are suppressed depending on the additional mass
$m_{s0}$ during the process of neutron freeze out.

We can observe this effect in figure 11 showing the evolution of $X_{n}(T)$
for $m_{\phi}=10^{-22}\,\mathrm{eV}$, $m_{s0}=250\,\mathrm{eV}$ and
$b=10^{-17}$ compared to the evolution in the pure ADR case with the same $b$
value.

Figure 11: As figure 4 comparing the pure ADR case with $b=10^{-17}$ with the
combined ADR and ALP scenario.

After adding the scalar field to the pure ADR model, one can no longer
distinguish the SM from the model curve within the given margin of error
$\sigma_{{}^{4}\mathrm{He}}=0.004$. This holds for all chosen ADR parameter
configurations with low $m_{\phi}$ and high $m_{s0}$, c.f table 7.

Table 7: Estimated helium abundances for different ADR parameters $b$ and $m_{\phi}\sim 10^{-22}\,\mathrm{eV}$ and $m_{s0}\sim 250\,\mathrm{eV}$ compared to the standard value $Y_{{}^{4}\mathrm{He}}^{\mathrm{SM}}=0.227\pm 0.004$. $b$ | $0$ | $10^{-17}$ | $10^{-15}$ | $10^{-12}$
---|---|---|---|---
$Y_{{}^{4}\mathrm{He}}$ | $0.229$ | $0.229$ | $0.229$ | $0.229$

Hence, in the combined scenario the deviation from the SM value is even
smaller than experimental uncertainties $\Delta Y_{{}^{4}\mathrm{He}}\sim
0.002<\sigma_{{}^{4}\mathrm{He}}$.

This effect gets weaker if we choose a lower $m_{s0}$ or increase the mass of
the scalar field up to values of $m_{\phi}\gg 10^{-18}\,\mathrm{eV}$. For
higher scalar masses, the ALP condensate already oscillates at times before
neutron-proton interactions freeze out and hence active-sterile oscillations
are not suppressed anymore regardless of the value of $m_{s0}$.

In order to demonstrate this effect, we again consider the benchmark points
$p_{1}$ and $p_{2}$ as in the $\Delta N_{\mathrm{eff}}$ analysis. For the
first point with lower $m_{s0}$, we again expect a less efficient mixing
suppression within the whole integration interval. This is indeed what we get
after integrating the QKEs for the neutrino density matrix and the neutron
fraction. The final helium abundance for $p_{1}$ amounts to

$\displaystyle Y_{{}^{4}\mathrm{He}}(p_{1})$ $\displaystyle\approx
0.234>0.229\,,$ (41)

which is larger than the value for $m_{s0}=250\,\mathrm{eV}$ and would be
observable in experiments.

For the second benchmark point with lower scalar mass, we again obtain a lower
value

$\displaystyle Y_{{}^{4}\mathrm{He}}(p_{2})$ $\displaystyle\approx 0.229\,.$
(42)

By adopting even smaller scalar mass parameters, we expect this value to
approach the pure ADR scenario. As a consequence, we would end up with more
helium-4.

This effect can be seen in figure 12.

Figure 12: Relative deviations of the $\nu_{e}$ and $\nu_{s}$ phase space
distribution $\Delta f_{\nu_{\alpha}}$ at $p_{2}$ and
$p_{\mathrm{ref}}=(10^{-22}\,\mathrm{eV},250\,\mathrm{eV},10^{-17})$ for
different fixed momentum values (black, red, green). The distributions are
plotted from the onset of scalar field oscillations $T_{\nu}(t\sim
m_{\phi}^{-1})$ down to the final temperature $T_{\nu,1}=0.02\,\mathrm{MeV}$,
with $m_{\phi}=10^{-12}\,\mathrm{eV}$.

As expected, the $\nu_{e}$ density decreases after $t\sim
m_{\phi}^{-1}=10^{14}\,\mathrm{eV}^{-1}$ while the sterile neutrino density
increases, compared to the benchmark point $p_{\mathrm{ref}}$ with
$m_{\phi}=10^{-22}\,\mathrm{eV}$ and all other parameters equal to those of
$p_{2}$. As a consequence, the higher $\nu_{s}$ density leads to an earlier
departure of $X_{n}$ from equilibrium, while the depleted $\nu_{e}$ density
causes an even bigger discrepancy between the SM curve and the corresponding
model curve due to the smaller $n-p$ reaction rates.

For higher $m_{\phi}$ this happens even earlier leading to a more significant
increase in $f_{\nu_{s}}$ (and a corresponding decrease in $f_{\nu_{e}}$).

##  5 Conclusions

In this paper we have analyzed the impact of altered dispersion relations
(ADRs) and couplings to an axion-like scalar field on cosmological bounds for
sterile neutrinos. Both effects have the potential to ameliorate such bounds,
depending on the concrete choice of parameters. In particular, we conclude
that ADR parameters in the range needed to give an explanation for short
baseline experiments, i.e $b=\mathcal{O}(10^{-17})$, alone are not sufficient
to suppress $\nu_{s}$ population in the early universe. We show this by
calculating the effective number of additional light degrees of freedom for
these parameters and by estimating the amount of helium produced during BBN.
This estimate results in the values

$\displaystyle\Delta N_{\mathrm{eff}}(b=10^{-17})$ $\displaystyle\approx
1.36\,,$ (43) $\displaystyle Y_{{}^{4}\mathrm{He}}(b=10^{-17})$
$\displaystyle\approx 0.235\,.$ (44)

Both quantities adopt values higher than allowed by experimental observations,
i.e. $\Delta
N_{\mathrm{eff}}\stackrel{{\scriptstyle\mathrm{Planck}}}{{\leq}}0.33,0.57$ and
$Y_{{}^{4}\mathrm{He}}=0.227_{\text{SM pred.}}\pm 0.004$.

In contrast, much larger ADR parameters, like $b\gtrsim 10^{-6}$, can indeed
suppress the $\nu_{s}$ population sufficiently to make light sterile neutrinos
compatible with early universe cosmology. For this range of parameters most
momentum modes experience mixing suppression because propagation and flavor
eigenstates almost coincide with each other resulting in the absence of active
sterile oscillations. Hence, we obtain

$\displaystyle\Delta N_{\mathrm{eff}}(b\gtrsim 10^{-6})$
$\displaystyle<0.04\,,$ (45) $\displaystyle Y_{{}^{4}\mathrm{He}}(b\gtrsim
10^{-6})$ $\displaystyle=0.227\,.$ (46)

Considering $b$-values differing from the short baseline anomaly scenario may
be, for example, motivated by effects making the ADR parameter dependent of
the cosmic evolution.

Moreover by adding the influence of an axion-like scalar field $\phi$ changing
the sterile neutrino mass $m_{ss}\to m_{ss}+m_{s0}\eta(m_{\phi},t)$ via a
Yukawa coupling, even very small ADR parameter values can be brought into
agreement with experiment. For a scalar field mass
$m_{\phi}=10^{-22}\,\mathrm{eV}$ and additional mass amplitude
$m_{s0}=250\,\mathrm{eV}$, we obtain

$\displaystyle\Delta
N_{\mathrm{eff}}(m_{\phi}=10^{-22}\,\mathrm{eV},m_{s0}=250\,\mathrm{eV},b=10^{-17})$
$\displaystyle\approx 0.26\,,$ (47) $\displaystyle
Y_{{}^{4}\mathrm{He}}(m_{\phi}=10^{-22}\,\mathrm{eV},m_{s0}=250\,\mathrm{eV},b=10^{-17})$
$\displaystyle\approx 0.229\,,$ (48)

which is compatible with observations.

While this analysis was carried out for a two neutrino generations framework,
i.e. one active and one sterile neutrino, our findings are expected to hold
even in scenarios involving greater numbers of generations. This expectation
is justified as it has been found that the active-sterile decoupling can be a
generic effect of the model [21] at energies higher than the resonance
energies. Also the additional sterile mass from the $\nu_{s}$-$\phi$ coupling
is expected to have the same effect if more neutrino generations are present.
Increasing the diagonal elements of the mass matrix corresponding to the
sterile species makes it dominantly diagonal resulting in suppressed active-
sterile mixing.

Furthermore, in this analysis we have neglected parameter configurations
leading to an early equilibrated species, i.e. at around
$T=\mathcal{O}(100\,\mathrm{MeV})$, for which the integration of the QKEs had
to be started much earlier at $T=\mathcal{O}(1\,\mathrm{GeV})$ as well as
finite temperature QED corrections. The former is justified since the
corresponding parameter configurations are not of interest to us since they
violate cosmological bounds by definition and would be excluded anyway.
Moreover, due to our findings described in section 3.2 and 4, we don’t expect
that our conclusions will be different if such an analysis is carried out.

Finite temperature QED corrections are important for precision predictions of
the total number of ultra relativistic degrees of freedom [52] since they can
lead to an increase in $\Delta N_{\mathrm{eff}}$ on the order of magnitude of
$\Delta N_{\mathrm{eff}}=\mathcal{O}(0.1)$. Here we were solely interested in
the impact of light sterile species on the number of additional neutrino
generations being mainly influenced by the oscillation Hamiltonian. Finite
temperature QED corrections only have a subleading influence on sterile
neutrinos because they do not interact directly with the electromagnetic
plasma.

In summary we find that the ADR-only scenario can only explain cosmological
observations if one assumes $b\gtrsim 10^{-6}$. The ALP only scenario works
well for $m_{s0}\gtrsim 100\,\mathrm{eV}$ and $m_{\phi}\lesssim
10^{-14}\,\mathrm{eV}$, whereas in the combined case also ADR parameters
compatible with SBL anomalies ($b\sim 10^{-17}$) can be brought into agreement
with experimental data for the same $(m_{\phi},m_{s0})$ configuration as in
the ALP only case.

Thus, if sterile neutrinos are discovered at future experiments, ADR effects
or a Yukawa coupling to a scalar condensate provide a promising explanation
why they did not reach thermal equilibrium in the early universe. By choosing
sufficiently high ADR parameters, high $m_{s0}$ and low $m_{\phi}$ these
effects lead to a suppression of $\nu_{s}$ population regardless of the
strength of vacuum mixing between active and sterile neutrinos.

##  Appendix A Evolution of the total Neutrino Density after Decoupling

To discuss the behavior of the total neutrino energy density after neutrino
decoupling, we first prove that the trace of the neutrino density matrix, i.e.
the sum of neutrino phase space distributions, is unaffected by neutrino
oscillations:

$\displaystyle(\partial_{t}-pH\partial_{p})\mathrm{Tr}(\varrho(t,p))$
$\displaystyle=-i\mathrm{Tr}([\mathcal{H}(t,p),\varrho(t,p)])+\mathrm{Tr}(\mathcal{C})[t,p,\varrho]$
(49)
$\displaystyle\Leftrightarrow(\partial_{t}-pH\partial_{p})\mathrm{Tr}(\varrho(t,p))$
$\displaystyle=\mathrm{Tr}(\mathcal{C})[t,p,\varrho]\,,$ (50)

where we made use of the fact, that the trace of a commutator always vanishes.
Thus, the sum of neutrino distribution functions is only changed by the
expansion of the universe and by neutrino scatterings with the background
plasma. After the freeze out of weak interactions, only the expansion of the
universe dilutes the total neutrino density and we get

$\displaystyle\partial_{t}\rho_{\nu}=-4H\rho_{\nu}=-3H(\rho_{\nu}+P_{\nu})\,,$
(51)

where $\rho_{\nu}$ is the sum of all neutrino energy densities. Thus, we
conclude that for a frozen out neutrino sector the total energy density obeys
the cosmic continuity equation, whereas before neutrino freeze out this
equation only holds for the total energy density of all radiation species.

##  Appendix B Choice of the Initial Scale Factor $x_{0}$

In section 3.2, we state our choice of the initial dimensionless scale factor
to be

$\displaystyle x_{0}=0.01\,,$ (52)

corresponding to an initial neutrino and photon temperature of

$\displaystyle T_{\nu}=T_{\gamma}=100\,\mathrm{MeV}\,.$ (53)

At this temperature, interactions between neutrinos and the eletromagnetic
plasma are sufficiently rapid such that $T_{\nu}=T_{\gamma}$ holds.
Furthermore, we know that the plasma at this temperature is mainly comprised
of the following particles

$\displaystyle
e^{-},e^{+},\\{\nu_{\alpha}\\}_{\alpha=e}^{s},\\{\bar{\nu}_{\alpha}\\}_{\alpha=e}^{s},\gamma\,,$
(54)

since the next heavier particles, i.e. muons and pions, are already non
relativistic and mainly decayed into these lighter particles. This
dramatically simplifies the collision terms we have to consider for the
neutrino Boltzmann equations. However, in order to verify the validity of the
chosen initial time, we have to integrate the Boltzmann equations at much
higher temperatures $T\gg m_{\mu}$ for a set of representative parameter
configurations.

In the following, we present an argument why it is still sufficient to
consider only the reactions listed in section 3.1 for these test runs. To show
that, we decompose the Boltzmann equation for the density matrix

$\displaystyle xH\frac{\partial\rho}{\partial x}(x,y)$
$\displaystyle=-i[\mathcal{H}(x,y),\rho(x,y)]+\mathcal{C}[\rho,x,y]\,,$ (55)

using the $SU(2)$ generator basis of $\mathbb{H}(2)$, i.e. the space of
Hermitian $2\times 2$ matrices. Thus, we expand all appearing matrices using
the extended set of scaled Pauli matrices

$\displaystyle\tau_{0}$ $\displaystyle=\frac{1}{2}\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\,,\quad\tau_{1}=\frac{1}{2}\begin{pmatrix}0&1\\\
1&0\end{pmatrix}\,,$ (56) $\displaystyle\tau_{2}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}0&-i\\\
i&0\end{pmatrix}\,,\quad\tau_{3}=\frac{1}{2}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\,.$ (57)

Before applying this decomposition, we rewrite the collision term as the
difference of gain and loss terms

$\displaystyle\mathcal{C}[\rho,x,y]$
$\displaystyle=\\{\Gamma^{+},\mathbb{I}-\rho\\}-\\{\Gamma^{-},\rho\\}\,,$ (58)

where $\Gamma^{+}\in\mathbb{H}(2)$ is the matrix valued collision rate for
gaining a neutrino from the plasma, while $\Gamma^{-}\in\mathbb{H}(2)$ is the
rate for losing a neutrino to the background (or to another momentum state).
Further using the linearity of the anticommutator yields

$\displaystyle\mathcal{C}[\rho,x,y]$
$\displaystyle=\\{\Gamma^{+},\mathbb{I}\\}-\\{\Gamma^{+},\rho\\}-\\{\Gamma^{-},\rho\\}$
(59) $\displaystyle=2\Gamma^{+}-\\{\Gamma^{+}+\Gamma^{-},\rho\\}\,.$ (60)

Now, we have to find the components of $\Gamma^{\pm}$ in our basis. To achieve
this, we note that sterile neutrinos do not interact with the plasma at all
which is why $\Gamma^{\pm}$ takes the form

$\displaystyle\Gamma^{\pm}=\underbrace{\gamma^{\pm}}_{\geq
0}\mathbb{P}_{a}\,,$ (61)

where $\mathbb{P}_{a}$ is the active neutrino projector. For one active and
one sterile neutrino flavor, it is given by

$\displaystyle\mathbb{P}_{a}=\begin{pmatrix}1&0\\\
0&0\end{pmatrix}=\tau_{0}+\tau_{3}\,.$ (62)

This automatically gives us the components of $\Gamma^{\pm}$

$\displaystyle\gamma_{0}^{\pm}=\gamma^{\pm}\,,\quad\gamma_{1}^{\pm}=0\,,\quad\gamma_{2}^{\pm}=0\,,\quad\gamma_{3}^{\pm}=\gamma^{\pm}\,.$
(63)

In the following, we employ the notation

$\displaystyle\mathbf{a}:=(a_{0},\vec{a}):=(a_{0},a_{1},a_{2},a_{3})\,,$ (64)

for the vector of components $(a_{k})_{k=0}^{3}$ of a Hermitian matrix $A$.
Moreover, the components of such a matrix are obtained by taking the scalar
product defined on $\mathbb{H}(2)$ of the matrix itself and the corresponding
basis matrix, i.e.

$\displaystyle a_{k}=\langle
A,\tau_{k}\rangle:=2\mathrm{Tr}(A\cdot\tau_{k})\,.$ (65)

The density matrix and Hamiltonian components are denoted as follows

$\displaystyle\rho$ $\displaystyle=\sum_{k=0}^{3}\varrho_{k}\tau_{k}\,,$ (66)
$\displaystyle\mathcal{H}$ $\displaystyle=\sum_{k=0}^{3}h_{k}\tau_{k}\,.$ (67)

Substituting this decomposition into eq. (55) yields the 4 coupled equations

$\displaystyle xH\frac{\partial\varrho_{0}}{\partial x}$
$\displaystyle=2\gamma^{+}-(\gamma^{+}+\gamma^{-})(\varrho_{0}+\varrho_{3})$
(68) $\displaystyle xH\frac{\partial\vec{\varrho}}{\partial x}$
$\displaystyle=\vec{h}\times\vec{\varrho}-(\gamma^{+}+\gamma^{-})\vec{\varrho}+(2\gamma^{+}-(\gamma^{+}+\gamma^{-})\varrho_{0})\vec{e}_{3}\,,$
(69)

with $\vec{e}_{3}=(0,0,1)^{T}$.

Now, we are able to argue that taking into account less reactions into our
collison terms yields an upper limit for the starting temperature888An upper
limit on $T_{\nu}$ corresponds to a lower limit for the starting scale factor
$x_{0}$.. This becomes appearant by considering the evolution equations for
the real and imaginary parts of the off-diagonal density matrix element
$\varrho_{1},\varrho_{2}$ as well as of the sterile neutrino component
$\rho_{ss}=(\varrho_{0}-\varrho_{3})/2$

$\displaystyle xH\frac{\partial\varrho_{1}}{\partial x}$
$\displaystyle=(\vec{h}\times\vec{\varrho})_{1}-(\gamma^{+}+\gamma^{-})\varrho_{1}\,,$
(70) $\displaystyle xH\frac{\partial\varrho_{2}}{\partial x}$
$\displaystyle=(\vec{h}\times\vec{\varrho})_{2}-(\gamma^{+}+\gamma^{-})\varrho_{2}\,,$
(71) $\displaystyle xH\frac{\partial\rho_{ss}}{\partial x}$
$\displaystyle=-(\vec{h}\times\vec{\varrho})_{3}\,.$ (72)

At $T\gg 100\,\mathrm{MeV}$, the total neutrino interaction rate
$\gamma^{+}+\gamma^{-}$ becomes sufficiently large leading to an exponential
dampening of the off-diagonal matrix elements. Consequently, the sterile
neutrino distribution remains constant at these times since
$\varrho_{1,2}\approx 0$ implies $\partial\rho_{ss}/\partial x\approx 0$.
Thus, if we start with an unpopulated sterile species at high temperatures it
will only start being populated as soon as the interaction rate of $\nu_{e}$
becomes weak enough such that oscillations aren’t suppressed anymore.

If we neglect some of the processes contributing to $\gamma^{\pm}$ the
temperature where $\gamma^{+}+\gamma^{-}$ is big enough to suppress
oscillations gets shifted to higher values, i.e. smaller $x$. Therefore, it is
safe to estimate a suitable starting temperature (scale factor) using an
incomplete set of interactions since the actual proper estimate will be lower
(higher). Of course, in case this procedure yields a starting temperature
$T_{\nu}>100\,\mathrm{MeV}$, we need to include the full set of reactions for
the actual integration process or deliver good arguments why $x_{0}=0.01$ is
still a good choice. In the following, we present the results of our test
runs.

Solving the set of Boltzmann equations in the high temperature regime is much
simpler because active neutrinos are in thermal equilibrium with the
electromagnetic plasma. Therefore, we use

$\displaystyle T_{\nu}=T_{\gamma}=\frac{\mathrm{MeV}}{x}\,,$ (73)
$\displaystyle f_{\nu_{e}}(x,y)=f_{\mathrm{eq}}(y)\,,$ (74)

to precalculate and interpolate the collision terms. Furthermore, we don’t
have to solve the evolution equation for the photon temperature.

In order to ensure that $\gamma^{+}+\gamma^{-}$ is large enough such that
oscillations are suppressed, we choose $T=10\,\mathrm{GeV}$ as our initial
temperature and integrate until $T=100\,\mathrm{MeV}$. To simplify the
procedure, we use the findings from above and only consider reactions with the
particles given in (54).
In figure 13, we show all components of the density matrix at $x_{0}=0.01$ for
a set of different models. Inspecting figure 13(a) shows that the electron
neutrino phasespace distribution perfectly aligns with its equilibrium value
as expected. Furthermore, we observe that for models with high ADR parameters
the sterile neutrino distribution is very close to 0 (cf. 13(d)) such that
numerical oscillations below machine precision take over. For much smaller ADR
parameters $\rho_{ss}$ can get very close to equilibrium especially for small
momenta already at $x_{0}=0.01$ indicating that $\rho_{ss}(x_{0})\equiv 0$
isn’t valid. Turning on the effective sterile mass contribution from the
coupling to the scalar field fixes this issue and $\rho_{ss}$ is close enough
to 0 for our initial assumption to hold. The same conclusions can be drawn
from 13(b) and 13(c) showing the off-diagonal density matrix elements which
quantify the correlations between $\nu_{e}$ and $\nu_{s}$.

(a) $\rho_{ee}$

(b) $\mathrm{Re}(\rho_{es})$

(c) $\mathrm{Im}(\rho_{es})$

(d) $\rho_{ss}$

Figure 13: Density matrix components at $x_{0}=0.01$ plotted against the
dimensionless momentum $y$ for different models displayed in the legends.

In figure 14, we compare the final sterile neutrino distribution from the test
run at $x_{0}=0.01$ to the distribution obtained in the full run at
$x=0.010001$ using the assumption that the sterile neutrino density vanishes
at $x_{0}$.

Figure 14: Comparison of the sterile neutrino densities obtained from the test
run at $x_{0}=0.01$ (red) and from the full run after 2 steps at
$x=x_{0}+10^{-6}$ (orange).

Here we can see that the sterile density already approaches and even surpasses
the values obtained in the test run at some momenta, hence we conclude that
the error introduced by using the wrong initial value is small.

In conclusion, we see that for sufficiently high ADR parameters or effective
sterile masses the initial condition

$\displaystyle\rho(x_{0}=0.01,y)\approx\begin{pmatrix}f_{\mathrm{eq}}(y)&0\\\
0&0\end{pmatrix}\,,$ (75)

is a good approximation.

Only for situations where $b$ is too small or mixing suppression through the
additional sterile mass is not sufficient, we obtain significant sterile
neutrino densities already at $x_{0}=0.01$. Thus, in order to achieve a
reliable estimate of the actual sterile neutrino density we need to start the
integration much earlier. However, since these models are already excluded for
the starting point $x_{0}$ (cf. 4) choosing an earlier $x_{0}$ would even lead
to higher or equal sterile neutrino densities. This is because oscillations
lead to an equilibration of the $\nu_{e}$ and $\nu_{s}$ number, while
collisions are pumping new electron neutrinos into the system as long as their
number deviates from the thermal equlibrium values. Therefore, we consider the
obtained sterile neutrino distributions as lower bounds of the actual values.
Moreover, it can be seen by comparing the density matrix components from the
test run at $x_{0}$ to the components obtained in the full run after a few
steps that they almost coincide. Hence, the estimate of $\rho(x,y)$ obtained
using $x_{0}=0.01$ for the problematic models is adequate for our purposes.

##  Appendix C Numerical Methodology

We solve this set of differential equations using C++ and the gnu scientific
library [53] (GSL). In detail, a fixed step size Burlisch-Stoer integrator is
employed and we take $\mathcal{O}(10^{4})$ steps between the initial
($x_{0}=0.01$) and final time ($x_{1}=50$).

Moreover, the calculation of all numerical integrals is carried out using the
closed Newton-Cotes formula of order 6, which is why we choose $N_{y}=109$. If
the density matrix is evaluated at momenta between the nodes defined by
$\tilde{\Omega}_{y}$, we interpolate it with a cubic Steffen Spline999The
Steffen Interpolation Method guaranties monotonicity of the cubic function
between two neighboring data points. [54].

##  Appendix D Thermal Interactionrates of Neutron-Proton Processes

The thermal reaction rates needed for the solution of the Boltzmann equation
for the neutron abundance are given by [55]

$\displaystyle\lambda_{np}(x)$ $\displaystyle=\lambda_{n+e^{+}\rightarrow
p+\bar{\nu}_{e}}+\lambda_{n+\nu_{e}\rightarrow p+e^{-}}+\lambda_{n\rightarrow
p+\bar{\nu}_{e}+e^{-}}\,,$ (76) $\displaystyle\lambda_{pn}(x)$
$\displaystyle=\lambda_{n+e^{+}\leftarrow
p+\bar{\nu}_{e}}+\lambda_{n+\nu_{e}\leftarrow p+e^{-}}+\lambda_{n\leftarrow
p+\bar{\nu}_{e}+e^{-}}\,,$ (77)

with the partial rates

$\displaystyle\lambda_{n+\nu_{e}\rightarrow p+e^{-}}(x)$
$\displaystyle=A\int\limits_{0}^{\infty}p_{\nu}^{2}p_{e}E_{e}[1-f_{e}(E_{e},T(x))]f_{\nu_{e}}(x,p_{\nu}a(x))\,\mathrm{d}p_{\nu}\,,$
(78) $\displaystyle\lambda_{n+e^{+}\rightarrow p+\bar{\nu}_{e}}(x)$
$\displaystyle=A\int\limits_{0}^{\infty}p_{e}^{2}p_{\nu}E_{\nu}[1-f_{\nu_{e}}(x,p_{\nu}a(x))]f_{e}(E_{e},T(x))\,\mathrm{d}p_{e}\,,$
(79) $\displaystyle\lambda_{n\rightarrow p+e^{-}+\bar{\nu}_{e}}(x)$
$\displaystyle=A\int\limits_{0}^{p_{0}}p_{e}^{2}p_{\nu}E_{\nu}(1-f_{\nu_{e}}(x,p_{\nu}a(x)))[1-f_{e}(E_{e},T(x))]\,\mathrm{d}p_{e}\,,$
(80)

where $p_{\nu}$ is the electron neutrino momentum, $p_{e}$ and $E_{e}$ are the
electron momentum and energy, respectively, and the common constant $A$ is
determined via the mean lifetime of the neutron in its rest frame $\tau_{n}$,
i.e.

$\displaystyle A\approx\left\\{0.0157\cdot\tau_{n}(\Delta
m)^{5}\right\\}^{-1}\,,\quad\mathrm{with}\quad\Delta m=m_{n}-m_{p}\,.$ (81)

The interaction rates for the back reactions are obtained from eqs. (78) to
(80) by replacing $f\rightarrow 1-f$ for all phase space distributrion
functions. Furthermore, the upper integration bound in eq. (80) is given by
$p_{0}=\sqrt{\Delta m^{2}-m_{e}^{2}}$ and the following energy conservation
relations are imposed:

  1. 1.

eq. (78): $E_{e}=E_{\nu}+\Delta m\,,$

  2. 2.

eq. (79): $E_{\nu}=E_{e}+\Delta m\,,$

  3. 3.

eq. (80): $E_{\nu}=\Delta m-E_{e}\,.$

## References

  * [1] A. A. Aguilar-Arevalo et al. [MiniBooNE], Phys. Rev. Lett. 121 (2018), 221801 doi:10.1103/PhysRevLett.121.221801 [arXiv:1805.12028 [hep-ex]].
  * [2] A. A. Aguilar-Arevalo et al. [LSND], Phys. Rev. D 64 (2001), 112007 doi:10.1103/PhysRevD.64.112007 [arXiv:0104049 [hep-ex]].
  * [3] G. Mention, M. Fechner, Th. Lasserre, Th. A. Mueller, D. Lhuillier, M. Cribier and A. Letourneau, Phys. Rev. D 83 (2011), 073006 doi:10.1103/PhysRevD.83.073006 [arXiv:1101.2755 [hep-ex]].
  * [4] C. Giunti and M. Laveder, Phys. Rev. C 83 (2011), 065504 doi:10.1103/PhysRevC.83.065504 [arXiv:1006.3244 [hep-ph]].
  * [5] P. Abratenko et al. [MicroBooNE], Phys. Rev. Lett. 128 (2022) no.24, 241801 doi:10.1103/PhysRevLett.128.241801 [arXiv:2110.14054 [hep-ex]].
  * [6] P. Abratenko et al. [MicroBooNE], Phys. Rev. D 105 (2022) no.11, 11 doi:10.1103/PhysRevD.105.112003 [arXiv:2110.14080 [hep-ex]].
  * [7] P. Abratenko et al. [(The MicroBooNE Collaboration)* and MicroBooNE], Phys. Rev. D 105 (2022) no.11, 112004 doi:10.1103/PhysRevD.105.112004 [arXiv:2110.14065 [hep-ex]].
  * [8] P. Abratenko et al. [MicroBooNE], Phys. Rev. D 105 (2022) no.11, 11 doi:10.1103/PhysRevD.105.112005 [arXiv:2110.13978 [hep-ex]].
  * [9] P. B. Denton, Phys. Rev. Lett. 129 (2022) no.6, 061801 doi:10.1103/PhysRevLett.129.061801 [arXiv:2111.05793 [hep-ph]].
  * [10] M. A. Acero, C. A. Argüelles, M. Hostert, D. Kalra, G. Karagiorgi, K. J. Kelly, B. Littlejohn, P. Machado, W. Pettus and M. Toups, et al. [arXiv:2203.07323 [hep-ex]].
  * [11] P. Adamson et al. [MINOS+], Phys. Rev. Lett. 122 (2019) no.9, 091803 doi:10.1103/PhysRevLett.122.091803 [arXiv:1710.06488 [hep-ex]].
  * [12] K. Abe et al. [T2K], Phys. Rev. D 99 (2019) no.7, 071103 doi:10.1103/PhysRevD.99.071103 [arXiv:1902.06529 [hep-ex]].
  * [13] P. Adamson et al. [NOvA], Phys. Rev. D 96 (2017) no.7, 072006 doi:10.1103/PhysRevD.96.072006 [arXiv:1706.04592 [hep-ex]].
  * [14] M. A. Acero et al. [NOvA], Phys. Rev. Lett. 127 (2021) no.20, 201801 doi:10.1103/PhysRevLett.127.201801 [arXiv:2106.04673 [hep-ex]].
  * [15] K. Abe et al. [Super-Kamiokande], Phys. Rev. D 91 (2015), 052019 doi:10.1103/PhysRevD.91.052019 [arXiv:1410.2008 [hep-ex]].
  * [16] M. G. Aartsen et al. [IceCube], Phys. Rev. Lett. 117 (2016), 071801 doi:10.1103/PhysRevLett.117.071801 [arXiv:1605.01990 [hep-ex]].
  * [17] H. Päs, S. Pakvasa and T. J. Weiler, Phys. Rev. D 72 (2005), 095017 doi:10.1103/PhysRevD.72.095017 [arXiv:hep-ph/0504096 [hep-ph]].
  * [18] D. Marfatia, H. Päs, S. Pakvasa and T. J. Weiler, Phys. Lett. B 707 (2012), 553-557 doi:10.1016/j.physletb.2012.01.028 [arXiv:1112.0527 [hep-ph]].
  * [19] E. Aeikens, H. Päs, S. Pakvasa and P. Sicking, JCAP 10 (2015), 005 doi:10.1088/1475-7516/2015/10/005 [arXiv:1410.0408 [hep-ph]].
  * [20] E. Aeikens, H. Päs, S. Pakvasa and T. J. Weiler, Phys. Rev. D 94 (2016) no.11, 113010 doi:10.1103/PhysRevD.94.113010 [arXiv:1606.06695 [hep-ph]].
  * [21] D. Döring, H. Päs, P. Sicking and Th. J. Weiler, Eur. Phys. J. C 80 (2020), 1202 doi:10.1140/epjc/s10052-020-08720-2 [arXiv:1808.07460 [hep-ph]].
  * [22] D. Jang, M. Kusakabe and M. K. Cheoun, AIP Conf. Proc. 1947 (2018) no.1, 030001 doi:10.1063/1.5030829
  * [23] G. A. Barenboim, P. Martínez-Miravé, C. A. Ternes and M. A. Tórtola, JHEP 03 (2020), 070 doi:10.1007/JHEP03(2020)070 [arXiv:1911.02329 [hep-ph]].
  * [24] D. Döring and H. Päs, Eur. Phys. J. C 79 (2019), 604 doi:10.1140/epjc/s10052-019-7122-8 [arXiv:1808.07734 [hep-ph]].
  * [25] N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [Erratum: Astron.Astrophys. 652 (2021), C4] doi:10.1051/0004-6361/201833910 [arXiv:1807.06209 [astro-ph]].
  * [26] A. Berlin, Phys. Rev. Lett. 117 (2016) no.23, 231801 doi:10.1103/PhysRevLett.117.231801 [arXiv:1608.01307 [hep-ph]].
  * [27] Y. Zhao, Phys. Rev. D 95 (2017) no.11, 115002 doi:10.1103/PhysRevD.95.115002 [arXiv:1701.02735 [hep-ph]].
  * [28] Y. Farzan, Phys. Lett. B 797 (2019), 134911 doi:10.1016/j.physletb.2019.134911 [arXiv:1907.04271 [hep-ph]].
  * [29] J. M. Cline, Phys. Lett. B 802 (2020), 135182 doi:10.1016/j.physletb.2019.135182 [arXiv:1908.02278 [hep-ph]].
  * [30] G. y. Huang and N. Nath, JCAP 05 (2022) no.05, 034 doi:10.1088/1475-7516/2022/05/034 [arXiv:2111.08732 [hep-ph]].
  * [31] G. y. Huang, M. Lindner, P. Martínez-Miravé and M. Sen, Phys. Rev. D 106 (2022) no.3, 033004 doi:10.1103/PhysRevD.106.033004 [arXiv:2205.08431 [hep-ph]].
  * [32] G. Krnjaic, P. A. N. Machado and L. Necib, Phys. Rev. D 97 (2018) no.7, 075017 doi:10.1103/PhysRevD.97.075017 [arXiv:1705.06740 [hep-ph]].
  * [33] V. Brdar, J. Kopp, J. Liu, P. Prass and X. P. Wang, Phys. Rev. D 97 (2018) no.4, 043001 doi:10.1103/PhysRevD.97.043001 [arXiv:1705.09455 [hep-ph]].
  * [34] J. Liao, D. Marfatia and K. Whisnant, JHEP 04 (2018), 136 doi:10.1007/JHEP04(2018)136 [arXiv:1803.01773 [hep-ph]].
  * [35] F. Capozzi, I. M. Shoemaker and L. Vecchi, JCAP 07 (2018), 004 doi:10.1088/1475-7516/2018/07/004 [arXiv:1804.05117 [hep-ph]].
  * [36] M. M. Reynoso, O. A. Sampayo and A. M. Carulli, Eur. Phys. J. C 82 (2022) no.3, 274 doi:10.1140/epjc/s10052-022-10228-w [arXiv:2203.11642 [hep-ph]].
  * [37] H. Davoudiasl and P. B. Denton, [arXiv:2301.09651 [hep-ph]].
  * [38] S. Hannestad, R. S. Hansen and T. Tram, Phys. Rev. Lett. 112 (2014) no.3, 031802 doi:10.1103/PhysRevLett.112.031802 [arXiv:1310.5926 [astro-ph.CO]].
  * [39] B. Dasgupta and J. Kopp, Phys. Rev. Lett. 112 (2014) no.3, 031803 doi:10.1103/PhysRevLett.112.031803 [arXiv:1310.6337 [hep-ph]].
  * [40] A. Mirizzi, G. Mangano, O. Pisanti and N. Saviano, Phys. Rev. D 91 (2015) no.2, 025019 doi:10.1103/PhysRevD.91.025019 [arXiv:1410.1385 [hep-ph]].
  * [41] Y. Tang, Phys. Lett. B 750 (2015), 201-208 doi:10.1016/j.physletb.2015.09.018 [arXiv:1501.00059 [hep-ph]].
  * [42] F. Forastieri, M. Lattanzi, G. Mangano, A. Mirizzi, P. Natoli and N. Saviano, JCAP 07 (2017), 038 doi:10.1088/1475-7516/2017/07/038 [arXiv:1704.00626 [astro-ph.CO]].
  * [43] M. Archidiacono, S. Gariazzo, C. Giunti, S. Hannestad and T. Tram, JCAP 12 (2020), 029 doi:10.1088/1475-7516/2020/12/029 [arXiv:2006.12885 [astro-ph.CO]].
  * [44] G. Sigl and G. Raffelt, Nucl. Phys. B 406 (1993), 423-451 doi:10.1016/0550-3213(93)90175-O
  * [45] A. Vlasenko, G. M. Fuller and V. Cirigliano Phys. Rev. D 89 (2014), 105004 doi:10.1103/PhysRevD.89.105004 [arXiv:1309.2628 [hep-ph]].
  * [46] L. Stodolsky, Phys. Rev. D 36 (1987), 2273 doi:10.1103/PhysRevD.36.2273
  * [47] P. F. de Salas and S. Pastor, JCAP 07 (2016), 051 doi:10.1088/1475-7516/2016/07/051 [arXiv:1606.06986 [hep-ph]].
  * [48] D. N. Blaschke and V. Cirigliano, Phys. Rev. D 94 (2016), 033009 doi:10.1103/PhysRevD.94.033009 [arXiv:1605.09383 [hep-ph]].
  * [49] S. Gariazzo, P. F. de Salas and S. Pastor, JCAP 07 (2019), 014 doi:10.1088/1475-7516/2019/07/014 [arXiv:1905.11290 [astro-ph.CO]].
  * [50] J. Bernstein, Cambridge University Press (1988), doi:10.1017/CBO9780511564185
  * [51] E. Aver, K. A. Olive and E. D. Skillman, JCAP 07 (2015), 011 doi:10.1088/1475-7516/2015/07/011 [arXiv:1503.08146 [astro-ph]].
  * [52] J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor and Y. Y. Y. Wong, JCAP 04 (2021), 073 doi:10.1088/1475-7516/2021/04/073 [arXiv:2012.02726 [hep-ph]].
  * [53] M. Galassi et al., Network Theory Ltd. (2009), ISBN:978-0954612078
  * [54] M. Steffen, Astronomy and Astrophysics 239 (1990), 443
  * [55] S. Weinberg, John Wiley and Sons (1972), ISBN:978-0-471-92567-5