Inflation, (P)reheating and Neutrino Anomalies: Production of Sterile Neutrinos with Secret Interactions

A number of experimental anomalies involving neutrinos hint towards the existence of at least an extra (a very light) sterile neutrino. However, such a species, appreciably mixing with the active neutrinos, is disfavored by different cosmological observations like Big Bang Nucleosynthesis (BBN), Cosmic Microwave Background (CMB) and Large Scale Structure (LSS). Recently, it was shown that the presence of additional interactions in the sterile neutrino sector via light bosonic mediators can make the scenario cosmologically viable by suppressing the production of the sterile neutrinos from active neutrinos via matter-like effect caused by the mediator. This mechanism works assuming the initial population of this sterile sector to be negligible with respect to that of the Standard Model (SM) particles, before the production from active neutrinos. However, there is fair chance that such bosonic mediators may couple to the inflaton and can be copiously produced during (p)reheating epoch. Consequently, they may ruin this assumption of initial small density of the sterile sector. In this article we, starting from inflation, investigate the production of such a sterile sector during (p)reheating in a large field inflationary scenario and identify the parameter region that allows for a viable early Universe cosmology.


Introduction
Several anomalies from different experiments measuring neutrino oscillations have hinted towards the existence of an additional sterile neutrino species. While LSND [1,2] and MiniBooNE [3] reported an excess inν µ →ν e and the latter have also indicated an excess of ν e in the ν µ beam. Within a 3+1 framework, MiniBooNE result hints towards the existence of a sterile neutrino with eV mass at 4.8σ significance, which raises to 6.1σ when combined with the LSND data. Further, Daya Bay [4], NEOS [5], DANSS [6] and other reactor experiments [7][8][9] probed the ν e disappearance in thē ν e →ν e channel, whereas Gallium experiments [10][11][12] like GALLEX [13], SAGE [14] have performed similar measurements in the ν e → ν e channel. Theν e disappearance data also hints in favour of sterile neutrinos at 3σ level.
However, there are significant tension among different neutrino experiments. In particular, observed excess in the experiments measuring ν µ (ν µ ) → ν e (ν e ) appearance (i.e. LSND and MiniBooNE) are in tension with strong constraints on ν µ disappearance, mostly from MINOS [15] and IceCUBE [16], while attempting to fit together using a 3+1 framework [17]. Thus, the existence of a light (O(1) eV) sterile neutrino within a simple 3+1 framework, as a possible resolution to the ν e appearance anomalies, remains debatable. 1 However, such a light additional sterile neutrino, with mixing sin θ O(0.1) with the active neutrino species, can be consistent with constraints from various terrestrial neutrino experiments.
In the early Universe, production of a light sterile neutrino, if exists, can be significant. Thanks to its sizable mixing with the active neutrinos, it can be produced via Dodelson-Widrow (DW) mechanism [19,20]. Further, inflaton decays during reheating, or any other heavy scalar particle can possibly decay into sterile neutrinos [21][22][23]. Due to its sizable mixing with the active neutrinos, thermalization with the SM particles are also ensured. However, several cosmological constraints disfavor the viability of such a scenario. In particular, constraints from Big Bang Nucleosynthesis (BBN) [24][25][26][27] restricts the effective number of relativistic degrees of freedom. This is because it would enhance the expansion rate at the onset of BBN ( 1 MeV). Depending on its mass, such a species can contribute as either matter or radiation during matter-radiation equality epoch (for more explanation see [28]). Additional non-relativistic neutrinos can also affect late-time expansion rate. Thus, it can affect the position of the acoustic peaks. Further, such light species can lead to slow down of Dark Matter clustering and thanks to their large free-streaming length [29], can wash-out small-scale structure. Thus, Cosmic Microwave Background (CMB), together with Baryon Acoustic Oscillation (BAO) [30][31][32] and Ly-α measurements [33,34] put forward significant constraint on the total neutrino mass Σm ν as well as on the number of relativistic degrees of freedom N eff [35]. Both of these constraints impact the viability of an additional light sterile neutrino species.
Planck [36], assuming three neutrinos with degenerate mass, Fermi-Dirac distribution and zero chemical potential, constrains the properties of neutrinos. If they be-come non-relativistic after recombination, they mainly affect CMB through the change of angular diameter distance, which is degenerate with H 0 . So, the cleanest signal is through lensing power spectra that in turn affect the CMB power spectra. Since neutrino mass suppresses lensing whereas CMB prefers higher lensing, neutrino mass is strongly constrained by CMB lensing data. Neutrinos with large mass, that become non-relativistic around recombination, can produce distinctive features in CMB (such as reducing the first peak height) and are thus ruled out. Planck constrains Σm ν < 0.12 eV and N eff = 2.99 +0. 34 −0.33 for the 2018 dataset [36] Planck TT + TE + EE + lowE + Lensing + BAO at 95% confidence level. It also constrains the effective mass of an extra sterile neutrino m eff ν,sterile < 0.65 eV with N eff < 3.29 for the same dataset and same confidence level (though this value depends on chosen prior). 2 While within the paradigm of standard cosmology all these constraints impact on the viability of a light sterile neutrino with sizable active-sterile mixing, it has been shown that the CMB constraints can be partly relaxed going beyond the standard cosmology modifying the primordial power spectrum at small scales [37], within the paradigm of modified gravity [38], within BSM physics with time-varying Dark Energy component [39], by light dark matter [40], from large lepton asymmetry [41]. However, it has been also pointed out that additional interactions in the sterile neutrino sector can also render such a scenario viable [42][43][44][45]. The presence of such interaction leads to the suppression of the active-sterile mixing angle and thus delay the DW production. This significantly reduce the sterile neutrino abundance in the early Universe. In addition, it provides a mechanism to cut-off the free-streaming length of the sterile neutrinos at late time, and opens up an annihilation channel for the same. However, suppression of DW production alone does not suffice to constrain the energy density in the sterile neutrino sector. In addition, one also needs to assume that post-inflationary production of sterile neutrino and the light mediator, at least from the inflaton decay, remains small. During the re-heating epoch, considering perturbative decay of the inflaton, this can be ensured by simply assuming that the branching ratio of the inflaton into the sterile sector particles remain insignificant compared to the Standard Model (SM) particles. However, even this additional consideration does not serve the purpose when a bosonic mediator is invoked. The reason is that post-inflationary particle production can be significant during preheating [46]. While light fermions, which couples to the inflaton, are not produced in abundance during this epoch, the same does not hold for bosons. A (light) boson, which couples to the inflaton (via a quartic and/or tri-linear coupling, say) can be copiously produced during this non-perturbative process, thanks to the large Bose enhancement. Thus, while attempting to suppress the DW production of the sterile neutrino with secret interactions, the possibility of producing the bosonic mediator during (p)re-heating, and therefore, that of the sterile neutrinos can not be ignored.
In this article, we have considered a minimal renormalizable framework, consisting of an inflaton, the Higgs boson, and the light mediator (interacting with the sterile neutrinos) as only scalar particles to explore issues which are essential to make such a sterile neutrino sector cosmologically viable, starting from inflation. Within this framework all renormalizable terms are sketched out and their roles have been explored. Generally, the inflaton couples to the light mediator which can give rise to large effective mass during the inflationary epoch. Consequently, this prevents the light field to execute jumps of order H 2π , H being the Hubble parameter during inflation [47]. Thus, it ensures any additional contribution to the energy density from the light scalar remains negligible, which in turn, evades stringent constraint from non-observation of iso-curvature perturbation by CMB missions [48]. However, the same term can lead to the possibility of production during the pre-heating epoch. While several studies have considered the presence of a light scalar field with negligible coupling to the inflaton, and have put forward stringent constraint on the quartic coupling of the light field [49], aspects of the non-perturbative production, especially with a small quartic selfcoupling, during the pre-heating epoch have not been considered in details. This may lead to serious issues which may destroy inflationary cosmology altogether. In this article, we have explored the production of the scalar mediator during (p)reheating and subsequent production of ν s , explicitly stressing on regions of the parameter space where the production of the light mediator, and consequently, the light sterile neutrino can be significant at the on-set of BBN, and also the regions where such a sterile sector may be viable. We then discuss about some benchmark parameter values elaborating the same. 3 The paper is organized in the following order. In section:2 the model has been discussed. In the following section:3 constraints on the relevant model parameters from inflation, stability of the potential has been described. Further, the constraints already present in the literature on such secret interactions has been sketched. Subsequently, in section:4 we discuss the production of the light pseudoscalar during pre-heating and estimate the abundance of the sterile neutrinos in details. Finally, in section:5 we summarize our findings.

Construction of the minimal framework
We shall begin with a simple potential and introduce terms by justification from phenomenological reasons so as to arrive at a minimal model that can serve our purpose. To start with, we keep a pseudoscalar χ (to have a secret interaction between sterile 3 Note that, since within the SM, Higgs boson vacuum becomes metastable at an energy scale 10 11 GeV. A rather light Higgs boson, during inflation, can be displaced by H 2π , and can end up in a different vacuum. Several articles have considered this issue, and it has been shown that thermal corrections to the Higgs potential can possibly resolve this. We, for simplicity, would use the Higgs quartic coupling as a parameter during (p)reheating, and assume that, new physics (e.g. the presence of a heavy scalar, suitably coupled to the Higgs) contributions can raise the quartic coupling at such energy scale. neutrinos), and introduce a real scalar φ for the inflaton (χ or Higgs can not be used as the inflaton as will be discussed in Section:3.1.1). So, the potential in this model looks like where H is the SM Higgs boson and V inf is the inflation potential. Although this model does not seem to be phenomenologically economical due to the presence of multiple d.o.f, (multiple scalar fieds) but we have the independent parameters λ φH and λ φχ , λ H and λ χ which determine the energy flow to χ sector and visible sector during preheating without giving effective mass term to χ. The vacuum stability is an important requirement for successful preheating; the negative Higgs quartic coupling can make the potential unstable during the preheating stage. In the SM, the Higgs quartic term λ H becomes negative at around 10 10 − 10 11 GeV [50]. In our model, the energy scale during preheating is 10 13 GeV . Moreover, the energy flow to a sector depends on its quartic coupling, as we will discuss in Section: 4. Ameliorating the Higgs vacuum stability by coupling Higgs to inflaton requires λ Hφ ∼ O(10 −1 ) [51], thus ruining the flatness of inflaton potential via radiative corrections (as reflected in Eq: 3.3). In that scenario, one needs to use large non-minimal coupling to gravity to make the inflaton potential flat. On the other hand, keeping λ φH small does not help to improve Higgs vacuum stability. Addressing the stability through large λ χH ∼ O(10 −1 ) will invariably thermalise the χ sector with the SM and lead to a mass term to the χ field, m χ ∼ √ .1 × 246 √ 2 GeV ∼ 55 GeV after electroweak phase transition (EWPT). It is also constrained by the invisible decay width of the Higgs. Therefore, for the simplest case we assume that Higgs potential has been stabilized by another scalar and treat the Higgs quartic coupling as a free parameter during preheating.
Generically, as observed in numerical simulations, after preheating, the energy density stored in inflaton (kinetic and gradient energy) is of the same order with that of the other fields coupled with it. The final energy fraction transferred to a sector χ (considering a λ φχ φ 2 χ 2 interaction) is weakly dependent on the λ φχ coupling if the coupling is greater than a certain threshold value, although the initial rate of energy transfer is larger for a sector with larger λ φχ . Further, a quartic interaction term like χ 4 of a field can block the energy flow to that sector, as will be discussed in Section:4. Similar statements are also applicable for energy transfer to Higgs sector. So, during preheating, energy flows to SM and χ sectors equally along with itself (considering only φ 2 χ 2 and φ 2 H † H like terms, with additional control of the energy flow governed by the self quartic terms). Since back scattering φ particles to χ is not much effective for energy transfer from the energy density left in inflaton to some other sector [52], the possible way to flow this energy is to consider trilinear term(s). Inflaton-sterile neutrino coupling (of the form φν sνs ) does not help in the case we are concerned about, because it results in the total energy of the inflaton flowing into the sterile neutrino sector making the energy density in sterile sector and SM sector comparable, which ruins the N ef f bound at BBN. So we need the trilinear term with scalars (Higgs and χ).
Although it may cause instabilities in the potential, this instability can be taken care of by the quartic term of that field (as will be see in Section: 3.1.2). So the first step towards extension of the potential (2.1) is given by, where we have added the the trilinear mixing terms σ φχ and λ φH with dimension as that of mass. The inflaton decay rates arising due to the trilinear terms are straight-forward and are given by, Energy flow to any sector i by decay of the inflaton depends on the branching ratio defined by, 3 Cosmology and Light Sterile Neutrinos: Initial Constraints

Parameters of the scalar potential
As observed in Section: 2, there are independent parameters in the minimal potential required for a viable cosmological scenario. Though, apparently, there seems to be a lot of independence in the parameter space, in this section, we argue that the parameters, in no way can be of arbitrary values. Rather, they are tightly constrained to certain regions imposed by the inflationary paradigm as well as phenomenological requirements which in turn make the framework a minimal one.

Inflation, quantum corrections and threat to flatness of potential
Commencing on the exploring the constraints from the inflationary paradigm, in this setup we have two scalar fields, namely, the SM Higgs and χ. Inflation with the Higgs field is a well studied subject [53]. It has been shown that Higgs as inflaton requires a large non-minimal coupling of order ξ ∼ 50000, which can result in so called unitarity violation. But if the Higgs vacuum becomes unstable i.e. λ H becomes negative at high field values (which is certainly the case during inflation), the non-minimal coupling is of no help to drive inflation. As shown in [54] successful Higgs inflation can take place even if the SM vacuum is not absolutely stable, taking into account effective renormalization of the SM couplings at the energy scale M P l /ξ and symmetry restoration due to high temperature effects after inflation, which leads to the temporary disappearance of the vacuum at Planck values of the Higgs field. Realistically it turns out, Higgs is not quite a natural choice for inflaton candidate. So, to avoid large non-minimal coupling, and for simplicity we shun the Higgs as inflaton. Using χ as inflaton is also problematic, as to have enough energy flow to the Higgs sector (and thus to SM sector) we need the field χ to couple to the Higgs field, but that consequently leads to a mass term of the χ field, thus χ not being light enough to evade Σm ν bound [44]. In order to bypass the above issues related to Higgs and χ fields, we introduce a separate inflaton field φ. While choosing the model for inflation, we keep in mind that if the inflaton gets a vev at the end of inflation, the χ field gets a mass term. This leads us to consider large field models of inflation where the inflaton does not obtain vev after inflation. For simplicity here we consider quadratic and quartic inflation. However, driving inflation to produce the right amount of seed perturbation requires the inflaton potential to be very flat. In particular, in the large field inflationary models, e.g. [48] constrain the scalar power spectrum tilt n s ∼ 0.9670 ± .0037, as well as the tensor perturbation via tensor to scalar ratio r ≤ 0.065 for the 2018 dataset Planck TT,TE,EE+low E+Lensing+BK14+BAO. It has been observed that although the observations marginally satisfy quadratic inflation predictions, it is not possible to satisfy the n s −r bounds with a potential like λ φ φ 4 alone, as for 60 e-foldings we get a point in n s − r space outside the 2σ region of the Planck data. However, these observations are usually satisfied by coupling the inflaton φ non-minimally to gravity with the potential, Inflation constraints then can be easily satisfied with small values of ξ ∼ 10 −3 to 1 with λ φ ∼ 10 −13 to 10 −10 [55]. Similar strategy can be used for quadratic potential of inflaton with ξ ∼ 10 −3 for m φ ∼ 10 −6 M pl [56]. In this work, we do not vary this two parameters from the values mentioned above and also keep the values of other parameters such that these constraints are satisfied. Now, for a potential like, the Renormalization Group Equations (RGE) of the quartic interaction coefficients in this Lagrangian for 1-loop look like [57,58], From Eq: [3.3], it may be concluded that if the value of λ φ is O(10 −14 ) (for quartic inflation) at inflationary scales, the terms in the RHS of the Eq: [3.3] is needed to be less than or of same order of the initial value of λ φ to keep λ φ close to that order during the entire inflationary period. This means that the value of λ φH and λ φχ to be of order O(10 −7 ) at inflationary scales as (λ φH does not evolve much with energy as clear from Eq: [3.3]). These constraints on λ φH and λ φχ can be weakened to ≤ O(10 −5 ) by increasing ξ coupling to O(1).
It is shown in Section: 4 that a positive value of λ H is required for a successful preheating phase, and the energy flow to the Higgs sector during preheating explicitly depends on its value during preheating. The point clearly manifesting from the RG running is that the Higgs vacuum stability can not be improved with the coupling with inflaton keeping the flatness of inflaton potential because then the large inflatonhiggs coupling can change the λ φ considerably. Also, λ χH cannot be used to improve Higgs stability as one needs to keep it small so that (i) χ does not get thermalized with the SM, (ii) Higgs vev does not induce high mass to χ and (iii) constraints from Higgs invisible decay are satisfied [59]. On the other hand, if the Higgs stability is addressed with inflaton by the modified RGEs, the value of the parameter λ φ can be stable throughout the period of inflation if it is at O(10 −4 ) during inflation. In that case a large non-minimal coupling is required to drive inflation with the observational constraints satisfied.
Along with these arguments mentioned above, it is interesting to study two different cases: • The Higgs stability is addressed by some other scalar field (see for example [60]) and the couplings λ φH , λ φχ are small to keep the λ φ stable throughout the period of inflation. Then we study the preheating of χ and H, inflaton decay into χ and H and produce ν s by back-scattering and oscillation from active neutrinos.
• The Higgs stability is addressed with inflaton by the modified RGEs (see for example [51]). To do that we need λ φ ∼ O(10 −4 ) and the inflation is driven with help of large non-minimal coupling to gravity. In this case the couplings λ φH , λ φχ can be large.
For the first case, this will result in a positive value for λ H at inflationary scale driven by the extra BSM scalar coupling to the Higgs. So this particle will be thermalized with Higgs and the SM sector and energy flow during preheating will be affected. But the ratio of ρχ ρ SM as shown in section: 4 will essentially remain the same due to extra BSM scalar field not coupling to χ and our analysis and results will remain the same. To summarize we assume that some Higgs vacuum is stabilised by one such mechanism and remains so during inflation. Thus the assumption of positive λ H is justified.
Taking a large non-minimal coupling changes the dynamics of the field considerably, so even near the minima of the inflaton potential we can not neglect the effects of the non-minimal coupling. So, it may change the results of the (p)reheating dynamics, which is dependent on the dynamics of inflaton near its minima, from the case with no non-minimal coupling. So, for simplicity of numerical calculation, in the present work we take the first route. We choose quadratic inflation for our scenario over quartic inflation for the reason discussed in Section: 3.1.2.

Handling the problem with vacuum stability
If inflation is driven by quartic potential along with a trilinear term involving inflaton and another field (e.g. χ) it gets vev at the end of the inflation. This is problematic for model-buidling perspective as the vev of inflaton and χ result in mass terms for χ and ν s respectively but we want the particles χ and ν s to be of small masses O(eV ). So, we would need extreme fine tuning in this case. On the other hand if the inflaton has a quadratic term then it is possible to have the minima of the potential at 0 field values. This is why we mainly consider quadratic inflation for our case. Nevertheless, even if we choose to ignore the quartic term at some energy scale (i.e. set it to 0), it will become nonzero at other energy scales due to RGE running. So, the final minimal potential that can serve our purpose is given by 4 : In order to have stability of the potential we only need to check if the potential is bounded from below or not. But, in the present scenario we also need the minima of the potential to be at (0,0,0) for small mass of the χ particle without any fine tuning. As the potential at (0,0,0) is 0, to have (0,0,0) as the minima, we need the potential to be non-negative and rewrite the potential as sum of positive terms: From this equation it is clear that the minima of the potential is at (0, 0, 0) in field space (here we are talking about scenarios just after inflation and so we neglect the EW vev of Higgs) under the conditions λ χ > given two arbitrary real constants α, β with α 2 + β 2 = 1. From these conditions, we can see that along with the quadratic inflation with m φ ∼ 10 −6 , an equally possible scenario with small m φ and small σ φH is possible with the inflation being a quartic one. In that case too, the energy left in the inflaton will have to be decayed, with additional benefit that the energy density left in the inflaton evolves as radiation and we do not have to be worried about the decay time scale of inflaton after preheating. But in the (p)reheating numerical simulations, we do not consider such cases and keep m φ ∼ 10 −6 and λ φ ∼ 10 −14 for all the simulations to understand the effects of the other parameters clearly.

Iso-curvature perturbations and stability of light fields during inflation
The relaxation time scale for a quantum fluctuation to roll back down to its minima is m −1 (during inflation, a field coupled to inflaton has mass depending on the expectation value of the inflaton φ 0 , if its bare mass is negligible i.e. light inflaton), whereas the time scale for the evolution of the universe is given by H −1 . So, if m > H, then the field is stable and the curvature perturbations due to that field may be neglected. This condition translates to ( for quartic or quadratic inflation. Keeping the couplings substantially larger (10 3 to 10 4 times) than λ φ or m 2 φ can stabilize those fields during inflation. This is a condition we keep in mind while choosing the parameter values.

Interaction parameters and bounds on m χ − g s plane
For the sake of completeness of our discussion on different constraints on the parameters coming out of different physical requirements, let us briefly summarise the existing constraints on the interaction parameters. After the (p)reheating is over, in principle the entire evolution of the spectrum(s) of the species is governed by the Boltzmann equation: Here L is the Liouville operator and C is the collision operator. We are interested in finding the spectrum of the sterile neutrinos through collision operators corresponding to χχ −→ ν s ν s and oscillation from active neutrinos. Solving the Boltzman equation with the entire spectrum f i is difficult to solve even numerically, so it is assumed that the χ particles should follow a thermal distribution, i.e. χ particles are thermalized among themselves. For this assumption to be true just we need a parameter region of λ χ estimated by the relation: where Γ is the interaction rate, H is the Hubble parameter, σ is the interaction crosssection, v mol is the Moller velocity of χ and n is the total number density. For our model, during radiation dominated epoch, for χ χ → χ χ scattering, n χ = 3 4 ζ(3) whereas the Hubble parameter is given by As temperature of any relativistic species goes down at same rate 1/a, even if the temperature of χ and SM are different, Γ eventually becomes lower than H and gets thermal distribution. An estimate of the thermalisation temperature shows λ χ ∼ 10 −8 gives T Ther ∼ 10 −16 M P l (assuming T SM ∼ T χ , as even much lower energy density means temperature difference of order (density) 1/4 ), which is far before BBN.
The thermalization process within some sector starts much before the interaction rate (calculated from scattering of the particles) becomes comparable to the Hubble rate. It is well known that at the start of the preheating epoch, modes with only some specific wave numbers gets excited exponentially governed by Mathieu equation. But, it has been observed from the LATTICEEASY simulation that, even if at the start of preheating stage, only some specific infrared momentum modes gets excited, as time progresses, the energy gets distributed to higher momentum modes. This observation can be interpreted as start of thermalisation process at the end of preheating [61].
The thermalisation of χ and ν s is governed by the interaction χ χ → ν s ν s , having . This means for g s ∼ 10 −4 the thermalisation happens at 1 GeV [43].
Bounds on m χ − g s plane (i) From BBN: The standard way to parameterize the radiation energy density (ρ R ) is like [62], where ρ γ is the photon energy energy density, N eff is the effective number of relativistic species. A universe with only active neutrinos lead to N eff = 3.046. Any extra radiation-like component (light sterile neutrinos), if present, will contribute to this N eff . BBN observations constrain this value of N eff , we use a conservative bound of N eff < 3.5 at 68% confidence level [28]. In the density matrix formalism the differential equations governing evolution of the neutrino density are (for details refer to [63][64][65][66]) where ρ is the 4 × 4 neutrino density matrix in the flavor basis with diagonal entries corresponding to physical densities, H is the full Hamiltonian, H m = U H 0 U † is a rotation of the free neutrino Hamiltonian in the mass basis H 0 = diag(E 1 , E 2 , E 3 , E 4 ), and ρ eq is the density matrix at thermal equilibrium, ρ eq = I 1/ 1 + e E/T and V eff is effective neutrino thermal potential in presence of the background χ field. The presence of V eff (see Appendix for detailed calculation) leads to suppression of the neutrino production through MSW-like effect via the χ particle through the effective mixing angle θ M . Note that we need a χ background to have the suppression, so the χ particles should not decay before the BBN. Though the decay time scale of such a particle into ν s or ν particles in the rest frame of χ is very small, the time dilation saves them from decaying when they are relativistic. To approximate an analytic solution for the evolution equation for f s , the phase space distribution of ν s , in terms of the interplay of oscillations and collisions, we start with the equation, From this equation, following [67], the contribution of ν s to ∆N BBN can be found to be We expect the analytic expressions to be less accurate than numeric solutions mainly due to the discrepancy in g * , which is kept fixed in the analytic solution. The results achieved from analytics nearly resemble the full numerical results by solving quantum kinetic equations as in [43]. The Blue region in Fig. 1 corresponds to the allowed region in m χ − g s plane from N eff BBN constraints.
It should be noted that if the sterile neutrino production is suppressed until the active neutrinos decouple from the SM soup, then production of sterile neutrinos from oscillation increases the d.o.f in the decoupled neutrino sector. This, in turn, decreases the energy density stored in this sector by a factor (g 1 /g 2 ) 1/3 = (21/32) 1/3 = 0.87. It is the opposite effect to that when the SM electron positron annihilates to photons resulting in the decrement of d.o.f in the sector and the energy density per unit comoving volume increases. So, in our case, the N eff decreases by a factor 0.87 as it is scaled to the energy density of the SM photons (which does not change due to the d.o.f change in the neutrino sector) as per Eqn. 3.8.
(ii) From CMB & LSS: Hitherto, the only requirement is to have suppression in the neutrino oscillation which consequently reduces sterile neutrino production until BBN from the early inflation and reheating epoch. However the portion of the parameter space allowed from BBN may not be compatible with CMB and LSS constraints. The physical condition for the observed power spectrum in CMB and LSS is not only to have the active neutrinos free-streaming, but also another sterile neutrino species (with O(eV ) mass) to be free-streaming. If this new species is of with similar number density as that of the active neutrinos, then there is much suppression in the power spectrum than that observed in CMB or LSS. So, to satisfy the CMB and LSS constraints of Σm ν (which basically quantifies the total mass of free-streaming species, assuming the same number density as that of active neutrinos), any massive species must annihilate or decay into lower mass particles to evade the mass constraints all together or they must interact among themselves or with some other species in order to cut off the free-streaming length scale. Ref. [68] followed the second route, i.e., they used strong interaction strength for the sterile neutrinos to self-interact via a secret mediator to cut off the free-streaming length. But, recent studies [69] have shown that this scenario generates interactions between the active neutrinos too, through flavor mixing (note that the suppression in oscillation is lifted off for most of the parameter space for eV scale and below), leading to a higher amplitude in power spectrum than the vanilla model itself. So, having a large interaction strength (with help of large coupling and/or small mediator mass) is perilous for the cosmological observables, if one relies only on cutting off the free-streaming length. Instead, [44] used the first principle to evade the mass bounds coming from CMB and LSS observations. They chose the mediator mass to be O(< 0.1 eV ) which meant that the interaction ν s ν s → χ χ goes only in the forward direction once the temperature of the Universe goes below O(1 eV ), i.e. the backward interaction becomes kinematically inaccessible due to Hubble expansion dominating over this process. Thus, for this suitable choice, the free-streaming length scale (as per CMB and LSS) need not be cut off since the sterile neutrinos annihilates into χ particles only, thereby not hurting the mass bound. The gray region to the right of m χ = 0.1 eV in Fig. 1 shows the bound from CMB and LSS.
(iii) From Supernova: Constraints from SNe observations [70,71] also does not allow the couplings to be g s > 10 −4 . So, this leaves us with a patch of parameter space at the lower left corner of the BBN allowed region as depicted in Fig. 1.

Production of Sterile Sector from (P)reheating
After the inflation is over, all the energy density is in the inflaton field, this energy density flows to other sectors by mechanisms like (p)reheating. In this epoch, there are two stages of evolution: • Initial stage that can be treated mostly analytically • Back-reaction dominated stage for which one needs to do a detailed numerical treatment.
In what follows we shall analyze the two stages separately in order to get a complete picture.

Initial stage of (P)reheating
In this section we try to understand the physics behind the energy flow from the inflaton to the other fields coupled to it by preheating. For simplicity we first consider the evolution of the inflaton (after inflationary period is over) neglecting the couplings to other fields. The evolution of the zero mode of inflaton is governed by, For quadratic potential it has sinusoidally oscillatory solution with decaying amplitude (due to Hubble friction term).The dynamics of the Fourier modes of a field χ coupled to the inflaton in FRW universe is given by [46], This can be written (neglecting expansion) as the well known Mathieu equation, , z = m φ t and prime denotes differentiation with respect to z. Mathieu equation has well known unstable exponential solutions for instability regions of A − q parameter space and hence for specific k values. Modes corresponding to these k values grow exponentially, which in physics is interpreted as exponential particle production in those modes.
The statements upto now is only true for the initial stages of preheating. The growth of the fluctuations give rise to a mass term λ φχ χ 2 to the e.o.m of the zero mode of the inflaton and as a result affects other modes through Eq: 4.1. This phenomenon is known as the back-reaction effect in the literature, after the back-reaction effect starts, the simple analytic expression Eq: 4.1 does not hold true. The back-reaction effect is usually estimated by the Hartree approximation, but still it does not take care of effects like re-scattering. So the only way to fully solve the e.o.m of the fluctuations throughout the preheating period is by lattice simulation. These simulations solve the classical field equations in lattice points numerically and give far accurate results than approximate analytic solutions.

Numerical evolution for Back-reaction dominated stage
To calculate preheating numerically on the lattice we use the publicly available code LATTICEEASY [72]. In this section we do case by case study beginning from the simplest potential. In each case we demonstrate some benchmark values of the parameters concerned, and study whether the case with the chosen values of parameters is viable with cosmological observations or not. As we proceed we add new terms to the potential and discover the implicit assumptions needed to be assumed to reconcile the models with cosmological observations. It is to be noted that, we neglect the effect of non-minimal coupling to gravity during the (p)reheating era. This is a logical assumption, as the potential remains unchanged near the minima of the inflaton, for small non-minimal coupling, and the Fourier modes of the species get exponentially excited only during the time the inflaton passes its minima.
We start with the potential, In this case we assume λ χH to be negligible, lest it thermalise the SM sector with that of the sterile. The energy density in the inflaton fluctuations, Higgs and χ field at the end of preheating is of the same order if λ H and λ χ are small with respect to λ φH and λ φχ . We have found from our simulation the energy-flow to a sector is not only dependent on the coupling of the field with inflaton (if the coupling is lower than a certain threshold value (see Ref. [73] for details), the energy flow to that sector is highly suppressed, whereas for values of that coupling over the threshold, the amount of energy flow is weakly dependent on the coupling , but also if the field has a quartic self-coupling, it blocks the energy flow to that sector, as clear from the plots (Fig:2). These salient features are also in agreement with those studied in Refs. [73,74]. The reason for this phenomenon is possibly because of the extra energy cost due to the potential term, which blocks the modes to grow. It can also be interpreted from a different point of view: as the Fourier modes of χ grow, the quartic behaves like an extra mass term 1 2 λ φχ χ 2 χ 2 , making it difficult for the particle to be produced. This piece of information is very vital for our calculation as we control the flow to χ sector by these quartic terms.
Significance of the trilinear coupling: If there is no trilinear coupling of the fields to inflaton, the inflaton cannot fully decay to other fields thus contributing a significant amount to the total energy density of the Universe. This is also shown in Ref. [52]. If a trilinear coupling of sterile neutrino to the inflaton, namely, g φνs φν s ν s is present, the total energy left in inflaton after preheating will be transferred to sterile neutrinos leading to similar energy densities in the SM and the sterile neutrinos. So this case is trivially not satisfied by the N ef f constraints. Therefore, in order to direct the desired flow the energy density stored in inflaton to other sectors, we introduce trilinear couplings of inflaton to other scalars. Below we will discuss the possible scenarios case by case.

Trilinear interactions of inflaton with Higgs only
First we discuss a case where inflaton has trilinear interactions with Higgs only. In this case the total energy left in inflaton after preheating will decay to SM. Even this is not be compatible with cosmology if after preheating the order of energy densities in inflaton, Higgs and χ is same. To have a viable scenario we must have a way so as to block the production of χ sector -we do it by using what we observed in the preheating section -a self-quartic coupling in χ sector is good enough for this matter. This is one of our main results where we show some benchmark values of mainly λ χ with respect to λ φχ with all the other values fixed, where the N eff bound described in Section: 3.2. We see that increasing the value of λ χ with respect to λ φχ increases the blocking of energy flow to χ sector, eventually resulting in lower energy density in this sector [Fig: 2]. Even if there is comparable energy density in the χ, Higgs and inflaton sectors after preheating, the final relative energy density of the sectors depends on the state of inflaton during its decay (that is, if it is relativistic or non-relativistic). In Fig:  4 we show some allowed benchmark values which satisfy the N eff constraints at BBN. Note that in this case we have trilinear coupling of inflaton only to the Higgs, making sure that decay of inflaton does not populate χ sector.
The values of λ χ and λ H are the parameters we may vary to have a grip on the energy flow fractions in different sectors. Taking the values of λ φH and λ φχ too large can ruin the flatness of inflaton potential as we saw in, whereas making them much smaller can make their contributions during preheating sub-dominant enough with respect to contributions from trilinear terms (as shown in the section: 3.1.1). We have checked that keeping the λ H and λ χ same whilst λ φH and λ φχ are unequal makes the energy flow to the Higgs and χ sector as expected, i.e. the sector with higher λ mix will get larger share of energy density [Fig: 6]. After preheating, the energy left in the inflaton is transferred to Higgs sector through its perturbative decay due to the σ φH coupling. The inflationary model, if chosen to be quartic, the energy density stored in the inflaton evolves as radiation. Whereas for the quadratic case, the process is complicated as the energy stored in the inflaton (behaving as a condensate) should behave as matter and this can dilute the energy densities produced in preheating. However it has been observed that, if a trilinear decay term is present, the equation of state behaves more like radiation than matter after the preheating stage [52]. The decay of the inflaton typically happens much after the preheating epoch and the dynamics of this stage is not well known. So, for simplicity we extrapolate the radiation-like feature of inflaton and consider that the inflaton becomes non-relativistic only after the temperature of the Universe becomes comparable to the mass of the inflaton. In our analysis, we keep the values of σ φH of the order ∼ 10 −10 M P l and/or ∼ 10 −8 M P l . The choice of σ φH ∼ 10 −8 M P l corresponds to the case when the T R ∼ m φ . Whereas, if we decrease σ φH , the decay of inflaton happens after it becomes non-relativistic, which means a matter dominated phase partly washes out the preheating contributions to the relativistic χ and Higgs and consequently, their final energy density, depending on the timespan during which the inflaton is non-relativistic. Then the final energy density tends to solely depend on the branching fraction of inflaton, which is a trivial result. On the other hand, as shown in section: 3.1.2, increasing σ φH and σ φχ means requires the increment of λ χ and λ H values as well, which subsequently will lead to more blocking of corresponding energy density flows during preheating, thereby making the final energy fractions dependent only on the branching ratios of the inflaton.

Trilinear interactions of inflaton with χ only
Inflaton with trilinear interactions with only χ is trivially not satisfied by the N eff constraints as in this case the total inflaton energy after preheating flows into the χ sector.

Trilinear interactions with both Higgs and χ
In this case inflaton has trilinear interactions with both Higgs and χ. At first during preheating both Higgs and χ are produced, and energy density ratios depend on the self quartic couplings of those sectors. Then depending on the mass of the inflaton and the trilinear couplings, the inflaton decays into Higgs and χ (and also possibly ν s ) according to their branching fractions. In Fig:7 we show the fraction of energy density of χ and Higgs during preheating when inflaton has same trilinear coupling to both χ and Higgs. We observe that for the parameter values we used, the energy flow fractions do not change much, even with the trilinear term present in the χ sector; this is due to the fact that the production during preheating is sub-dominant from to the trilinear term. But even if the fraction of χ energy is small after preheating, due to the trilinear term, the inflaton decay to χ channel is open and for the chosen σ φχ = σ φH = 10 −10 M P l value, the branching ratio is same for both the fields (neglecting m H and m χ to be negligible with respect to m φ ). This case is not allowed by N eff bounds of BBN. Varying the ratio of σ φχ and σ φH changes the branching fraction enabling the Higgs and the χ sectors to be populated unequally during inflaton decay. But, if we do not want an epoch when the inflaton behaves as a non-relativistic species (i.e., when the history of preheating dilutes away), then we would like to have large σ φH and σ φχ (O(10 −8 ) M P l ).
In Fig: 5, we show how the decay channel of the inflaton into χ can change the N eff values. However, in that case, the H and χ sectors will thermalise with each other resulting in a fully thermalised species with the SM, thereby trivially not respecting the N eff bounds of BBN.

Allowed benchmark parameter values and additional constraint on
In our scenario, since we assume m χ < 2m νs , which is certainly our scenario, the primary production channel of the ν s particles is from the χ via back-scattering particles since that from active neutrinos through oscillation is suppressed . The ν s production from backscattering depends on the n σv mol of the interaction in comparison to the Hubble parameter. The thermally averaged cross-section for the process χχ → ν s ν s is given by (in the relativistic limit) [75], So, the χ and ν s sectors thermalise at T ∼ 1 GeV for g s ∼ 10 −4 . The extra N eff coming from the entropy of χ sector, and from the ν s sector, both of which are now partly or fully thermalised with each other, changes the bound in m χ −g s plane from BBN which was previously considered in the case for production of ν s only from active neutrinos when inflation was not considered in the whole picture.. It is clear from Fig: 4, that in several regions of the parameter space, the N ef f bound is violated at the time of BBN, just from χ production during preheating. For the values where the bound is satisfied, we get additional constraints in the m χ − g s plane as shown in Fig: 8 as an example.

Conclusions and Outlook
In this article we investigated for the possibility of having an extra sterile neutrino in the particle spectrum. It is well-known that it is possible to reconcile the cosmological observations with this extra species, if we introduce a pseudoscalar interacting with the sterile neutrinos. As a background field, the pseudoscalar particle creates an effective thermal neutrino potential, which, due to its matter-like effect, suppresses the DW-like production of sterile neutrinos. Such a BSM scalar, if present in the early universe, should also be produced inadvertently during preheating due to its coupling to the inflaton. Moreover, Bose enhancement will make this production copious enough such that its primordial abundance and that of sterile neutrinos cannot be assumed to be negligible with respect to SM particles for most scenarios. This assumption was the cornerstone of the scalar or vector interactions that alleviate the sterile neutrino constraints from cosmology. Even though the production during reheating can be neglected by considering small trilinear coupling as the inflaton decay is straightforward and depends solely on the branching fraction, the production from preheating is non-trivial. Thus, it is important to consider the production right from preheating, which is a highly non-linear process after initial exponential growth; this was studied in exquisite details in this article. We made use of analytical arguments and numerical calculations using LATTICEEASY simulation to find the regions of the parameters where this production will be significant. Consequently, we arrive at the following conclusions: • λ χ needs to be kept large to suppress the preheating production.
• λ χH needs to be negligible to prevent SM from thermalizing with the BSM sector.
• λ φH and λ φχ cannot be large, or else the flatness of inflaton potential will be ruined due to RG running.
Thus, we built up a minimal model that can take into account possible issues of (p)reheating, commenting on the Higgs stability and, finally, re-laid the existing observational observational bounds from BBN, CMB and LSS in context to the parameters relevant to (p)reheating dynamics. It turns out that Higgs vacuum stability cannot not be improved by the fields and the configuration taken into consideration; we assumed it to be improved by some other scalar field or other mechanism in that context. It should be noted that the Higgs quartic coupling is required to be positive during preheating for a stable (p)reheating dynamics. We conclude that the pseudoscalar abundance from preheating can be high enough to violate the bounds of N eff from BBN, if the self quartic coupling is not high enough and/or there is no period before decay of inflaton when it is non-relativistic. It is to be noted that even if these two conditions are satisfied, the inflaton can still decay to the pseudoscalar along with Higgs and the BBN N eff bound can be at risk. We found benchmark values of the model parameters, for which the initial neutrino abundance of the pseudoscalars and sterile neutrinos are favorable for a viable cosmological scenario, and discussed the impact of the various parameters on the final abundances of the pseudoscalar and SM particles. It turns out that the for a small non-minimal coupling to gravity, an inflation mass (m φ of 10 12 GeV) compatible with recent n s −r observations of Planck 2018, and small self-quartic of inflaton (λ φ ∼ 10 −14 ), we have the suitable parameter space of λ χ ≥ 2 × 10 −5 for mixing parameters λ φχ and λ φH of the O(10 −6 ). The trilinear term in the potential, σ φH was kept of order O(10 −8 ) M P l chosen such that the inflaton decays when the temperature of universe is of the same order of the m φ , λ H of order O(10 −4 ). In the neutrino parameter plane, we placed strong constraints, that is, for coupling strengths g s ≤ 10 −5 and sterile neutrino mass m χ ≤ 10 −1 eV all the cosmological bounds are satisfied. Thus the sterile neutrino with the pseudoscalar "secretinteraction" model can be an viable possibility when all the early universe cosmology is considered on one hand, and, the existing neutrino anomaly is invoked on the other hand.
Future CMB missions like LiteBIRD [76], COrE [77], PIXIE [78], CMB S4 [79], CMB Bharat [80] aim at constraining the inflationary observables and the other cosmological parameters further and hence will be an important probe for this model. Results from the neutrino experiments MicroBooNE [81] may decidedly prove the existence of the sterile neutrino. Sterile neutrinos with secret interactions have been proposed to be looked for in IceCube experiment [82]. In CMB polarization observations of BICEP the sterile neutrinos may also be a relevant signature to look for as per Ref. [83]. Finally primordial gravitational waves from inflation and preheating will be a very important signal (See Ref. [84] for a review) for the model. Traditionally, sterile neutrino as a dark matter candidate has also been considered (see [85] for a recent review) and see [86] for a recent approach using SIMP [87,88]. A single model with precise predictions is able to explain some of the central problems today and be tested in laboratory experiments and cosmological observations. We hope to take up some of these issues in future.

Suppression of ν s production from active-sterile oscillation
The evolution of spectrum of ν s , f νs (E,t) is governed by, assuming equilibrium distribution for active neutrinos,

∂ ∂t
− HE ∂ ∂E f νs (E, t) = C χχ−→νsνs + 1 2 sin 2 (2θ M (E, t)Γ(E, t))f a (E, t) (7.1) where Γ is [89,90]: for sterile neutrino redistribution through Fermi-Dirac distribution. The latter can be neglected since this is very small. C χχ−→νsνs is the collision term corresponding to χ annihilation given by: 4 δ E δ p f eq q f eq q (7.4) Now there maybe two cases one where the particles are already thermalised, so that one may take f eq q = e − q T and in another case where they are not in thermal equilibrium and need to be numerically evolved from preheating dynamics. The second term in (7.4) corresponds to an oscillation term, θ M being the mixing angle which is suppressed by introduction of the hidden sector interaction through a effective potential V eff as, sin 2 (2θ M ) = sin 2 (2θ 0 ) cos(2θ 0 ) + 2E δm 2 V ef f 2 + sin 2 (2θ 0 ) (7.5) This phenomena is basically the neutrino oscillation while propagating through a thermal heat bath filled with χ particles. The sterile neutrino self-energy at one-loop in a thermal bath is given by: Here, m is the sterile neutrino mass, p is its 4-momentum and u is the 4-momentum of the heat bath and is taken to be u = (1, 0, 0, 0) in its rest frame. The energy dispersion relation in the medium becomes: in the UV regime, which gives us: with corresponding η f (p) and η f (b) are the Fermi-Dirac and Bose-Einstein distribution occupation numbers respectively. Using [64,92] first delta function integrals are done to do p 0 part. Then the 3-momenta p-integral is shifted to spherical co-ordinates which will reduce to the standard p integral which is to be performed numerically.
We give some analytic estimates for the integral for low temperature limit:

32E
(T χ , E >> m χ ) (7.14) And similarly, for tadpole diagrams, The origin of the Eqn. (7.5) comes from the secret "pseudoscalar" interaction which introduces a matter potential causing MSW-like effect for sterile neutrinos of the form [93,94]: where f φ is the Bose-Einstein distribution for the pseudoscalar and f s is the distribution for the sterile neutrinos. The potential V s (p s ) is basically the thermal contribution of the background field in the form of bubble diagrams; an order-of-magnitude estimate of it goes as: V s ∼ 10 −1 g 2 s T. (7.17) considering a common temperature T for the all the species. This is the central idea behind the phenomenology that this matter-effect would induce a mixing angle different from that of the standard ν s → ν a (2-flavor approximation) and stop it from thermalizing with the SM. This can be alternatively looked upon as a minute shift in the effective mass-difference between the neutrino states. For a 2-neutrino framework, the thermalization process can be treated easily by the density matrix formalism leading to solving the Quantum Kinetic Equation (QKE) in equilibrium: where f 0 is the Fermi-Dirac distribution function. The QKEs are noẇ P a = V x P y + Γ a [2 − P a ] , P s = −V x P y + Γ s 2 f eq,s (T νs , µ νs ) f 0 − P s , and the potentials are: where p is the momentum, G F is the Fermi coupling constant, M Z is the mass of the Z boson, and n νs = f s d 3 p/(2π) 3 is the number density of sterile neutrinos. The range of the values of coupling g s for which N eff varies from 1 to 0 is g s ∼ 10 −6 to 10 −5 [43].