Comprehensive constraints on heavy sterile neutrinos from core-collapse supernovae

Sterile neutrinos with masses up to $\mathcal{O} (100)$ MeV can be copiously produced in a supernova (SN) core, through the mixing with active neutrinos. In this regard the SN 1987A detection of neutrino events has been used to put constraints on active-sterile neutrino mixing, exploiting the well-known SN cooling argument. We refine the calculation of this limit including neutral current interactions with nucleons, that constitute the dominant channel for sterile neutrino production. We also include, for the first time, the charged current interactions between sterile neutrinos and muons, relevant for the production of sterile neutrinos mixed with muon neutrinos in the SN core. Using the recent modified luminosity criterion, we extend the bounds to the case where sterile states are trapped in the stellar core. Additionally, we study the decays of heavy sterile neutrinos, affecting the SN explosion energy and possibly producing a gamma-ray signal. We also illustrate the complementarity of our new bounds with cosmological bounds and laboratory searches.


I. INTRODUCTION
In the last decade, new theoretical ideas to address dark matter and other fundamental questions predict a dark sector composed of feebly interacting particles (FIPs) with sub-GeV masses and very feeble interactions with Standard Model (SM) particles [1][2][3][4].The most common approach to describe the interaction of the dark sector with the SM is through some portal.At this regard, the minimal portals mixing new dark sector states with gauge-invariant combinations of SM fields are: vector (dark photons), scalar (dark Higgs), fermion (heavy neutral leptons) and pseudo-scalar (axions) [1].These portals are subject of intense experimental investigations with interesting plans for the next years [1][2][3][4].
In this paper, we will focus on another class of FIPs, namely heavy neutral leptons and in particular a heavy sterile neutrino, ν 4 , mostly a flavour-sterile one (ν s ) with a (generally small) mixing with active neutrinos ν α , α = e, µ, τ .These states have been often introduced to explain the origin of neutrino masses [23][24][25][26].We remark that although the sterile neutrino scale considered here is not heavy for particle physics standards, it is so if compared to the current bounds on the mass scale in the active neutrino sector, i.e. m ν ≲ 1 eV.Therefore, we will use the adjective heavy in this sense.
This constraint has been recently re-evaluated in the free-streaming regime in Ref. [38], considering weaklymixed sterile neutrinos that escape the SN without interacting with stellar matter.However, recent develop-ments in SN simulations and new proposals to improve FIP constraints from SNe suggest that the heavy-sterile neutrino limits can be significantly strengthened.It is worth noting that recent works, as Ref. [38], have considered the scattering of active neutrinos as the dominant channel for sterile neutrino production, neglecting the neutral current interactions with nucleons.The latter channel was expected to be suppressed, due to the Fermiblocking associated with nucleon degeneracy in the SN core.Nevertheless, in the few cases where nucleon scattering was considered, as in the seminal papers [42,43], the corresponding bounds were stronger than the ones obtained, e.g., in Ref. [38].However, since the treatment of these processes in Refs.[42,43] was cursory, it seems to us important to revisit and clarify this issue.Additionally, from recent SN simulations [50,51] it emerges that a population of muons is present in the core and neutrinos interact with them through charged current interactions.These interactions are especially relevant in enhancing the production of sterile neutrinos mixed with muon neutrinos, allowing for an improvement of the previous bounds.
Furthermore, it is possible to constrain the ν 4 parameter space by studying the energy deposited inside a SN via the electromagnetic decays of sterile neutrinos.This is relevant for massive sterile neutrinos, where various decay channels are possible.For example, the decay ν 4 → π 0 ν µ would deposit at least 135 MeV of energy inside the SN [32,40].At this regard, it has been recently shown in Ref. [15] that in order not to exceed the explosion energy observed in low-energy SNe, strong constraints can be placed on energy deposition induced by FIP decays.This argument has been applied to the heavy sterile neutrino case in Ref. [52].Finally, the flux of daughter particles produced outside the SN, especially e + and γ, may lead to strong bounds (see Ref. [53] for a seminal study on ν 4 ) similarly to the ones recently discussed in Refs.[17,[54][55][56][57][58] for the case of heavy axion-like particles.
Given these motivations, we devote this work to strengthen the existing bounds on heavy sterile neutrinos from SNe exploring different aspects: • including neutral current interactions of ν 4 with nucleons; • including charged current interactions of ν 4 with muons; • characterizing the trapping regime of ν 4 verified at large mixing angles, adopting the the so-called "modified luminosity criterion", (see Refs. [13,19,59]).This recently proposed recipe allows one to extend the SN energy-loss bounds also to the regime where ν 4 are strongly interacting with matter; • considering (non)radiative decays of heavy neutrinos, we strengthen the cooling including the constraint from excessive energy deposition, following the method proposed in Ref. [15,52], and from an observable gamma-ray signal.
The plan for this paper is as follows.In Sec.II we recall the heavy neutrino production in SNe and summarize the relevant production and absorption processes.Then in Sec.III we discuss the different arguments presented in the literature to constrain FIPs from SNe and we apply them to the case of heavy ν 4 .In Sec.IV we combine all our bounds and compare them with the other laboratory and cosmological constraints in the same mass range.We conclude in Sec.V.In App.A we discuss the details of the evaluation of the charged and neutral current interactions involving ν 4 .

II. STERILE NEUTRINO PRODUCTION
We limit ourselves to heavy sterile neutrinos with masses 10 MeV ≲ m 4 ≲ 600 MeV [60,61], to avoid any possible resonant production which usually happens in the sub-MeV range [33,34,62].In this mass range since the mixing of a sterile neutrino with electron neutrino is very constrained (see, e.g., Ref. [3]), we assume that the sterile neutrino is mixed dominantly with one active neutrino ν α , with α = µ, τ , such as where ν ℓ and ν 4 are a light and the heavy mass eigenstate, respectively, U is the unitary mixing matrix, linking mass and flavour states, and the most interesting parameter space corresponds to |U α4 | 2 ≪ 1, i.e. ν ℓ is mostly active and ν 4 is mostly sterile.
In the SN core, sterile neutrinos are produced via the processes listed in Tab.I. We characterize these processes closely following Ref.[38].We have neglected the bremsstrahlung process N N → N N ν ν since it is always sub-leading in the interesting parameter space.Indeed, as a production channel, the computed luminosity according to the rate of [9] is inferior to the one associated to the other processes in Tab.I.As an absorption channel, it is suppressed compared to N ν → N ν for obvious phase-space reasons.
The production rate of sterile neutrinos per unit volume and energy can be written as where E 4 and p 4 are energy and momentum of the sterile neutrino, f i is the distribution function of i-th particle involved in the process, |M| 2 12↔34 is the sum of the squared amplitudes for collisional processes 1 + 2 ↔ 3 + 4 relevant for the sterile production/absorption, reported in Tab.I. Given recent SN simulations including muons [50,51], V V are reported in App. A. In the µ N ↔ N ′ ν4 matrix element we neglected the terms of higher order in the nucleon momenta, thus working at leading order in a non-relativistic approximation.Here, N and N ′ represent the different nucleons involved in the interaction in the initial and final states.Finally, for all the processes we have considered their corresponding charged conjugate for the ν4 production (the last two processes are irrelevant because of the absence of µ + in the SN core) and absorption.
here we consider for the first time reactions involving muons, also reported in Tab.I.Moreover, we also include the neutral current interaction between neutrinos and nuclei, which results to be one of the main channels for the production and absorption of the heavy state.Despite the fact that its possible relevance was already pointed out in Ref. [64], in most literature it has been neglected, similarly to the process ν e − → ν e − , due to the large assumed fermion degeneracy and Pauli blocking effect.The abundance and degeneracy of nucleons in the SN core can be assessed by considering that the nuclear medium is described in a relativistic mean-field (RMF) picture [65], according to which the nucleon distribution function f N is given by [65] f N (p) = 1 where m * N = m N + Σ S is the effective nucleon mass, with Σ S the so-called nuclear scalar self-energy, and and µ * N the effective or kinetic chemical potential, defined as [65] µ where µ N is the nucleon chemical potential including the nucleon rest mass and Σ V is the RMF vector selfenergy.Thus, nucleons have Fermi-Dirac distribution functions equivalent to a non-interacting system with effective chemical potentials µ * N and particle masses m * N and their degeneracy can be estimated by introducing the degeneracy parameter η N defined as On the other hand, leptons in the SN simulations are described by the usual Fermi-Dirac distributions with m l their bare mass and µ l their chemical potential, leading to the degeneracy parameter We mention here that electrons in the plasma acquire an effective mass that for typical SN conditions (µ e ≫ T ≫ m e ) can be written as m * e 2 = e 2 (µ 2 e + π 2 T 2 )/8π 2 [66].This expression leads to m * e ≲ O(10) MeV ≪ µ e in the SN core.Thus, using m e or m * e marginally affects the evaluation of η e in Eq. ( 7) [67] and, consistently with our benchmark SN model described in the following, we neglect m * e in our analysis.Particles i in the plasma are non-degenerate if η i < 0, while they are fully degenerate for η i ≫ 1 and only partially degenerate for intermediate values of η i [9].
We compute the sterile neutrino production using as a benchmark an 18 M ⊙ progenitor mass (roughly consistent with Sanduleak-69 202, the progenitor of the SN 1987A) obtained using a 1D spherically symmetric and general relativistic hydrodynamics model, based on the AGILE BOLTZTRAN code [68,69], including muons.While we expect that these simulations capture the basic physics of the phenomenon, differences of a factor of a few can be associated to the implementation scheme of the neutrino microphysics, general relativistic effects, multi-dimensionality, etc.We think that this constitutes the dominating systematic error in the derived bounds.
We show in Fig. 1 the thermodynamical conditions for our benchmark SN model in the inner core (r ≲ 20 km) at the post-bounce time t pb = 1 s.The upper panel shows the temperature T (solid black line), with a peak T ∼ 40 MeV at r ∼ 10 km, and the matter density ρ (dashed black), with a maximum ρ ∼ 4 × 10 14 g cm −3 at the center and decreasing at larger radii.The central panel shows the fermion degeneracy parameters η α for nucleons α = n, p and leptons α = e, µ.In the very inner core (r ≲ 5 km), neutrons (solid black line) are degenerate (η n ≈ 5) and protons (dashed black line) are partially degenerate (η p ≈ 2).For larger radii (r ≳ 10 km), the nucleon degeneracy decreases, implying non degenerate protons (η p < 0) and only partially degenerate neutrons (η n ≲ 1).On the other hand, throughout the SN core, electrons (dotted black line) are highly degenerate (η e ≳ 5) and muons (dot-dashed black line) are nondegenerate (η µ < 0).This implies that the ν e − → νe − process is suppressed by the electron degeneracy, while neutral current interactions with nucleons cannot be neglected, at least in the outer layers of the core.Moreover, we checked that the conclusions concerning the η e parameter are unchanged (with a discrepancy lower than 5%) even considering the effective electron mass due to QED at finite temperature and density, yielding m * e ∼ O(10) MeV at r ≲ O( 10) km [66,67].Finally, in the lower panel we present the electron Y e (solid line) and muon Y µ (dashed line) abundance with respect to the nucleon one, with Y α = n α /n B , where n α is the density per unit volume for the particle α = e, µ and n B is the baryon number density.We realize that the muon abundance around the peak of the temperature can be O(10%) of the nucleon one.Therefore, for definiteness, we evaluate the sterile neutrino production by taking into account also processes involving muons.
We compute the production rate for sterile neutrinos by reducing the nine-dimensional integral in Eq. ( 2) to a three-dimensional one following the procedure in Ref. [70].As an example, in App.A, we show how it is possible to write the interaction matrix elements for the charged current process process µ N ↔ N ν 4 and the neutral current interaction ν α N ↔ N ν 4 using the formalism in Ref. [70].The same procedure can be applied to evaluate the interaction matrix elements for the other processes we consider.

III. SN CONSTRAINTS ON HEAVY STERILE NEUTRINOS A. Cooling bound
From the observation of the νe neutrino burst from SN 1987A [44][45][46][47], it is possible to infer the temporal evolution of the neutrino lightcurve.Despite the sparseness of the data, the duration of the neutrino burst ex- Lower panel : Radial profiles of the electron fraction Ye (solid line) and muon fraction Yµ (dashed line).All panels refer to the post-bounce time t pb = 1 s.tending over 10 s is in agreement with the expectations from the SN cooling via neutrinos [9].Therefore, from the SN 1987A neutrino data there is no evidence of a dominant emission of FIPs, that would have significantly shortened the duration of the neutrino burst [8,9]. 1n order to avoid a significant shortening of the observed neutrino burst due to FIP emission, one should require that the luminosity of the exotic particles should be less then the one carried by neutrinos.Namely, for our fiducial model at t pb = 1 s one has [8,9,59] Our goal is to use this constraint to exclude values of the ν 4 mixing with muon neutrinos (|U µ4 | 2 ) and tau neutrinos (|U τ 4 | 2 ).We do not consider the case of mixing with the electron flavour, since in this case the parameter space is overconstrained.We adopt the "modified luminosity criterion" [13,19,59] to smoothly interpolate between the regimes in which sterile neutrinos are so weakly interacting that they freely escape from the SN (i.e.weak-mixing regime with |U α4 | 2 ≲ 10 −5 , see Fig. 2), also known as free-streaming regime, and a regime of stronger interactions with matter (i.e.strongmixing regime with |U α4 | 2 ≫ 10 −5 , see Fig. 2), when they are trapped in analogy with active neutrinos.In this formalism, the ν 4 luminosity is [13,19,59,72] where α is the lapse factor to account for the gravitational redshift and the exponential suppression e −τ takes into account the possibility of ν 4 absorption inside the SN.In particular, ⟨ e −τ ⟩ is a directional average of the absorption factor [59,[72][73][74]] where λ is the sterile neutrino mean-free path (mfp), s µ is the ν 4 redshifted energy, µ = cos β and β is the angle between the outward radial direction and a given ray of propagation along which s is integrated.We emphasize that, lacking self-consistent SN simulations including the feedback due to the emission of sterile neutrinos, as the rest of the literature (implicitly) does, we also resort to an extrapolation whenever the extra neutrino luminosity is comparable with or larger than the luminosity of the species ν α mixed with ν 4 (see dashed lines in Fig. 2).The results at values much larger than the active neutrino luminosity are however only nominal, and not essential in obtaining the bound.Yet, it is conceivable that this limitation may introduce a factor of a few uncertainty in the limits from the cooling argument.
Sterile neutrinos with sufficiently strong interactions with ordinary matter are trapped in the SN via interactions with active neutrinos, electrons, positrons and neutrons.The considered processes are listed in Tab.I.Moreover, depending on their mass, sterile neutrinos may decay in different particles after their production, through the processes shown in Tab.II.In Fig. 3, we show the branching ratios for the relevant decay channels.When evaluating the ν 4 absorption mfp, we have to consider absorptions and decays separately.In the former case, the mfp is defined as where n is the density of targets and σ is total absorption cross-section.Since all absorption processes in Tab.I are 2 → 2 scatterings, it is possible to write the following TABLE II: Decay channels up to m4 ≲ 250 MeV, for a ν4 mixed with νµ where fπ = 135 MeV, gL = − 1 2 +sin 2 θW , gR = sin 2 θW , and the electron mass is neglected [38,75].The decay mode into two muons is neglected since it is characterized by a small branching ratio for m4 < 500 MeV, as shown in Fig. 3.The decay processes for the sterile neutrinos mixed with ντ via |Uτ4| 2 are the ones not involving a single muons in the final states, i.e. all but ν4 → νe e + µ − and ν4 → µ − π + .expression for the cross section with a suitable choice of Fermi-Dirac distribution functions f i for i = 1, 2, 3, and of the matrix element |M| 2 , taken from Tab. I.In this context, the mfp can be explicitly evaluated by employing the procedure discussed in Ref. [70].
Regarding the sterile neutrino decays, the mfp is defined as where Γ tot is the sum of the decay widths Γ of all the relevant decay processes (see Tab. II), γ = (1−β 2 ) −1/2 is the Lorentz factor and the ν 4 velocity is β = p/E 4 .Note that Dirac neutrinos are implicitly assumed throughout our paper; for Majorana states, the rates for exclusive processes are the same, but L-conjugated processes, e.g. in decays, are also allowed, thus doubling the inclusive rates.By combining Eq. ( 11) and Eq. ( 13), we obtain the total mfp as which is used to evaluate the sterile neutrino luminosity in Eq. ( 9) and impose the constraint in Eq. ( 8).Here it is important to mention that for an unstable particle, the absorption factor in Eq. ( 10) would always be zero if the integration limit in the exponential were infinity (since unstable particles can decay in vacuum).However it is possible to fix the upper integration limit in Eq. ( 10) to a far radius R far = 100 km [19,76,77], much larger than the protoneutron star radius R PNS ≈ 10 km where sterile neutrinos are produced.This choice for R far is an order-of-magnitude estimation for the neutrino gain radius [77].In this way, we assume that only sterile neutrinos reaching the gain radius contribute to the SN cooling.Uncertainties related to the choice of R far will not affect the conclusions of this work, since the cooling bound is always subdominant with respect to other constraints, as discussed in the following.
The obtained bound is shown in Fig. 4. The contour area delimited by the solid line refers to the mixing with ν µ , while the dashed line represents the bound for ν τ mixing.This criterion excludes a region between O(10 −9 ) < |U α4 | 2 < O(10 −2 ) for m 4 ∼ O(10) MeV, probing masses up to ∼ 400 MeV for |U α4 | 2 ∼ O(10 −7 ).Our bound agrees at the order of magnitude level with the bound estimated in the seminal papers [42,43].In the same figure, the dot-dashed black line shows the lower bound obtained neglecting the ν N → N ν 4 interaction.Notice that this constraint is in agreement with the one SN cooling ) for α = µ (solid line) and α = τ (dashed line).The dot-dashed black line is the bound for |Uτ4| 2 which would be obtained neglecting the νN → N ν4 process as in Ref. [38].
obtained in Ref. [38] under the same assumption 2 .Thus, we confirm the crucial importance of the inclusion of neutral interaction processes with nucleons in obtaining the lower exclusion bound on |U α4 | 2 .Let us also highlight that the bound for α = µ is a factor ∼ 2 stronger than α = τ because of the larger sterile neutrino luminosity caused by the presence of muons.Henceforth, bounds in the literature ignoring this effect tend to be too conservative.Note that in principle, following a calculation similar to the one in [38], a related bound might be obtained by the non-observation of a high-energy neutrino flux, from sterile neutrino decay, in the experiments that detected the neutrino signal from the SN 1987A.Our estimates suggest that while comparable or slightly better than the cooling bound, this would not be competitive with other constraints discussed below, and will not be considered further in this article.

B. SN explosion energy bound
As we can see from Tab. II, all the sterile neutrino decay channels except for ν 4 → ν µ ν α να (with α = e, µ, τ ) produce photons, leptons or pions.If sterile neutrinos have a decay length between the core radius of about 10 km and the progenitor star radius of about 10 13 cm, they decay inside the SN envelope, depositing at least part of their energy inside the star.This phenomenon allows us to use SNe as efficient calorimeters.As proposed 2 As a corollary, the magnitude of the signatures associated to benchmarks discussed in [38] remain roughly the same if the benchmark mixings adopted there are scaled-down by the same amount as the tightening of the bounds in presence of neutral current interactions with nuclei.
in Refs.[15,21], there is an upper limit on the amount of energy that can be deposited inside a SN by FIP decays without producing too energetic explosions that would be incompatible with observations of low-energy SNe.This constraint requires that E e.m.FIP ≲ 10 50 erg , where E e.m.FIP is the energy released in the electromagnetic sector by sterile neutrino decays.
In the decays, we assume that the daughter particles are emitted with an appropriate fraction of the energy in the center-of-mass frame, depending on the channel.Following Ref. [38], it is possible to write the deposited energy as E e.m.
where the index i runs over the decay processes under consideration, B i is the branching ratio of the i-th process, Ē and E are the daughter particle energies in the center-of-mass and in the laboratory frame, respectively, E 4 is the sterile neutrino energy, R * = 2.5 × 10 13 cm is the stellar radius and [15, 38] with d 2 N 4 /dE 4 dt accounting for the fraction of sterile neutrinos escaping from the core, with R p = 40 km.
In this way, we take into account only the energy carried out from the core and deposited by sterile neutrinos decaying in the SN envelope.We expect the bound in Eq. ( 15) to set a constraint on |U α4 | 2 two orders of magnitude more stringent than the SN cooling bound, for sufficiently high ν 4 masses.At lower masses, the longer lifetime and larger boost factors imply that decays are not efficient and this constraint is relaxed.Indeed, in Fig. 5 we see that for The bump in the constraint around 135 MeV reflects the opening of an extra decay channel for sterile neutrinos, ν 4 → ν α π 0 .Recently, the authors of Ref. [52] have obtained a SN explosion energy bound without considering the νN → ν 4 N interaction.Similarly to the SN cooling case, neglecting the neutral current interactions with nucleons leads to a constraint two orders of magnitude weaker than the one obtained in Fig. 5.

C. 511 keV bound
Sterile neutrinos escaping the SN envelope and decaying in the interstellar medium give rise to a diverse phe-Expl.energy nomenology, depending on the considered decay products.Here, we focus on the positrons produced by a portion of the ν 4 -decay channels.
As extensively discussed in Refs.[55,56,78,79], this exotic injection of positrons in the Galaxy would originate a distinctive soft gamma-ray signal.Precisely, positrons emitted by sterile neutrino decays are trapped in the Galaxy by its magnetic field.While traveling on scales smaller than O(1) kpc from the decay point, positrons lose energy by Bhabha scattering on the galactic electron population.This thermalization process lasts between 10 3 and 10 6 yrs, depending on the electron density.This long time-scale explains why the positron injection, caused by SNe during the history of the Galaxy, can be assumed continuous.Once positrons are almost at rest, ∼ 25% of them form a parapositronium bound state with an electron, before decaying in two back-toback photons, each one with an energy of 511 keV, determined by the electron rest mass [80].
A Galactic 511 keV line, at least partially explained by standard positron emission mechanisms, is prominently observed from the direction of the Galactic bulge [81].The contribution to this signal induced by sterile neutrinos can be calculated as where dΩ = dl db cos b, with −π ≤ l ≤ π being the longitude and −π/2 ≤ b ≤ π/2 being the latitude in the Galactic coordinate system (s, b, l), with s distance from the SN to the Sun.Moreover, k ps = 1/4 accounts for the fraction of positrons annihilating through parapositronium.According to Ref. [82], we fix Γ cc = 2.30 SNe/century as the Galactic SN rate.Finally, n cc is the SN volume distribution [83] in the Galactocentric coordinate system (r, z, l), with r the galactocentric radius and z the height above the Galactic plane, connected with the Galactic coordinate system through the relations Here, we set the solar distance from the Galactic center to d ⊙ = 8.5 kpc.Requiring the photon flux in Eq. ( 18) to be smaller than the observed signal in the range l ∈ [28.(20) This is the most conservative limit obtained by the comprehensive analyses of Refs.[78,79], taking into account different SN distribution models and diffusive smearing effects.This upper bound on N pos corresponds also to the constraint placed by XMM-Newton observations of the Galactic X-ray background [84].Indeed, an excess of electron/positron injection in the Galaxy would source a diffuse X-ray signal via inverse Compton scattering on the stellar background light.
In order to apply the constraint in Eq. ( 20), we calculate the number of injected positrons as with the average number of positrons produced in a sterile neutrino decay.Moreover, following Ref.[20] we fix for the envelope radii of Type II and Ib/c SNe, while according to Ref. [82], we take as average fractions of SNe of Type II and Ib/c In Fig. 6 we show the calculated N pos as a function of the mixing angle.At low mixing, sterile neutrinos are not efficiently produced and, therefore, the number of positrons produced in the decay is smaller than the limiting value, represented by the dotted line.Then, as the sterile neutrino production increases, given that almost the totality of neutrinos decay inside the Galaxy, the injected positrons can be a sizable number.As it can be seen in Fig. 6, light sterile neutrinos with m 4 = 20 MeV (solid lines) can produce up to ∼ 10 58 positrons per SN.A smaller number is obtained by more massive neutrinos (dashed lines), since their production is Boltzmann suppressed.For small values of the mixing (e.g., |U α 4 | 2 ≲ 10 −6 for m 4 = 20 MeV), we notice a relatively small difference between sterile neutrinos mixed with muon neutrinos (black lines) or tau neutrinos (red lines), due to the larger production of the former ones induced by charged current interactions of muons with nucleons.For larger values of the mixing (e.g., |U α 4 | 2 ≳ 10 −6 for m 4 = 20 MeV), the number of positrons is exponentially suppressed and the different production and absorption processes lead to an even smaller difference in the positron production.The bound obtained with this approach is expected to exclude relatively light ν 4 , with masses above a few tens of MeV, and it can be extended to small couplings because it is a cumulative diffuse flux.Indeed, we can see from Fig. 7 that the obtained lower bound is |U α4 | 2 ∼ O(10 −11 ) for m 4 < O(100) MeV.

D. SN 1987A gamma-ray bound
As discussed in the previous subsection, sterile neutrinos decaying after escaping the SN envelope lead to peculiar signatures.One of the most powerful constraints is given by the non-detection of a gamma-ray signal in coincidence with the neutrino burst of SN 1987A, as studied in the seminal work of Ref. [53].The Gamma-Ray Spectrometer of the Solar Maximum Mission places an upper limit of [85] on the photon flux at energies between 25 MeV and 100 MeV for 232.2 s after the first neutrino arrival.This upper limit translates into a constraint on |U α4 | 2 , since the radiative decay of massive ν 4 would give rise to a gamma-ray signal in coincidence with a SN explosion.From Tab.II, we notice that the only decays contributing to this signal are ν 4 → ν α γ and ν 4 → ν α π 0 , because of the successive decay π 0 → γγ.The spectrum of photons originated by ν 4 decay directly into photons is [53] dN where the average energy, in the center-of-mass frame, of the daughter particle j, from the decay of the parent which is larger or equal to when expressed in the laboratory frame.In addition, the fraction of sterile neutrinos decaying outside the SN envelope is Similarly to Eq. ( 26), we can evaluate the energy spectrum of neutral pions produced by sterile neutrinos as In a second step, the gamma-ray spectrum from the almost immediate pion decay is obtained as [53] dN In conclusion, the expected gamma-ray flux can be written as where We set stringent constraints on the sterile neutrino properties by integrating Eq. ( 32) over the observation time of 232.2 s and comparing the result with the limit in Eq. (25).The constraint obtained in this way becomes particularly relevant as soon as ν 4 is heavier than the pion, opening up the pion decay channel.Our results are reported in Fig. 8, showing that the lower bound strengthens from |U α4 | 2 ≳ O(10 −10 ) for m 4 < m π , m π being the pion mass, to |U α4 | 2 ≳ O(10 −12 ) for m 4 > m π due to the decay channel ν 4 → ν α π.As argued in [16] for axionlike particles, we mention that for sterile neutrinos with masses of a few 10 MeV the decay-product photons may create a fireball, making part of the "SN 1987A gamma-rays" bound not valid.However, the fireball would produce a gamma-ray flux with energy of a few MeV and the non-detection of such a signal in coincidence with the SN 1987A burst by Pioneer Venus Orbiter (PVO), constrains again this region, which in turn is already excluded by other bounds such as the 511 keV one.

E. Diffuse gamma-ray bound
The same phenomenology discussed above can be applied to evaluate the cumulative gamma-ray flux induced by SN ν 4 during the history of the Universe.This would constitute a diffuse, isotropic and constant gamma-ray flux at a few tens of MeV.
The gamma-ray spectrum for a single SN is calculated as in Eq. (32), redshifted in energy and integrated over the SN explosion rate as (see Ref. [86] for calculation details in the case of the SN diffuse neutrino spectrum and Ref. [87] for the axion case) where z is the redshift, R SN (z) is the SN explosion rate taken from [88], with a total normalization for the corecollapse rate with the cosmological parameters fixed at H 0 = 67.4kms −1 Mpc −1 , Ω M = 0.315, Ω Λ = 0.685 [89].The flux in Eq. ( 34) is imposed to be smaller than extracted from measurements of Fermi-LAT of the diffuse gamma-ray background [54].The advantage of this constraint, reported in Fig. 9, is that it extends to smaller masses, where the decay rate is less efficient, excluding |U µ4 | 2 ≳ 2 × 10 −11 and |U τ 4 | 2 ≳ 3 × 10 −11 for m 4 ≲ 100 MeV.

IV. COMBINATION OF DIFFERENT BOUNDS
In Fig. 10 we combine all the bounds obtained in the previous Sections for ν 4 mixed with ν µ (upper panel) and ν τ (lower panel).

A. Laboratory bounds
Constraints on heavy sterile neutrinos decaying into leptons and pions are set by the long-baseline neutrino oscillation experiment T2K [97].A beam of 30 GeV protons produced a large amount of kaons in their scattering FIG.10: Overview of the bounds from SNe (green region), cosmology [63,90,91] (blue region) and experiments [92] (gray region) for sterile neutrinos mixed with muon neutrinos (upper panel) and tau neutrinos (lower panel).The dashed lines represent the sensitivities of the future experiments DUNE [93,94] (black), SHiP [95] (red) and MATHUSLA [96] (brown).The hatched area below the dotted blue line represents the region of the parameter space in which sterile neutrinos are produced in the early universe out of equilibrium, whose viability is more model-dependent.
on a graphite target at J-PARC.Then, kaons might produce sterile neutrinos in their decay.A detector placed at a baseline of 280 m was used to reveal the decay of sterile neutrinos.The constraints obtained by T2K complement and improve the results of CHARM [98] and PS191 [99].Current experimental bounds are shown as the gray shaded area in Fig. 10.
The aforementioned searches will be improved by future experiments, whose projected sensitivities are shown as dashed lines in Fig. 10.In particular, current experimental constraints on sub-GeV sterile neutrinos considered in this work will be strengthened by DUNE [93,94], probing however regions not excluded by SN arguments only for masses m 4 ≳ 400 MeV, as shown by the dashedblack line in Fig. 10.Moreover, the future beam-dump experiment SHiP [3] is designed to probe exotic long-lived particles produced by a 400 GeV proton beam from the Super Proton Synchrotron at CERN, allowing the exploration of a much larger region of the parameter space for sterile neutrinos mixed with muon neutrinos [95], as represented by its sensitivity (dashed-red line) in the upper panel of Fig. 10.On the other hand, as shown by the dashed-brown line in the lower panel of Fig. 10, a currently unexplored region of the parameter space of sterile neutrinos mixed with tau neutrinos will be probed by MATHUSLA [96], another CERN experiment planned to study sterile neutrinos by searching for displaced vertex signatures near the LHC interactions points.

B. Cosmological bound
Constraints on heavy sterile neutrinos from cosmological observations emerge considering that their decay, after the active neutrino decoupling, generates extra neutrino radiation and entropy production in the Early Universe.Therefore, they alter the value of the effective number of neutrino species N eff , measured by the cosmic microwave background (CMB), and affect primordial nucleosynthesis (BBN), notably 4 He production, which is reflected in the Y p value.Using the latest measurements of the Planck collaboration [89,100], it is possible to obtain cosmological constraints, see [63,90,91].These arguments exclude up to |U α4 | 2 ∼ O(10 −1 ) for m 4 ≈ 10 MeV, as represented by the blue region in Fig. 10 labelled as "COSMO".For heavier sterile neutrinos, with mass m 4 > m π , the strongest impact on BBN is induced by the meson-driven p ↔ n conversion, which significantly increases the helium abundance and constrains sterile neutrinos with lifetimes larger than 0.02 s [91].The dotted blue line delimiting from below the blue region in Fig. 10 corresponds to the limit of validity of the assumptions used to obtain cosmological bounds.Indeed, these constraints are derived by considering only sterile neutrinos thermally produced and sufficiently short-lived so that they do not change the nuclear reaction framework by their meson decay products [91].In this context, we mention that cosmological constraints may be extended also to the region of the parameter space in which sterile neutrinos are produced non-thermally (the hatched region below the dotted blue line, labelled as "COSMO OUT OF EQUILIBRIUM"), as discussed for instance in Ref. [101].In this region, however, the bounds are more strongly dependent on the assumptions on the early universe and the extra couplings of these states.In the case of ν 4 mixed with ν τ , the situation of the bounds is qualitatively similar.Factor ∼ 2 differences are due either to the extra production processes for sterile neutrinos associated with charged current interactions with muons, or to the extra decay channels, present only for sterile neutrinos mixed with ν µ .

V. CONCLUSIONS
In this work we revised and improved current bounds on heavy sterile neutrinos mixed with the active ones.In particular, we considered the cooling bound derived from neutrino observations from SN 1987A.We also studied the decays of heavy sterile neutrinos, affecting the SN explosion energy and possibly producing a gamma-ray signal.We improved the characterization of sterile neutrino neutral current interactions of ν 4 with nucleons.We also include charged current interactions of ν 4 with muons, which is relevant for sterile neutrino production mixing with ν µ .Contrary to consolidate belief, it results that the dominant channel for sterile neutrino production is associated with neutral current interactions.Furthermore, we extended the bounds to the trapping regime of ν 4 verified at large mixing angles, adopting the the socalled "modified luminosity criterion".We also strengthened the SN cooling bounds considering (non)radiative decays of heavy neutrinos, and characterizing their effect on excessive energy deposition in the SN envelope, and the observable gamma-ray signal when decays occur outside the SN.The combination of all the SN bounds (together with laboratory ones) allows one to exclude values |U α4 | 2 ≳ 2 − 3 × 10 −11 for m a ≲ 100 MeV.At larger masses the bound tightens to |U α4 | 2 ≲ 10 −13 − 10 −14 in the range 200 ≲ m 4 ≲ 500 MeV.It is worthwhile to mention that another possible astrophysical bound on sterile neutrinos comes from the observation of binary neutron star merger events (see, e.g., Ref. [102] for related constraints on axionlike particles).We reserve this analysis for future work.The most interesting region that SN bounds leave open for future laboratory searches (such as DUNE, SHiP and MATHUSLA) is the range 10 −6 ≲ |U µ4 | 2 ≲ 10 −9 for m 4 ≳ 500 MeV.In the case of mixing with ν τ , DUNE would also have the potential to robustly probe the range of masses down to ∼ 10 MeV at large mixings, where the overlap between SN and laboratory experiments is minimal or absent.
We conclude with two remarks.Having a synoptic view of the bounds following from SN arguments reveals that in most of the parameter space, at least a couple of arguments lead to constraints of similar strength.Since they suffer from different systematics, this is reassuring in supporting the overall reliability of such indirect limits.For instance, diffuse gamma-ray and 511 keV bounds rely on average properties of SN, such as their rate, rather than the single SN 1987A event.Also note that one does not have to rely on the cooling argument, which has been repeatedly criticized in recent years, to derive the strongest bounds from SN for heavy sterile neutrinos.
A similar remark applies on the relation between SN and cosmological bounds.It is reassuring that the bulk of the excluded parameter space overlap.The underlying assumptions in deriving the two classes of bounds are indeed very different.For instance, in non-standard cosmological scenarios with low-reheating temperatures, the BBN bounds can be lifted [103,104].Since in astroparticle physics one cannot control experimental conditions, the accumulation of independent ways to probe a certain type of new physics is essential for a broad acceptance of the robustness of the derived bounds.

FIG. 1 :
FIG.1: Upper panel : Radial profiles of the temperature T (solid black line) and density ρ (dashed line) in the SN core.Middle panel : Radial profile for the degeneracy parameters of neutrons ηn (solid black line), protons ηp (dashed black), electrons ηe (dotted black) and muons ηµ (dot-dashed red).Lower panel : Radial profiles of the electron fraction Ye (solid line) and muon fraction Yµ (dashed line).All panels refer to the post-bounce time t pb = 1 s.

α=μ, m 4 2 L 4 (erg s - 1 )FIG. 2 :
FIG.2: Sterile neutrino luminosity as a function of |Uα4| 2 for α = µ, m4 = 10 MeV and m4 = 120 MeV (black and red lines, respectively) and the same for α = τ (orange and blue lines, respectively).The horizontal dotted line corresponds to the limit value of Lν = 3 × 10 52 erg s −1 .Dashed lines are used for values of the mixing where the sterile neutrino luminosity exceeds the luminosity of the species να mixed with ν4.

FIG. 3 :
FIG.3: Branching ratios for the relevant decay channels listed in Tab.II as a function of the sterile neutrino mass m4 for a mixing with νµ (upper panel) and with ντ (lower panel).

For ν 4
mixed with ν µ , SN arguments lead to the lower limit |U µ4 | 2 ≳ O(10 −11 ) for m a ≲ 100 MeV, notably due to the 511 keV line argument (see Sec. III C) and the diffuse gamma-ray flux (see Sec. III E).At larger masses, the SN bound tightens to |U µ4 | 2 ≲ 10 −14 in the range 200 MeV ≲ m 4 ≲ 500 MeV, due to the absence of gamma-rays in coincidence with SN 1987A (see Sec. III D).The existing laboratory bounds nicely complement the SN ones, excluding the parameter space all the way to large mixing angle, and overlapping with SN bounds here dominated by the explosion energy argument (see Sec. III B).Future laboratory experiments are expected to charter new parameter space only for m 4 ≳ 500 MeV, probing mixing angles 10 −6 ≲ |U µ4 | 2 ≲ 10 −9 .

TABLE I :
[63]red matrix elements for sterile neutrino scattering processes (assuming mixing with the species α, and β ̸ = α), summed over initial and final states, where gL = − 1 2 + sin 2 θW , gR = sin 2 θW[63].The symmetry factor S = 1/2! is already included when two identical particles are present in the same state (second row).The particles involved in each reaction are enumerated as 1 + 2 ↔ 3 + 4. The last two processes are valid only in the case of mixing with the muon neutrino, α = µ.The terms |M| 2 AA , |M| 2 V A and |M| 2