Sterile Neutrino Dark Matter from Generalized $CPT$-Symmetric Early-Universe Cosmologies

We generalize gravitational particle production in a radiation-dominated $CPT$-symmetric universe to non-standard, but also $CPT$-symmetric early universe cosmologies. We calculate the mass of a right-handed"sterile"neutrino needed for it to be the cosmological dark matter. Since generically sterile neutrinos mix with the Standard Model active neutrinos, we use state-of-the-art tools to compute the expected spectrum of gamma rays and high-energy active neutrinos from ultra-heavy sterile neutrino dark matter decay. We demonstrate that the sterile neutrinos are never in thermal equilibrium in the early universe. We show that very high-energy Cherenkov telescopes might detect a signal for sterile neutrino lifetimes up to around 10$^{27}$ s, while a signal in high-energy neutrino telescopes such as IceCube could be detectable for lifetimes up to 10$^{30}$ s, offering a better chance of detection across a vast landscape of possible masses.


I. INTRODUCTION
Dark matter (DM) continues to stand as one of the greatest mysteries at the interface of particle physics and cosmology [1]. Since no particle in the standard model (SM) can account for the bulk of the cosmological dark matter, one must extend the particle content of the theory to include at least an additional particle species, or a bound state thereof, that can account for a matter species making up around five times the baryonic matter. Adding right-handed neutrino (RHN) states is a well-motivated and oft-followed route that can at once address the origin of active neutrino masses, explain the matter-antimatter asymmetry, and provide a viable, arguably minimal, DM candidate [2].
Recently, Ref. [3,4] considered a peculiar and interesting possibility for the generation of the cosmological dark matter: while in flat space-time the "natural" vacuum is the one respecting the isometries of Minkowski space, in a general, curved space-time this is no longer true. The choice of vacuum is generally observer-dependent (for instance to inertial observers at different points in the space-time), and different observers will define different, inequivalent vacua. As a result, the zero-particle state according to one observer can actually be populated by particles according to a different observer (this is well known in the context of Hawking radiation from black hole evaporation, or the Unruh radiation experienced by an accelerated observer) [5,6]. In particular, this holds generally for Friedmann-Robertson-Walker (FRW) cosmologies. Ref. [3,4] realized that assuming time-reversal invariance allows to identify uniquely a CP T invariant vacuum in a FRW background, enabling a well-defined calculation of the number density of gravitationally produced particles as seen by a late-time observer.
Gravitational particle production in an FRW background without inflation, similar to the scenario under consideration here (albeit, here without the additional requirement of CP T invariance), has been studied extensively in the past (for reviews, we refer the Reader to Ref. [7,8]). The key theoretical results are that particle production (specifically fermions, or conformally coupled scalars) occurs at the cosmological epoch when particle masses are comparable to the Hubble expansion rate; the particle number density n ∼ m 3 , and particles are created with the equation of state of dust; finally, the fractional relic density of these particles X at the time of radiation-matter equality was found to be Ω X ∼ (m X /10 9 GeV) 5/2 [9]. Further studies considered similar scenarios, but in the presence of inflation (see e.g. [10] and references therein).
If the early universe can be described by a conformally flat FRW metric, the early-time and the late-time (or the 'in' and 'out') vacuum states are not equivalent due to the lack of time-translation symmetry in the metric. Calculating the number density of particles then involves performing a Bogoliubov transformation on the quantum field modes and finding the appropriate Bogoliubov coefficients. Remarkably, Ref. [3] found that the standard in/out vacuum states do not respect CP T , and constructing a CP T invariant vacuum state leads to a certain unique form of the relevant Bogoliubov coefficient [3,4]. Given the number density produced by this expansion mechanism, if this species were to account for the total dark matter (DM) density ρ DM , then its mass would have to be around 4.8×10 8 GeV.
This calculation was performed in the context of an early universe cosmology described by a radiation-dominated cosmology, with a scale factor a(τ ) ∝ τ , where τ indicates the conformal time, and for the most natural and minimal DM candidate, the right-handed neutrino, assumed to be absolutely stable.
Here, we note that the requirement of CP T symmetry is not unique to radiation domination: as we explain in this work, the universe could have been dominated, at early times, by some other species with a different energy density, and yet satisfy CP T symmetry. Note that the gravitational production of particles discussed here does not touch upon the question of whether the issues of homogeneity and isotropy are addressed or not, and only relies on CP T symmetry at positive conformal times [3,4].
The possibility of an alternate cosmological history prior to Bug Bang Nucleosynthesis has received considerable attention (for a recent review see Ref. [11]), and here, as far as our motivation, we build largely on the literature on the topic. The key observational fact for radiation domination is the synthesis of light elements; prior to that, there is no observational evidence that underpins radiation domination. In fact, the possibility of an early matter-dominated phase is very well justified and has been widely considered [12]; other causes of nonstandard early-universe expansion history include single or multi-field post-inflation reheating , multi-step or thermal inflation (in fact the presence of any extra stages of inflation, see e.g. Ref. [13] for a review), heavy particles and dark sectors [14], moduli fields [15]; the resulting structure of the cosmological expansion history impacts, among other things, gravitational particle production, as was noted in a recent study [16] that considered the case of quintessential cosmology in the early universe, with a formalism not dissimilar to the one employed here. Note, however, that the present setup is uniquely different from previous work where gravitational production is connected to inflation.
A non-standard FRW background leads to different equations of motion for the right-handed neutrino field, which, we note here, can be approximated using the WKB method. Once the appropriate Bogoliubov coefficient is computed from the WKB solutions in the asymptotic future, one can follow a similar analysis to Ref. [3] and calculate the resulting particle number density at late times. The final result of this exercise is a DM mass as a function of a parameter which indicates the scaling of the energy density with the cosmological scale factor and, as we discuss below, a second parameter that depends on when the universe transitions from the non-standard early-times cosmology to radiation domination, needed from cosmological observations to occur prior to the synthesis of light elements [1].
In this work, besides generalizing the results of [3,4] to any CP T -symmetric early-universe cosmology, we relax a key assumption in the dark matter sector: that the right-handed neutrino dark matter be absolutely stable because of a discrete symmetry. Relaxing such assumption produces the interesting prospect of testing experimentally this production mechanism by searching for the decay products of massive right handed neutrinos, with a lifetime related to the mixing angle between sterile and active neutrinos. In turn, this calls for addressing the question of whether the gravitationally-produced neutrinos can ever thermalize, which we address in detail in sec. IV.
The remainder of this work is structured as follows: in the following section II we describe in detail the calculation of the dark matter number density generated in arbitrary CP T symmetric cosmologies in the early universe; the following sec. III discusses the right-handed neutrino sector and the decay modes thereof; sec. IV explores the issue of thermalization of right-handed neutrinos; sec. V addresses constraints and prospects for the detection of decaying massive right-handed neutrinos as predicted in the present scenario; finally, sec. VI concludes.

II. STERILE NEUTRINO AS DARK MATTER IN GENERALIZED CP T SYMMETRIC COSMOLOGIES
Consider a cosmology dominated by a species with equation of state P = wρ, and thus with an energy density-scale factor relation Recalling that, for w = 1, time t ∼ a 3(1+w)/2 , and that conformal time is defined as the relation between conformal time and scale factor for an equation of state parameter w is The requirement for CP T symmetry is that a(τ ) = −a(−τ ); this, in turn, requires As a result, we obtain the equation of state of radiation (w = 1/3) for N = 1, −1/3 < w < −1/9 equations of state for positive N = 2 (we disregard the w = 0 case), and −1 < w < −1/3 for negative N . We do not concern ourselves with the issue of, and the specific model for, how the transition between the era when this species dominates in the early universe and radiation domination occurs, since this does not yield any effect on what we concern ourselves with here, i.e. gravitational particle production. We hereafter generalize the particle production in a radiation-dominated CP T -symmetric universe discussed in Ref. [3,4] to the non-standard cosmologies just described above. Ref. [3,4] considered a heavy sterile neutrino essentially decoupled from the thermal bath of the early universe and constructed the corresponding quantum field coupled to the expanding Friedmann-Robertson-Walker (FRW) metric that is asymptotically flat at early (τ → −∞) and late (τ → +∞) times. The Lagrangian for the spinor field in the expanding background is given by the usual with equation of motion Following [3], we posit that the right handed neutrino mass term originates from the vacuum expectation value of a field Φ whose equation of motion at early times, prior to the electroweak phase transition, is simply We call τ 1 the conformal time at which the universe switches from being dominated by a non-standard energy density and thus: and hence, the effective mass is: To find the constants in the equation above, first ask that the asymptotic value in the large τ limit give the "right" radiation domination answer, i.e. µ(τ ) = y(2ρ 1 ) 1/2 τ = γτ , with ρ 1 the value at the "beginning" of radiation domination, i.e. for us at conformal time τ 1 , and y = M/μ: this sets C 2 = 1; second, demand continuity of both Φ(τ ) (no constraints) and of its derivative, which imposes C 1 τ N −1 Given the symmetries associated with the FRW background, different observers disagree on what is a positive-or a negative-frequency mode: each observer will have their own corresponding mode functions, operators, and vacuum states. Ref. [3,4] label these the ingoing '−' and outgoing '+' states, for observers in the asymptotic, conformally invariant past, and in the non-conformally invariant late-time universe. The fields expansion for these modes are given by In the equation above, the operators a ± (p, h) and b ± (p, h) act on their respective vacuum states of the asymptotic observers, which are defined by In order to relate the field expansions and operators of the two asymptotic observers, one employs the Bogoliubov transformation. The key result of this transformation is that mode functions and operators for early-time observers can be expressed as a linear combination of the late-time observer's mode functions and operators, which form a complete set of states, as where α(p, h) and β(p, h) are the so-called Bogoliubov coefficients. If we apply the number operator of the ingoing observer N − = a † − a − to the outgoing observers vacuum state, the result is Integrating over the phase space gives the mean particle number density at late times, one gets the desired result The equation above implies that particle number is observer-dependent, which is a well-known fact in the context of the theory of quantum fields in expanding backgrounds. However, Ref. [3,4] where this quantity gives the mean particle density given the assumed CP T invariant vacuum state of the early universe. The Bogoliubov coefficient β(p) in the inverse sine is calculated considering the ingoing and outgoing vacuum states.
In the context of the generalized conformally invariant cosmologies we consider here, the Bogoliubov coefficient for the ingoing and outgoing states is corrected by an effective mass that now reads µ(τ ) = γτ N /τ N −1 1 . Since the fields mode functions form an orthonormal basis, one can follow the normalization condition the Bogoliubov coefficient β(p) in (19) can be then calculated using the inner product of ingoing and outgoing mode functions, The next task is to solve the equation of motion (7) and find the asymptotic form. For a free Weyl invariant scalar field defined on an expanding FRW metric with an effective mass µ the Lagrangian is given by with the equation of motion Taking the solutions to momentum space: ϕ(x) = u(p, τ )e ipx turns the equation of motion into that of a harmonic oscillator equation with frequency ω(τ ), To see how the mode solutions vary with conformal time τ , we numerically solved the equation for various N values, and compared with an approximate solution from the WKB formalism [17]. For the ingoing and outgoing observers the WKB solutions are given by where τ 0 is arbitrary and can be fixed for convenience. The integral can be carried out analytically, and the result is given by This expression simplifies considerably in the large τ limit. Here, we assume that particle production occurs at conformal times smaller than the time τ 1 at which the universe becomes radiation dominated. Thus, the large τ limit here means that τ 1 yT 1 1, i.e. τ 1 T 1 1 y = M P M , with M the right-handed neutrino mass. Using we have where H(a) H 0 √ Ω w a −(N +1)/N at a < a 1 by assumption, where Ω w is the relative energy density of the exotic component, We thus find Finally, we find that the condition does not depend on a 1 , and is always satisfied since, putting in the numbers (with H 0 ∼ 10 −42 GeV and Ω r ∼ 5 × 10 −5 ), we have implying M 10 −5 eV. In order to be the dark matter, the RHN must be heavier than at least a few keV, so the large-τ condition is always satisfied.
In the large τ limit, the first term in Eq. (26) reads Upon integration, one gets for that first term: For the hypergeometric function in Eq. (26), in the same limit, one can make the replacement The approximate asymptotic behavior of u ± (p, h, τ ) and υ ± (p, h, τ ) in the large τ limit is then found to be where, in the equation above, For convenience, above we indicate with γ = γ τ N −1 1 . Plugging this expression into the expressions above, we get the Bogoliubov coefficient β(p) as where we introduced the quantity ξ ≡ τ 0 /τ 1 . Notice that ξ depends on the (in principle arbitrary, as long as preceding Big Bang Nucleosynthesis) conformal time at which the universe transitions to being radiation-dominated, and on the parameter τ 0 , which is chosen so that Eq. (29) is valid (see [3] for details on this). Plugging this into (19) gives the generalized number density In the equation above, I is a dimensionless constant defined as Our result recovers the number density found in [3] for the radiation dominated background, N = 1. To extract a prediction for the dark matter (right handed neutrino) mass out of this number density, we follow the same procedure as in [3] (See their section 5.1.2). With n dm , we define the dark matter yield with s the entropy density in the universe, which we assume to be conserved at early times and in particular at the transition to radiation domination. After the reference conformal time τ 1 , we have the usual radiation-domination relations From these relations, solving for temperature, one finds Combining all expressions for n dm , γ and s, we can express Y dm in terms of the particle mass M dm where g * = 106.75 is the number of the degrees of freedom in the standard model, andμ = 5.19 × 10 8 GeV is the same right-handed neutrino reference mass scale chosen by [3]. Notice that the formula, remarkably, does not depend on the conformal time τ 1 at which radiation domination sets in (there could have in principle been a dependence both in the dark matter number density and via the energy density ρ 1 ).
We can now match the mass of the dark matter particle and the number density to the present day dark matter energy density ρ  Fig. 1 shows, as a function of the parameter ξ, the mass of the right-handed neutrinoν giving the right dark matter amount for a given value of the parameter N .
This setup does not imply values for parameters such as the Majorana and Dirac masses of the neutrino sectors, and thus the lifetime of the dark matter particle is not a prediction of the model. We describe ways to detect sterile neutrino decay products in this scenario in the following section, in a model-independent way.

III. STERILE NEUTRINO DECAY WIDTH
For simplicity, we assume the existence of a single, Weyl, sterile neutrinoν. We take the sterile neutrino/SM interactions to stem from a Yukawa interaction between the RH neutrino and the SM lepton doublet. In this scenario, the Lagrangian density is given by: In the expression above,m is the Majorana mass of the RH neutrino, y ν is a Yukawa vector with three components and φ, andL are the Higgs and lepton doublets (note that all the above fermion fields are written in two-component notation). The doublets are given by: with G + and G 0 being the charged and neutral Goldstone bosons, h being the standard model Higgs, v ∼ 246 GeV being the Higgs vacuum expectation value,ν i being the left-handed neutrinos andˆ i being the left-handed charged leptons. Expanding out the neutrino Lagrangian and gathering all the mass terms, we find the following: where the neutrino vectorν = ν 1ν2ν3ν T and the (4 × 4) neutrino mass matrix is given by: Diagonalizing the neutrino mass matrix is done through Takagi diagonalization with a unitary matrix Ω, satisfying Ω T M ν Ω = diag(m νe , m νµ , m ντ , mν). While it is possible to analytically compute Ω for general y ν , we choose to simplify our model by assuming only a single non-zero Yukawa coupling. That is, we take y k ν ≡ y and y k =0 ν = 0 with k = 1, 2 or 3. in this case, the unitary Takagi matrix is given by an orthogonal O(2) matrix times a diagonal unitary matrix: The angle satisfies with the light and heavy neutrino masses m ν and mν are given by: The mass eigenstates for the LH and RH neutrinos, ν k andν are related to the interaction eigenstates via Ω ν kν T , or more explicitly:ν k = −i cos θν k + sin θν (47) Using these transformations, we find that the Lagrangian density containing the interactions between the RH neutrino and the standard model relavent to this study is given by: where the · · · are terms not containing interactions between the RH neutrino and the SM or they contain irrelevant vertices (for example, interactions with the Goldstones.) The partial decay widths forν → W + + , W − + † , Z + ν and h + ν are given by: Taking mν m H , we find: The total decay with of the RH neutrino is roughly (for θ 1 mν m H ): and the corresponding lifetime τ = 6.6 × 10 −28 s θ 2 1000 GeV mν meaning that where we indicate with τ U the age of the universe.

IV. STERILE NEUTRINO THERMALIZATION
A necessary condition for the validity of the dark matter mass predictions discussed above is that the right handed neutrinos not thermalize: should they reach chemical equilibrium with the Standard Model thermal bath, the corresponding thermal relic density would be unrelated to the abundance from gravitational production. We must thus ensure that the parameter space under consideration be consistent with the absence of thermal equilibrium for the neutrino.
At T > T 1 , where T 1 is the temperature at which the energy density of radiation equals the energy density that dominates in the early universe during gravitational dark matter production, Matching ρ rad (T 1 ) with the expression above gives Right handed neutrinos can be in thermal equilibrium only when relativistic, which allows us to express the corresponding number density as nν ∼ T 3 (T > mν). At large-enough temperatures, larger than both mν and the electroweak scale, the cross section responsible for neutrino thermalization is approximately σν −SM ∼ θ 2 /T 2 , from decay and inverse decay into e.g. a gauge boson and a lepton.
In the radiation domination phase (T < T 1 ), when the Hubble rate is approximately H ∼ T 2 /M P , thermalization happens at a temperature such that Since by assumption T < mν, a sufficient condition to prevent thermalization is that mν < M P θ 2 . From the discussion above, Eq. (58), for TeV and heavier right handed neutrinos, M P θ 2 10 −28 GeV, so thermalization can never occur during radiation domination.
In the non-standard phase, T > T 1 , the Hubble rate is approximately thus to ensure that thermal equilibrium never be attained we simply need to require that Below the Planck scale, nνσν −SM < M P θ 2 , implying that equilibrium can only be attained for T 1 ∼ θM P which is much smaller than the scale of Big Bang Nucleosynthesis T BBN ∼ 1 MeV, at which the universe must be radiation dominated (in other words, T 1 > T BBN . Thus, thermal equilibrium is never attained at any point below the Planck scale.

V. CONSTRAINTS ON STERILE NEUTRINO LIFETIME
Tight constraints on the sterile neutrino lifetime can be computed using experiments such as the Fermi Large Area Telescope (LAT) [18], IceCube [19], and HAWC [20].
The gamma-ray and neutrino differential fluxes detectable on Earth from the decay of a sterile neutrino are given by: where τ is the sterile neutrino life-time, dN γ,ν / dE γ,ν is the gamma-ray/neutrino spectrum per decay and J is the so-called J-factor, the integral along the line of sight, and averaged on the telescope's angular aperture, of the dark matter density, which depends on the target being observed. Unlike for annihilation, the calculation of the J factor for decay is prone to significantly less severe uncertainties from the slope of the inner density of the dark matter profile. Here, we use the same J factors as in [21,22]. Computing the gamma-ray/neutrino spectrum for a sterile neutrino mass mν ∼ 10 8 GeV is a non-trivial task. Since mν m W , m Z , m H , the electro-weak bosons emitted from the sterile neutrino decays are highly-relativistic and are prone to undergo large numbers of splittings into more electro-weak bosons. For this reason, we use a newly developed software package, HDMSpectra, specifically developed for computing spectra from the decay of very heavy dark matter (m DM m H ) [23]. Fig. 2 shows the gamma-ray and neutrino spectra from the decay of a sterile neutrino with a mass of mν ∼ 4.8 × 10 8 GeV.
To compute the constraints from HAWC data, we use the differential fluxes provided in Ref. [20]. In that study, the HAWC collaboration measured the gamma-ray fluxed from the Andromeda galaxy, an ideal target to look for dark matter decay debris. They measured fluxes of photons with energies ranging from a TeV to 100 TeV and were able to place constraints on E 2 dN γ / dE γ 10 −12 − 10 −11 TeV cm −2 s −1 . Using these results, we find a constraint of the sterile neutrino lifetime of roughly τ 7 × 10 24 s (see Fig. (3)).
To calculate constraints for IceCube, we use the Galactic center neutrino upper limit fluxes shown in Fig. 2 of Ref. [19]; similarly, for the constraints from the Fermi Large Area Telescope (LAT) [18], we use the sensitivity discussed in Ref. [18].
We show the upper limits to the lifetime of a decaying right-handed neutrino in fig. 3. By far the best limits arise from the neutrino telescope IceCube, which in the mass range of typical interest, between 10 3 and 10 8 GeV, achieves a sensitivity in excess of 10 29 s, two to three orders of magnitude better than gamma-ray telescopes. This lifetime sensitivity translates into probing mixing angles as small as θ ∼ 10 −59 .  : Upper limits on the right handed neutrino lifetime as a function of the particle mass, from IceCube [19] (blue line), HAWC [20] (red line), and Fermi LAT [18] (yellow line).

VI. SUMMARY, DISCUSSION, AND CONCLUSIONS
We generalized gravitational particle production in the early universe in CP T symmetric cosmologies beyond the case of radiation domination discussed in Ref. [3,4]; our key result is that in non-standard, CP T symmetric universes the mass of the dark matter candidate, a "sterile" right-handed neutrino, depends on the transition from the nonstandard cosmological epoch to radiation, and on the dependence of the energy density on conformal time at early times. As a result, unlike the case of radiation domination, the dark matter mass in largely unconstrained and could range between the electroweak and the Planck scales. Note that we did not consider cases where the expansion history of the universe, while respecting CP T invariance, changes at multiple points at early times: particle production in that case would happen during one particular phase, which would fall within the cases considered here.
We relaxed the assumption that the dark matter be protected from decay by a Z 2 symmetry, and discussed the possible detection of massive, right-handed neutrinos. We verified that the neutrino is never in thermal equilibrium in the early universe. We found that the optimal way to search for the neutrino's decay products is to employ highenergy neutrino telescopes such as IceCube, which at present rules out neutrinos with lifetimes in excess of 10 29 s up to masses as large as 10 8 GeV.
Future tests of this mechanism for the generation of the cosmological dark matter will depend upon observational indications for the lack of an inflationary period, as well as in favor of a CP T symmetric early phase; the detection of the debris from the decay of a superheavy species compatible with that species being a right-handed neutrino would also be a smoking-gun signature for this scenario.