Light sterile neutrinos in the early universe: effects of altered dispersion relations and a coupling to axion-like dark matter

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.

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 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 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 eff := 3+∆N 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), Y4 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 ϕ 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 [37][38][39][40][41][42][43] 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 ν s -ϕ coupling is small and consequently, the induced matter potential and non-linear self interactions of ν 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 ϕ.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 ∆N eff and the helium fraction Y4 He , we solve the quantum kinetic equations (QKEs) for the neutrino density matrix ϱ.This defines the main task of this work, since from ϱ 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 ∆N 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

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 dimensions 1 leads to new terms in the propagation Hamiltonian.Employing the usual ultra-relativistic expansion of the neutrino dispersion relation yields where p is the average neutrino momentum, b is the so called ADR parameter controlling the strength of the ADR effect, P s is the sterile neutrino projector and M is the neutrino mass matrix.
In the flavor basis, the last two quantities read for the 2 × 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 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] 1 The 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.
with the Fermi coupling G F , the electron and neutrino energy densities ρ α 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 and find the following relation for the mixing angle The mixing becomes maximal as soon as θ = π/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 The momentum p 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.

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] Here, J 1 4 is the regular Bessel J function of fractional order 1/4, m ϕ is the mass of the scalar and t is the cosmic time.Moreover, we allow for a feeble 2 coupling to sterile neutrinos via a Yukawa interaction term in its Lagrangian, Since we assume the coupling constant λ to be very small, we neglect possible interactions between ν s and ϕ quanta and interactions of ν s with itself mediated by ϕ.Nevertheless, we account for an additional time dependent sterile neutrino mass term modifying the mixing between electron and sterile neutrinos.Although this additional mass term is proportional to λ, it is indeed significant if the amplitude ϕ 0 of the classical component of ϕ is sufficiently large.Assuming ϕ to contribute a substantial amount to the dark matter density, one can deduce [29] ϕ 0 ∼ 10 15 GeV 10 −15 eV m ϕ 1 4 .
Thus, for scalar masses on the order of m ϕ ∼ 10 −10 eV one can still obtain λ • ϕ 0 ∼ 100 eV for λ ∼ 10 −21 .At early times, the mass contribution remains approximately constant leading to a constant sterile neutrino mass m eff ≈ m ss + λϕ 0 .If ϕ 0 is sufficiently large, the effective mixing of electron and sterile neutrinos is negligible since tan(2θ) is suppressed by the large mass gap m 2 eff − m 2 ee in the denominator of eq. ( 9) with m s s replaced by m 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 ∼ 1/m ϕ ) the effective sterile neutrino mass approaches m ss , allowing for significant oscillations again.Therefore, the neutrino oscillation behavior depends on the mass parameter m ϕ , i.e. the smaller m ϕ the longer ϕ(t) will remain approximately constant and active-sterile mixing will be suppressed by m s0 .

Central Quantities and Parameter Space
In the following analysis, we focus on two major quantities: 1.The effective number of additional neutrino generations ∆N eff (t) : = 8 7 11 4 Here, ρ k (t) are the neutrino energy densities at time t and ρ γ (t) is the corresponding photon density 3 .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.

The helium mass fraction
with the helium and baryon number densities n He , n 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 ∆N 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 Y4 He has to be computed shortly after BBN.It is, however, sufficient to evaluate ∆N eff directly after e ± annihilation since from there on it remains constant.This is because after e ± annihilation the total neutrino energy density only changes due to the expansion of the universe 4 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 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.The scalar field mass m ϕ determining when the scalar field starts to oscillate.
2. The amplitude m s0 := λϕ 0 of the additional mass contribution for the sterile neutrino.

The ADR parameter b.
In the following, we assume a range for m s0 ∈ [10 eV, 250 eV] in accordance with ref. [29].Furthermore, for m ϕ we choose the interval of possible values to be [10 −22 eV, 10 −14 eV] since for m ϕ ≤ 10 −22 eV the scalar field starts oscillating so late that the effective sterile neutrino mass remains constant during the considered temperatures.Thus, for m ϕ ≲ 10 −22 eV our scenario becomes independent of m ϕ and yields the same constraints on m s0 as for m ϕ ∼ 10 −22 eV.
At values larger than m ϕ = 10 −14 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, ∞), i.e. we consider anything between zero ADR potential and an arbitrarily large ADR effect.A special point in this parameter range is b ∼ 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.

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 (⃗ p) of neutrino mass eigenstates, i.e. [44,45] (2π where the thermal average of an operator Ô is defined using the density operator5 Π of the thermal system as usual ⟨ Ô⟩ := Tr( Π Ô).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 ν 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 ν k momentum states on average are occupied in the system.Moreover, by inspecting the off-diagonal elements of ϱ, we obtain information about the average correlation between ν j (p) and ν 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 ϱ jk into this basis using the neutrino mixing matrix U from eq. ( 8) via From now on, we will work with ϱ f only and drop the superscript f .

Quantum Kinetic Equations for the Density Matrix
The time evolution of ϱ in an expanding, homogeneous and isotropic universe is governed by a Boltzmann-like, quantum kinetic equation (QKE) [44,45,[47][48][49]] where p is the modulus of the neutrino momentum, t is the cosmic time, H is the Hubble rate, H is the Hamiltonian from eq. ( 6) and C is the collison operator.The convectional derivative operator, ∂ t − pH∂ p , on the left hand side includes the effect of the expansion of the universe redshifting the neutrino momentum as p ∝ 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 C. While the commutator part governs the evolution of ρ 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.
Table 1: All relevant processes considered in the neutrino collision terms.
k process The collision operator is the sum of the individual collision operators corresponding to a scattering process involving neutrinos where the set of all processes P contains the interactions given in table 1.Here we neglect neutrinonucleon scattering processes due to strong Boltzmann suppression of the nucleon distribution functions at temperatures of O(100 MeV) and smaller.For example, the collision term for the process νe − ↔ νe − is given by where d 3 ⃗ π j := (2E j (2π) 3 ) −1 d 3 ⃗ p j denotes the Lorentz invariant phase space measure, g L = 1/2 + sin 2 (θ W ), g R = sin 2 (θ W ), θ W is the Weinberg angle and the statistical factor ϕ νe − νe − is given by with )]P e , ϱ(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 γ of the electromagnetic plasma, we need to solve the continuity equation of the universe where ρ 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 equilibrium 6 , this equation can be reformulated into a differential equation for T γ .

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 where we choose m 0 = 1 MeV.Therefore, x represents the dimensionless scale factor and y is a dimensionless momentum variable not being redshifted over time, since p ∝ a −1 .Moreover, x takes the role of the reciprocal of the neutrino temperature which is equal to T γ at early times but deviates from it after neutrino decoupling and electron positron annihilation.Transformed to these new variables eq. ( 20) assumes the form with ρ 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. ρ → ϱ.
In order to integrate eq. ( 27), in principle we had to start at x 0 = 0 and set ϱ ik (x 0 , y) ≡ 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.Active neutrinos are in thermal equilibrium with the electromagnetic plasma.
2. Quarks and gluons are bound into hadrons.
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 ϱ is given by We found that x 0 = 0.01, i.e.T γ,0 = 100 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: 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 Ω y = [0, ∞) 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 min , y 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 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 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 The discretized version of Ω y then reads 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 solved 7 .

Calculating the helium Abundance
In this section, we present how the 4 He mass fraction Y4 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.
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 7 For the numerical details see appendix C because the baryon number in a comoving volume, N B ≈ a 3 (n n + n p ), is conserved within all relevant processes shortly before neutron freeze out.These processes are n + e + ↔ p + νe (32) n + ν e ↔ p + e − (33) The differential equation governing the evolution of X n reads [50] dX in terms of x = m 0 a(t).Here λ np is the thermal interaction rate of all processes converting neutrons to protons, while λ 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 MeV) where neutrons and protons are still in thermal equilibrium up until x(T = 0.07 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.
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.Here the abbreviations TT, TE, EE, lowE, lensing, BAO refer to different measurement techniques / features of the CMB data (i.e.TT = intensity (temperature) only, TE = temperature + curl free polarization data, EE = curl free polarization data only, lowE = curl free polarization data only at low multipole moments, lensing = grav.lensing measurement, BAO = 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. σ4 He = 0.004.

∆N eff and Y4 He in the pure ADR Scenario
The values of ∆N eff obtained after the full integration of QKEs for different ADR parameters are shown in table 2.
Here, we see that the resulting ∆N 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 ∆N 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 ≳ 10 −6 ) leads to a decrease in ∆N 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.
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 res for multiple ADR parameters in the temperature range T ν ∈ T ν := [3, 100] MeV.For larger ADR parameters the resonance curve passes through the relevant y region from above leading to resonantly enhanced  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 ∆N 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 |ϱ es | := Re(ϱ es ) 2 + Im(ϱ es ) 2 for y = 5.
The off-diagonal element ϱ es is important since it contains information about the energy transfer from ν e to ν s .Therefore, the plot shows that shortly after T = O(100 MeV) for models with vanishing or small ADR parameters a dip occurs in ϱ es leading to a significant enhancement of ρ ss , afterwards.Considering, the curve for b = 10 −15 , we see that the resonance around O(7 MeV) leads to a significant impact in |ϱ 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 ≥ 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, ∆N 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 γ = O(3 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 ν 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 ν 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 Y4 He ≃ 2X n at T γ ≈ 0.07 MeV.
The final helium abundances for different ADR parameters are shown and compared to the expected standard value of Y std 4 He ≈ 0.227 in table 3. Here, we observe the same consistent picture as for our ∆N eff observable.Small ADR parameters lead to a deviation ∆Y4 He = O(0.01)from the SM expectation much larger than the experimental uncertainty σ4 He ∼ 0.004 of the observable Y4 He .On the other hand, very large ADR parameters, i.e. b ≳ 10 −6 lead to discrepancies much smaller than σ4 He and hence would be in agreement with experiment.
The argument here is exactly the same as before since the depletion of ν 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 .Around the temperature of O(1 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 ν s and much higher interaction rates λ np , λ pn compared to the previous case caused by the non-dilution of the ν e density.

∆N eff and Y4 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 → m ss + m s0 η(t), where η is given by eq. ( 12).
We  5 and the tan(2θ) 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 ϕ to reconcile the existence of the sterile species with experimental observations.These conjectures are supported by figures 7a and 7b showing the off-diagonal element |ϱ 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.Now, we consider the obtaind helium abundances for the different benchmark points.In table 5 the resulting values for Y4 He are shown.Here the same pattern arises as for ∆N eff with the difference that our BBN observable is more sensitive to the ALP mass.For masses m ϕ ≳ 10 −18 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 ν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 ϕ ≳ 10 −16 eV and m s0 ≳ 100 eV are within the uncertainty  around the standard value.
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.As we have seen, a bigger sterile mass matrix element due to the coupling of a scalar field suppresses the mixing of ν e and ν s such that the resulting helium fraction and effective degrees of freedom become compatible with experimental bounds.Next, we consider the combined ADR and

∆N eff and Y4 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 ν s is suppressed.Choosing smaller ADR parameters leads to a strong population of ν s and hence large corrections to N 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 ϕ from section 2.2.This expectation is further substantiated from inspecting figure 9, which shows a (on average) smaller |ϱ es | than for small b values in the ADR only case.We confirm this expectation by choosing m ϕ ∼ 10 −22 eV and m s0 ∼ 250 eV.For this value of m ϕ 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 ∆N eff again for b ∈ {0, 10 −17 , 10 −15 , 10 −12 } with the addition of the scalar field.
We see that the values of ∆N 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 eff are in agreement with all bounds from Planck, i.e.Of course this suppression effect disappears if one allows for higher scalar masses such that the asymptotic sterile mass at O(1 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 where we adopted a smaller value for m s0 for p 1 and a bigger value for m ϕ for p 2 , respectively.Integrating the QKEs for the first configuration yields which is slightly smaller than the original result ∆N eff (b = 10 −17 ) ≈ 1.36 in the pure ADR case, but is still much larger than for m s0 = 250 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 which is slightly higher than the value we obtained for m ϕ = 10 −22 eV but still in agreement with three out of four Planck bounds.This is because after neutrino freeze out ∆N eff remains constant and ϕ(t) only starts oscillating shortly before.Thus, we obtain a similar result for ∆N eff (p 2 ) as for After having discussed the combined scenario for ∆N 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 ϕ ≲ 10 −18 eV the oscillatory behavior of the scalar field starts after t BBN ∼ 300 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 ϕ = 10 −22 eV, m s0 = 250 eV and b = 10 −17 compared to the evolution in the pure ADR case with the same b value.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 σ4 He = 0.004.This holds for all chosen ADR parameter configurations with low m ϕ and high m s0 , c.f table 7. Hence, in the combined scenario the deviation from the SM value is even smaller than experimental uncertainties ∆Y4 He ∼ 0.002 < σ4 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 ϕ ≫ 10 −18 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 ∆N 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 which is larger than the value for m s0 = 250 eV and would be observable in experiments.
For the second benchmark point with lower scalar mass, we again obtain a lower value 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.As expected, the ν e density decreases after t ∼ m −1 ϕ = eV and all other parameters equal to those of p 2 .As a consequence, the higher ν s density leads to an earlier departure of X n from equilibrium, while the depleted ν 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 ϕ this happens even earlier leading to a more significant increase in f νs (and a corresponding decrease in f νe ).

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 = O(10 −17 ), alone are not sufficient to suppress ν 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 Y4 He (b = 10 −17 ) ≈ 0.235 .
In contrast, much larger ADR parameters, like b ≳ 10 −6 , can indeed suppress the ν 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 Y4 He (b ≳ 10 −6 ) = 0.227 .
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 ϕ changing the sterile neutrino mass m ss → m ss + m s0 η(m ϕ , t) via a Yukawa coupling, even very small ADR parameter values can be brought into agreement with experiment.For a scalar field mass m ϕ = 10 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 ν s -ϕ 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 = O(100 MeV), for which the integration of the QKEs had to be started much earlier at T = O(1 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 ∆N eff on the order of magnitude of ∆N eff = 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 ≳ 10 −6 .The ALP only scenario works well for m s0 ≳ 100 eV and m ϕ ≲ 10 −14 eV, whereas in the combined case also ADR parameters compatible with SBL anomalies (b ∼ 10 −17 ) can be brought into agreement with experimental data for the same (m ϕ , 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 ϕ these effects lead to a suppression of ν s population regardless of the strength of vacuum mixing between active and sterile neutrinos.

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: 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 where ρ ν 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.
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 corresponding to an initial neutrino and photon temperature of At this temperature, interactions between neutrinos and the eletromagnetic plasma are sufficiently rapid such that T ν = T γ holds.Furthermore, we know that the plasma at this temperature is mainly comprised of the following particles e − , e + , {ν α } s α=e , {ν α } s α=e , γ , 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 ≫ m µ 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 using the SU (2) generator basis of H(2), i.e. the space of Hermitian 2 × 2 matrices.Thus, we expand all appearing matrices using the extended set of scaled Pauli matrices Before applying this decomposition, we rewrite the collision term as the difference of gain and loss terms where Γ + ∈ H(2) is the matrix valued collision rate for gaining a neutrino from the plasma, while Γ − ∈ 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 = 2Γ + − {Γ + + Γ − , ρ} .
Now, we have to find the components of Γ ± in our basis.To achieve this, we note that sterile neutrinos do not interact with the plasma at all which is why Γ ± takes the form where P a is the active neutrino projector.For one active and one sterile neutrino flavor, it is given by This automatically gives us the components of Γ ± γ ± 0 = γ ± , γ ± 1 = 0 , γ ± 2 = 0 , γ ± 3 = γ ± .(63) In the following, we employ the notation a := (a 0 , ⃗ a) := (a 0 , a 1 , a 2 , a 3 ) , for the vector of components (a k ) 3 k=0 of a Hermitian matrix A. Moreover, the components of such a matrix are obtained by taking the scalar product defined on H(2) of the matrix itself and the corresponding basis matrix, i.e.
The density matrix and Hamiltonian components are denoted as follows Substituting this decomposition into eq.( 55) yields the 4 coupled equations xH ∂ϱ 0 ∂x = 2γ + − (γ + + γ − )(ϱ 0 + ϱ 3 ) (68) with ⃗ 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 temperature 8 .This becomes appearant by considering the evolution equations for the real and imaginary parts of the off-diagonal density matrix element ϱ 1 , ϱ 2 as well as of the sterile neutrino component ρ ss = (ϱ 0 − ϱ 3 )/2 At T ≫ 100 MeV, the total neutrino interaction rate γ + + γ − 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 ϱ 1,2 ≈ 0 implies ∂ρ ss /∂x ≈ 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 ν e becomes weak enough such that oscillations aren't suppressed anymore.
If we neglect some of the processes contributing to γ ± the temperature where γ + + γ − 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 ν > 100 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 f νe (x, y) = f eq (y) , 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 γ + + γ − is large enough such that oscillations are suppressed, we choose T = 10 GeV as our initial temperature and integrate until T = 100 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 13a shows that the electron neutrino phasespace distribution perfectly

1 .
Neutrino interactions are completely frozen out 2. All free neutrons are bound into light nuclei 3. The neutrino distribution functions have reached their asymptotic values 4. The relativistic approximation for the oscillation Hamiltonian and collisions is valid

1 Figure 5 :
Figure 5: Behavior of the normalized scalar vev η for different ALP mass parameters m ϕ at temperatures in the integration range of the QKEs.

𝑚 1 Figure 6 :
Figure 6: Effective mixing angle between ν e and ν s for different values of m s0 at temperatures in the integration range of the QKEs.Here we fixed m ϕ = 10 −22 eV so that η remains constant for all temperatures of interest.

Figure 7 :𝑋 1 Figure 8 :
Figure 7: As figure 3 for the ALP only scenario, for fixed m ϕ (left) and for fixed m s0 (right).Shown are three different parameter configurations per panel.

5 ( 1 Figure 9 :Table 6 :
Figure 9: Comparison of the evolution of |ϱ es (y = 5)| for the ADR only scenarios b ∈ {10 −17 , 10 −6 } versus the combined scenario with m ϕ = 10 −22 eV, m s0 = 250 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.

𝜈 (𝑦) 1 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.

6 𝛥𝑓𝛥𝑓 1 Figure 12 :
Figure 12: Relative deviations of the ν e and ν s phase space distribution ∆f να at p 2 and p ref = (10 −22 eV, 250 eV, 10 −17 ) for different fixed momentum values (black, red, green).The distributions are plotted from the onset of scalar field oscillations T ν (t ∼ m −1 ϕ ) down to the final temperature T ν,1 = 0.02 MeV, with m ϕ = 10 −12 eV. 10 14 eV −1 while the sterile neutrino density increases, compared to the benchmark point p ref with m ϕ = 10 −22eV and all other parameters equal to those of p 2 .As a consequence, the higher ν s density leads to an earlier departure of X n from equilibrium, while the depleted ν 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 ϕ this happens even earlier leading to a more significant increase in f νs (and a corresponding decrease in f νe ).

Table 2 :
Estimated additional light degrees of freedom ∆N eff at x = 50 for different ADR parameters b. b 0 10 −17 10 −15 10 −12 10 −6 10 −4 10 −2Now, 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 ∆N eff and the estimated helium-4 abundance Y4 He .These simulated values for ∆N eff are compared to bounds obtained by the Planck collaboration, i.e.

Table 3 :
Estimated helium abundances for different ADR parameters b compared to the standard value Y SM 4 He = 0.227 ± 0.004.

Table 4 :
4ntegrate the QKEs for different parameter values shown in table4.The obtained results for ∆N eff show the clear pattern that higher m s0 values and lower m ϕ values are favored by experimental observation.We can explain this by looking at the behavior of the time dependent part ∝ η of the sterile Figure 4: Evolution of the neutron abundance X n with respect to the photon temperature T γ for two different ADR parameters, i.e. b ∈ {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.Estimated additional light degrees of freedom ∆N eff at x = 50 for different scalar field parameters./ eV 10 −20 10 −16 10 −12 10 −20 10 −16 10 −12 10 −20 10 −16 10 −12