Medium effects in supernovae constraints on light mediators

In this article, we reevaluate supernovae (SN) constraints on the diffusion time of neutrinos for a family of extensions of the Standard Model that incorporate new light scalar and vector mediators. We compute the neutrino mean free path, taking into account medium effects in the neutrino-nucleon scattering cross-section, and a radial dependence of the density, energy, and temperature inside the proto-neutron star to determine the coupling strengths compatible with SN1987A constraints on the time duration signal of diffusing neutrinos. We show that medium effects can induce an order of magnitude enhancement in the neutrino mean free path with respect to the vacuum calculation. The increase is more significant when new physics terms dominate over the Standard Model contribution (that is, for small mediator mass and large couplings). Finally, we interpret these results as bounds on the parameter space of a vector $U(1)_{B-L}$ model and scalar lepton number conserving and lepton number violating scenarios, improving on previous results in the literature where medium effects were ignored. We show that SN constraints on the neutrino diffusion time lie within regions of the parameter space that are already ruled out by other experimental constraints. We also comment on potential limits due to changes in the SN equation of state or right-handed neutrino free-streaming, but argue that detailed numerical simulations are needed to improve the reliability of these limits.


I. INTRODUCTION
There is growing interest in particle physics constructions that feature new low-mass mediators as an extension of the Standard Model (SM). This interest has been fuelled in part by observed anomalies in low-energy observables, such as the long-standing discrepancy in the muon anomalous magnetic moment, and in part by the potential for these models to connect to an otherwise secluded sector that could account for the dark matter content of the universe. In the absence of any clear signal of new physics either from the high-energy regime explored at the LHC or from direct detection searches for canonical weakly-interacting massive particles, alternative search strategies are gaining attention, including those that favour low-energy precision measurements across a wide range of experimental techniques.
In general, particle models with light mediators include new interactions in the neutrino sector, which can often be interpreted in terms of an effective theory and described as non-standard interactions (NSI). This is a particularly interesting possibility, as it widens the range of experimental probes of these constructions. For example, data from dedicated neutrino detectors can be used to set constraints on the resulting neutrino-electron scattering, as no deviation has yet been found from the SM prediction. Experimental data from GEMMA [1] and Borexino [2] lead to competitive bounds on the parameter space of these models for mediator masses of the order of the MeV [3,4]. Likewise, new physics in the neutrino sector can also be probed through coherent elastic neutrino nucleus scattering (CEνNS). This rare process, which in the SM occurs at leading-order through the exchange of a Z boson [5], has been recently observed by the COHERENT collaboration, employing two different targets. Their original results on CsI [6] and the latter measurement in liquid argon [7] are in agreement with the SM prediction, which allows to set limits on the parameter space of new-physics scenarios with light mediators [4,[8][9][10][11][12][13]. The large exposure of projected future experiments at spallation neutrino sources [14][15][16] makes them well suited to probe mediator masses in the MeV range.
Dark matter direct detection experiments have recently joined the quest for light mediators. Although they were not specifically designed for this aim, these detectors can be sensitive to neutrino interactions with the electrons and nuclei of their target. In fact, due to the resemblance of this signal with that expected from a weakly-interacting massive particle, CEνNS is considered a background in the search for dark matter particles [17], and it can be interpreted as a 'neutrino floor' in the dark matter-nucleon scattering cross section [18]. New physics in the neutrino sector can alter the SM predictions for nuclear recoils [19,20], thereby raising the neutrino floor especially at low dark matter masses [21,22]. Likewise, the effects of new physics could also be felt in the electron recoil spectrum. So far, the absence of a confirmed signal in either electron and nuclear recoils has been translated into constraints on simplified models and scenarios with non-standard neutrino interactions [4,19,[23][24][25][26].
The presence of light new states can also have crucial implications in cosmological and astrophysical observations. In particular, a bound on the effective number of relativistic degrees of freedom can be derived from Big Bang nucleosynthesis, that sets stringent constraints on the abundance of new particles with masses below 1 MeV [27].
Finally, light particles can be produced in stellar interiors and alter different aspects of stellar evolution, in particular, SN explosions. The observed neutrino flux from SN1987A [28][29][30] permitted the study of neutrino properties (setting upper limits on their mass and lifetime) and set bounds on new physics in this sector [31][32][33][34][35][36][37][38][39][40][41][42][43]. Most of these constraints address the cooling effect of light weaklyinteracting particles, which would escape the star unimpeded (thus leading to a different neutrino luminosity than the one actually observed). This is applicable to a wide range of models that include axions (or axion-like particles) as the most representative example. It also affects new mediators in the neutrino sector, for example light gauge bosons (dark photons) [44][45][46][47][48] or scalar particles [49,50]. It should be noted that cooling constraints generally restrict a window of very small couplings (as the particle needs to escape the star and these bounds do not apply if it is reabsorbed after production), however, large couplings can affect the diffusive properties of neutrinos [51], decreasing their mean free-path and leading to an emitted flux which can exceed the observed ∼ 10 s. This has been used in Refs. [42,43] to derive upper bounds in the neutrino coupling to new scalar and vector mediators, showing that they can compete with limits from other experimental sources. As it was then shown in Ref. [21], these constraints are relevant, since they determine the maximum height of the neutrino floor in direct dark matter detection experiments.
Contributions to neutrino-nucleon scattering for a scalar mediator (left) and a vector mediator (right).
In this article, we carry out an improved computation of the neutrino mean free path to derive a more reliable upper limit on the mediator coupling to neutrinos. To do this, our analysis takes into account the radial dependence of the density, energy, and temperature in the proto-neutron star (proto-NS). More importantly, given that the density can be several times the nuclear saturation density, ρ 0 ∼ 2 × 10 14 g cm −3 , we include medium effects in the computation of the neutrino-nucleon scattering cross-section, which were not included in the analysis of Refs. [42,43]. Medium effects effects can increase the neutrino mean free path by approximately an order of magnitude [51][52][53][54][55] with pure SM interactions, and in this work we extend the analysis to new physics, showing that they can be even more important for low mass mediators. This can significantly alter the resulting neutrino diffusion time. We interpret our results as upper bounds on the neutrino couplings in simplified models, for which we choose a (lepton number violating and lepton number conserving) scalar and a U (1) B−L vector model. This article is organised as follows. In Section II we introduce simplified models that incorporate a new vector or scalar coupling to the SM and provide the new matrix elements that contribute to neutrino-nucleon scattering. In Section III we review the main constraints on light new mediators that can be derived from SN physics. We explain in detail how to compute the neutrino mean free path, for which we consider medium effects and a radial dependence of the density and temperature. In Section IV we show the medium effects on the neutrino mean free path and study the area of the two-dimensional parameter space that is excluded by SN constraints, putting it in context with bounds from other detection methods. Finally, our conclusions are presented in Section V.

II. SIMPLIFIED MODELS OF NEW PHYSICS WITH LIGHT MEDIATORS
We consider two classes of low-scale simplified models, which introduce a new light mediator (scalar or vector) that provides extra interactions between neutrinos and SM fermions (quarks and/or leptons) [56]. We do not address here how to realise these scenarios as complete UV models, as this might imply adding extra fields for anomaly cancellation. The new interactions induce corrections to neutrino-nucleon scattering as illustrated in Fig. 1, thereby altering the neutrino diffusion time. In this section we compute the new matrix elements that enter the computation of the neutrino mean free path.
In the following, we denote p µ = (E N , p ) and p µ = (E N , p) as the four-momentum for the outgoing and incoming nucleon of effective mass m * N , respectively, and k µ = (E ν , k ) and k µ = (E ν , k) the analogous for neutrinos of mass m ν . The transferred four-momentum is defined as q µ = (q 0 , q) = p µ − p µ = k µ − k µ . We will give more details about the concept of effective mass in Section III.

II.1. Scalar mediator
As a first example, we introduce a new light scalar field, φ, of mass m φ that interacts with SM fermions. The scalar must necessarily couple a left handed fermion field to a right handed one, which can lead to two different scenarios for neutrinos, which we study separately.
First, we can consider a lepton number violating (LNV) scenario with Majorana neutrinos. In this case, we assume that the masses of possible heavy neutrino states are much heavier than the typical energies inside the supernova and their interactions can be disregarded. Alternatively, assuming Dirac neutrinos, the scalar field might couple to right handed neutrinos in a lepton number conserving (LNC) model. In this case, different regimes are possible for the right handed neutrino mass, m ν R . In this work we consider that m ν R O(1 eV), and therefore the relevant cross sections for neutrino scattering mediated by the new LNC scalar in most terrestrial experiments are nearly identical to the LNV case 1 .
In the case of an LNV scalar, the new interaction terms in the Lagrangian take the form where C v , C q , C l and are the couplings to neutrinos, quarks, and leptons, respectively. This differs from the LNC scalar only in the first term, with Within SN, we consider neutrinos scattering coherently with nucleons. The quark contribution to the Lagrangian can therefore be rewritten in terms of whole nucleons as with where values of m q and f N q can be found in Ref. [42]. For universal Q q , as considered in our work, one obtains C N = 13.8 C q . The presence of the additional scalar provides a new channel for neutrinonucleon scattering. The squared matrix element is modified from its SM value by an additional term with In principle, C ν and C N are separate free parameters in our model. We can reduce the dimensions of the parameter space by placing constraints on the effective coupling Y = √ C ν C N . Any other process that involves neutrino-quark scattering can directly set bounds on Y , and so can be immediately compared with our constraints from SN physics. Limits on the individual couplings C ν and C q can also be combined to obtain a constraint on Y . However, applying bounds from other processes require us to make assumptions about the relation between the couplings to different particles.
Constraints from neutrino scattering (either CEνNS or neutrino-electron scattering) in terrestrial experiments are practically identical for both a LNC and a LNV scalar (with m ν R 1 eV). However, each species of scalar has its own unique constraints as well, stemming from limits on either the production of right handed neutrinos, or the rate of LNV processes. Furthermore, the rate of neutrino scattering within the dense interior of the SN differ for LNV and LNC scalars due to the different nature (Majorana and Dirac, respectively) of neutrinos and to the fact that LNV interactions induce changes in the chemical potentials of the particles in the medium, as it will be discussed in Section III.2.
, the neutrino-scattering cross section would be kinematically suppressed when the centre-of-mass energy of the collision is smaller than the ν R mass, and right handed neutrinos would not be produced in any of the scenarios that we consider. For the intermediate case where the ν R mass is of the order of the MeV, scattering would still occur at some or all of the energy scales we consider, with the presence of the massive ν R producing potentially measurable effects that could even allow us to determine mν R . This would produce a complex phenomenology which, while interesting, is beyond the scope of this work.
II.2. Vector mediator: Gauged U (1)B−L As a paradigm of vector-mediated scenarios, we consider a U (1) B−L model. The symmetry is spontaneously broken, leading to a new vector mediator, Z , with mass m Z , which, depending on the scale of the symmetry breaking, can be much lighter than the SM Z boson.
The new interaction terms in the Lagrangian take the form where all quarks have the same charge, C q = g B−L /3, and all charged and neutral leptons have charge C l,ν = −g B−L . The Z may also interact with SM fermions via kinetic mixing with either the SM Z or the photon. Here we assume that the couplings to SM fermions induced through kinetic mixing are subdominant to the direct couplings C i , and so we neglect kinetic mixing effects. Gauging the U (1) B−L symmetry in a theory with only SM fermions leads to non-zero anomalies, so additional particles must be added to ensure these anomalies are cancelled out. Although several nonminimal solutions exist, the simplest way is to include three right handed neutrinos [57]. In this case the anomalies from the left and right handed neutrino loops exactly cancel and the theory is anomaly free. As in the case of the LNC scalar considered in the previous section, neutrinos in this model are Dirac fermions and interactions with the U (1) B−L mediator do not violate lepton number. There is thus no mechanism for the efficient conversion of electron neutrinos to neutrinos and antineutrinos of other flavours which can deleptonize the star or impact properties such as the entropy or the temperature of the SN core. Unlike in that model, right handed neutrinos are not produced in t-channel neutrino scattering events mediated by a U (1) B−L mediator. We therefore assume that right handed neutrinos are heavy enough so as not to be produced inside the SN (this requires masses above the GeV scale, natural in most see-saw models). This choice has no effect on the CEνNS or neutrino-electron scattering cross sections in the U (1) B−L model.
Considering this scenario, the total matrix element contains an interference term between Z and the SM Z boson, Since neutrinos are Dirac, only left handed neutrinos and right handed antineutrinos couple to the Z and Z vector mediators. Remember also that we are only concerned with the light neutrino states, as the heavy ones are assumed to have masses well above the GeV and are therefore not produced in the SN. Both for the left handed neutrino nucleon interaction, νN → νN , and for the right handed antineutrino nucleon one,νN →νN , the pure BL term reads On the other hand, as the SM amplitude is different for νN → νN andνN →νN , the mixing term is different too. For νN → νN it reads and forνN →νN In these expressions,F N 1 (q 2 ) is the neutral current vector form factor andF N A (q 2 ) the axial one. Since the interactions considered take place inside the SN core with a temperature k B T m * N (where k B the Boltzmann constant that we set to k B = 1 from now on), for the scenario that we study |q 2 | 0.15 GeV 2 and we can approximateF Fig. 4.9 of [58], where the dependence of the form factors on |q 2 | is shown, to understand better this approach. This is equivalent to say that at these energies neutrinos do not see the quark structure but the nucleon as a whole. Therefore, for our calculations we takeF n

III. SUPERNOVAE CONSTRAINTS ON NEUTRINO INTERACTIONS
Neutrinos play a key role in the first minutes after the collapse of a massive star. After these stars explode as SN, neutrinos are copiously produced and trapped in a core of radius ∼ 10 km during the first ∼ 10 s after the explosion [59,60]. At this time, the proto-NS core reaches densities several times that of nuclear saturation density, ρ B ∼ (2 − 3)ρ 0 , and temperatures T 50 MeV [61]. These neutrinos (and antineutrinos) arise mostly from the electron capture on protons and flavour equipartition, i.e. electron (ν e ), muon and tau (ν τ , ν µ ), collectively ν X . Matter in the core is constituted of free nucleons, mostly neutrons, with which hot (thermalized) neutrinos interact many times before decoupling from matter inside the star. In this way, it is reasonable to think that new physics in the neutrino sector may alter SN dynamics in several ways. In particular, light mediators that couple to both neutrinos and nucleons could affect the diffusion time of neutrinos and their energy spectrum, which were first known from the measurements of the SN 1987A burst [28][29][30].
New neutrino interactions can also modify the equation of state (EoS) of the proto-NS core [42]. This is particularly relevant if the new couplings between neutrinos and nucleons violate lepton number, as these processes can remove ν e and convert them intoν e , ν µ ,ν µ , ν τ andν τ . It has been stated before [62][63][64] that this type of interaction would imply sizeable differences in physical conditions, in particular in the entropy, lepton fraction and temperature, but these effects do not affect the explosion, as it is explained in Ref. [64]. In order to fully assess the effects of changes in the EoS, one must resort to numerical simulations, which is beyond the scope of this work. We therefore do not implement these bounds, but we show their potential reach.
The effects described above rely on the determination of the neutrino mean free path. As we argue in this section, this is a density and temperature dependent quantity, a fact that cannot be ignored in order to derive reliable bounds on new physics. In the following, we explain how this has been implemented in our calculation.

III.1. Mean free path in a nuclear medium
The mean free path in a dense and hot medium must be obtained from a self-consistent treatment, taking into account the distribution functions of all the particles of the interaction in the incoming and outgoing kinematic phase space at finite density and temperature.
The standard differential cross-section for the scattering between (anti)neutrinos, (ν) ν, and nucleons, N , can be written as where the phase space volume element is The Fermi-Dirac distribution functions for the ith-type particle (protons and neutrons), f i (E i ) = 1/1 + e (Ei−µ * i )/T , are written in terms of the effective mass and chemical potential m i and µ * i , which differ from the naked values in vacuum by the presence of meson fields [65]. Here, E i is the nucleon energy and T the temperature of the medium. The distribution function for the neutrino, f ν (E ν ) = 1/1 + e (Eν −µν )/T , depends on its Majorana or Dirac nature. Note that the contribution from the scattering of neutrinos and antineutrinos off electrons is subleading [60] and is therefore not included in our analysis.
Given the high values of the density and temperature that can be reached in a SN core, matter effects cannot be ignored when computing the neutrino mean free path. Vacuum expressions do not capture the rich physics of the many-body effects, which can increase the mean free path by more than two orders of magnitude [51,[66][67][68]. In our calculation, we have included the Fermi-Dirac distribution functions, which partially restrict the outgoing phase space, and the effective values for the nucleon mass (being m N < m N ) and chemical potential, thereby improving the results of previous works [42,43].
The neutrino mean free path can be expressed as the inverse of the cross section per unit volume, λ = V /σ. The differential cross section in Eq. (12) has to be weighted with the phase space volume for the incoming nucleon sector, Neglecting the neutrino mass, m ν ∼ 0, and performing a partial integration (see e.g., Ref. [67] for more details), we obtain where F is the Fermi blocking term We denote θ as the polar angle between q and p, and θ ij and φ ij as the polar and azimuthal angles between the momentum of particle i and particle j (we take i = 1, 2 for the incoming neutrino and nucleon, respectively, and i = 3, 4 for the outgoing ones). The angle θ 0 is fixed by energy conservation, and the flux term can be expressed as The polar angle, θ 12 , between the incoming neutrino and nucleon is obtained by solving After doing this, Eq. (14) can be rewritten as with Note that the limits of the integral over | p| comes from imposing | cos θ 0 | ≤ 1.
The density profile of the proto-NS and the temperature are radially dependent functions. In the centre it is well above the nuclear saturation density and decreases towards the outer regions. As a consequence, the neutrino energy is also a function of the radius, and it changes within the hot and dense medium. A complete treatment should also include the time-dependence of these quantities, but this is only possible with a dedicated numerically simulation, which is beyond the scope of this paper. As a more modest approach, we use the results from Ref. [69], where the radial and time dependence of the proto-NS temperature, density, etc, were computed for an 18 M progenitor, considering the TM1 model [70,71] to describe the physics of nuclear matter. As we are most interested in processes during the Kevin-Helmholtz phase, we have used results at times of 1 and 5 seconds post-bounce. TABLE I. Values of neutron effective chemical potential, µ * n , proton effective chemical potential, µ * p , neutrino chemical potential, µν , and nucleon effective mass, m N , for the spherical shells (labeled by the index k and defined by an outer radius R) that we consider at the two time snapshots of 1 s and 5 s, with a baryonic density, nB, temperature, T and electron fraction, Ye, in our analysis for Dirac neutrinos. In the last two columns the resulting neutrino mean free path in the SM, λ SM νe , and antineutrino mean free path, λ SM νe , are also shown. Temperatures, densities and electron fraction are taken from Ref. [69]. Since the density decreases slowly within the core, we have considered spherical shells where temperature and density are taken to be constant, with values extracted from Fig. 7 of Ref. [69] as detailed in Tables I and II. We consider these shells from the centre of the star to the neutrino sphere radius, R ν , where neutrinos decouple from matter. The neutrino sphere radius is a function of the neutrino energy, E ν , flavour and thermodynamical state of the matter [60,69], and therefore depends on the time after bounce. Knowing the antineutrino spectrum from SN 1987A, which suggests that neutrinos decoupled with energies around a dozen MeV or, accordingly, temperatures of a few MeV (see Fig. 1 in Ref. [60]), the neutrino sphere can be determined from the following condition on the optical depth, τ ν ,

Majorana Neutrinos
Neutrinos with purely SM interactions (i.e., the electroweak force) already fulfil the opacity constraint and they remain diffusive while streaming out. Adding a new physics contribution to neutrino interactions increases, in general, the probability of scattering, and therefore the condition of Eq. (20) is satisfied. The addition of a new vector mediator, however, must be considered more carefully, as the destructive interference with the Z boson can dominate in certain regions of the parameter space. We have checked explicitly that neutrinos remain diffusive in all the parameter space.
Regarding the neutrino energy, we consider Dirac neutrinos to be in thermal equilibrium with an energy E ν = µ ν + πT [61], where µ ν is the neutrino chemical potential. Majorana neutrinos do not reach thermal equilibrium (and do not have a chemical potential associated to them). Their energy should be determined numerically, but at the time of writing this article, we are not aware of any hydrodynamic simulation of SN with Majorana neutrinos. Therefore, we consider two limiting cases for the Majorana neutrino energy, namely E ν = πT (a relation that has been used, e.g., in Ref. [63]) and E ν = µ D ν + πT , where µ D ν is the chemical potential for Dirac neutrinos. In order to determine the effective nucleon mass and chemical potentials in each shell, we consider the TM1 model (also used in Ref. [69] to obtain the radial temperature and density profiles that we use in Tables I and II). We use a relativistic mean field (RMF) approach where baryons are considered as Dirac quasiparticles moving in classical meson fields and the field operators, φ, are replaced by their expectation values, φ . The TM1 [70,71] model, widely used in current numerical simulations, is a representative example where the set of parameters used can smoothly connect low and high density regions in the dynamical stellar description. The presence of an effective nucleon mass and effective chemical potential is due to the non vanishing values of the Lorentz scalar meson, σ , Lorentz vector, ω µ , and vector-isovector, ρ µ , meson fields. The self-consistent solution in the RMF approximation of the meaningful components of these field values (σ, ω 0 , ρ 03 ) was obtained in Refs. [70,71]. In this way, an effective nucleon mass, m * N = m N − g σN σ , and effective chemical potentials, µ * i = µ i − g ωN ω − g ρN t 3i ρ , are obtained for each thermodynamical state. The quantities g σN , g ωN , and g ρN are dimensionless constants that couple nucleons to the σ, ω, and ρ mesons, respectively. t 3i is the third component of the isospin of the proton or the neutron, i = p, n. Relativistic baryonic energies are defined as E * i (k) = k 2 + m 2 N . This parameter set includes self-interaction terms from scalar, vector and vector-isovector mesons in non isospin symmetric nuclear matter at finite temperature. TM1 interaction terms are constrained by the nuclear masses, radii, neutron skins and their excitations. When applied to the derived proto-NS, the mass-radius diagram allows to fulfil the subsequent two solar mass constraint from recent observations of older objects.
By imposing baryonic charge conservation, fixed lepton fraction and charge neutrality, the RMF equation for the fields, σ, ω, and ρ are solved to obtain effective nucleon masses and chemical potentials at finite temperature [70]. Note that at finite temperature a non-vanishing positron fraction is available and a net lepton number (including the neutrino sector) arises. The resulting values of m N and chemical potentials are shown in Tables I and II. For reference, we also indicate the resulting neutrino mean free path in the SM, λ SM . For the Majorana case, we show the two values that correspond to the limiting cases for the energy window we consider, E ν = πT and E ν = µ D ν + πT . Finally, using these values of the effective masses and chemical potentials in Eq. (19), we compute the neutrino and antineutrino mean free paths for each spherical shell, λ k , to determine the diffusion time in each of them. We add these to estimate the total duration of the neutrino (antineutrino) emission as where R k is the outer radius of the k-th spherical shell, and R 0 = 0. Since the neutrinos detected in the signal from SN1987A were electron antineutrinos we compare the total duration of the antineutrino emission with the observed ∆t ∼ 10 s.

III.2. Effects on SN dynamics
During the core collapse, protons and electrons combine into neutrons, with each of the O(10 57 ) reactions producing an associated electron neutrino but also other flavours (muon and tau) and reactions come into play including electron-positron annihilation, or nucleon bremstrahlung [72]. Many of these neutrinos and antineutrinos escape the developing neutron star, carrying away the largest portion of the gravitational energy lost in the collapse. This neutrino burst was detected by the experiments Kamiokande II, IMB, and Baksan in coincidence with the SN 1987A event [28][29][30] and their measurement allows us to place constraints on models of new physics which would affect the propagation of the neutrino burst through the proto-NS. In particular, the diffusion time scale ∆t of Eq. (21) can be altered by new interactions between neutrinos and nucleons. One must impose that this quantity is of the order of the observed duration of the burst, ∆t 10 s.
When we apply Eq. (21) using pure vector-axial SM interactions, we obtain ∆t SM = 2.6 s for Dirac neutrinos and ∆t SM = 2.4 s for antineutrinos. For Majorana neutrinos, we obtain ∆t SM = 1.1 s (∆t SM = 2.7 s) for E ν = πT (E ν = µ D ν + πT ). For consistency, these numbers have been obtained using the configuration for t ∼ 1 s of Tables I and II 2 .
The mean free path for Majorana neutrinos is extremely sensitive to the incoming neutrino energy. If the same energy as for Dirac neutrinos is used, Majorana and Dirac neutrinos yield a very similar λ SM k . The main difference between these two scenarios comes from the energy distribution of neutrinos in the star, and if E ν = πT is used, the neutrino mean free path for Majorana neutrinos increases considerably (by up to a factor ∼ 4, depending on the SN shell). To obtain a more precise answer, neutrino energy should be determined by numerical simulations, especially for Majorana neutrinos, which do not reach thermal equilibrium.
We should point out that we have not considered correlations that typically enhance the mean free path, see for example Ref. [66], where random phase approximation correlations indicate that in-medium instabilities cause an increased mean free path. In addition, the inclusion of weak magnetism effects can alter these values by approximately 20%. A fully consistent treatment of the dispersion suffered by neutrinos would require the use of Boltzmann transport and energy conservation at each interaction step. Since this is not feasible in our calculation, we rely on semi-analytical estimates to include the size of corrections expected when new physics is involved. This, together with the choice of a pure axial-vector current for the SM interaction description (realistic simulations consider extra tensorial terms [66]) might explain why our results for the SM diffusion time is slightly shorter than 10 s.
After new physics contributions are included, we need to ensure that the total neutrino diffusion time, ∆t, remains consistent with the duration of the observed SN1987A neutrino burst, which leads to the condition ∆t 10 s, setting an upper limit on the strength of neutrino interactions with quarks. Since what experiments have observed is the burst of antineutrinos through the measurement of positrons coming fromν e n → pe + , we apply this constraint to the antineutrino diffusion times. Lacking a full simulation analysis, the systematic error of this numerical approach is difficult to estimate. Since the SM times that we obtain are slightly shorter than the 10 s observed, we also indicate the region in which the new physics contribution becomes as important as the SM one, which can be estimated as follows, For consistency, the bound given by Eq. (22) is computed using the t ∼ 5 s configuration of Tables I  and II, whereas for the one of Eq. (23) we must use the t ∼ 1 s configuration, so that the stellar properties correspond to those at half the amount of time that neutrinos typically spend in the supernova. If the interaction between neutrinos and quarks via a scalar mediator is LNC, right handed neutrinos with masses m ν R 1 eV can be produced in the scattering of left handed neutrinos off SM fermions. Once ν R are produced, they free stream out of the star since they do not interact via SM interactions, leading to a suppressed ν L flux and a much shorted cooling time. This can be prevented in two ways, either ν L do not interact with SM fermions via a new physics interaction while they are trapped in the proto-NS, i.e., with where λ NP k the mean free path computed only with new physics interactions. Alternatively, if the new physics scattering rate is sufficiently large, the ν R do not free stream out of the core. Instead, they scatter back into ν L , and the overall effect is that the neutrinos are constantly switching back and forth between ν L and ν R as they diffuse out of the SN. In line with the analysis of Ref. [42], we constrain this effect by requiring that at least 100 new physics interactions take place over the path the ν R would take when free streaming out of the SN, this is, where ∆R k = R k − R k−1 . Note that, in this case, neutrinos spend on average half of their time inside the star as ν R and not feeling electroweak interactions. Thus, the diffusion time is calculated as The condition in Eq. (26) must be understood as a qualitative statement, limiting the relative size of the new physics and SM contributions. 3 Finally, as we have mentioned above, if new neutrino physics happens via a LNV interaction, quantities as the entropy, lepton fraction and temperature can be greatly affected. Qualitatively speaking, one can estimate when neutrino interactions might have an impact on the EoS of SN matter by requiring that they scatter at least once as they stream out of the star. This minimal condition can be imposed as in the LNC case, requiring that Eq. (24) is fulfilled.
Changes in the EoS of SN matter would only indirectly affect the neutrino signal: either through a change of detected flavours or time correlation on Earth. Currently, preliminary microscopic simulations have been performed using one neutrino species allowing internal deleptonization with σ/σ SM 10 −3 for the reaction ν e +N →ν e +N [73]. It should be noted that this work, based on one single flavour LNV reactions, employ one-dimensional codes where luminosities are monitored for moderate times after the trigger of the stellar collapse. Because of their simplified microphysics input, they do not display all the physics features relevant to the dynamics and outcome of the observable signals at Earth, such as convection [74] and pre-collapse seed perturbations [75,76]. The results of Ref. [73] suggest some degree of insensitivity of the overall dynamical evolution of the models to rather dramatic modifications of the microphysics, but to obtain a reliable result, it is necessary to extend the simulation to larger times.
More recent dedicated studies have performed detailed numerical simulations that include other exotic species, such as sterile neutrinos with O(50) MeV masses that decay into SM neutrinos. It was found that these could have a significant impact on the dynamics of the core collapse. In fact, for some parameter choices, they can lead to too energetic explosions thus indicating that with the help of the astrophysics data it is possible to further constrain the parameter space of sterile neutrinos [77]. Given the proto-NS star conditions described in our contribution we anticipate that the early deleptonization could be very different from that predicted using SM neutrino transport. Further work on the computational side is needed in order to obtain the full picture regarding how the final spectra would look like as measured in a hypothetical neutrino detection on Earth.
Because of the arguments above, and in line with previous literature [42], we only indicate when the condition given by Eq. (24) is fulfilled for an LNV model in Fig. 4, but we do not consider it a solid constraint yet. In order to obtain robust results, a full hydrodynamic simulation would be needed that relates changes in the EoS with SN observables.
Previous works [42,43] have implemented the conditions of Eqs. (22) to (26) to constrain new physics models, but ignoring the Pauli blocking for the outgoing nucleons and their effective masses and chemical potentials. In our work, we use the results of the previous section, where medium effects have being taken into account, for the calculation of the mean free path in Eq. (19) to derive a more consistent estimate.

IV. DISCUSSION
In this section, we compute the neutrino mean free path for the new physics models of Section II and show the constraints on the corresponding parameter spaces.
First, in order to quantify the importance of medium effects on our results, we compare the values of the mean free path in the k-th shell of the star including medium effects (λ k , obtained using Eq. (19)) with those found neglecting both the Pauli blocking of the outgoing states and the effect of the dense nuclear medium on the nucleon mass (λ 0 k ). To compute λ 0 k one must include the density of nucleons in each shell. To mimic vacuum calculations, where λ vac ∼ (σn B ) −1 , we perform the integral over the phase space volume for the incoming nucleon sector, n B ∼ 2d 3 pf N (E N )/(2π) 3 , retaining the effective mass of the nucleon only in this factor for consistency with the baryonic density of each shell (note that the effective mass is related to the baryonic density through the RMF equations). Since we are keeping f N , λ 0 k is not truly a vacuum mean free path and it incorporates some in-medium effects. In Fig. 2, we plot the ratios λ k /λ 0 k as a function of the mediator mass for two representative choices of the coupling in both the U (1) B−L and LNC scalar models, and for the five radial shells at t = 5 s described in Table I. We can observe that medium effects are more important in the inner shells (where the temperature and density are higher) and decrease as we move outwards in the SN. The neutrino mean free path can increase by approximately one order of magnitude with respect to the vacuum estimate in the central region of the SN, an increase which is more pronounced (up to a factor 14) with new physics contributions. In each plot, the lines become horizontal in the limit where the SM dominates over the new physics contributions, which happens at small values of the couplings and large mediator masses.
The non-trivial behaviour of λ k /λ 0 k when new physics contributions dominate is due to the dependence of the scattering amplitudes of Eqs. (6) and (9) to (11) on the effective nucleon mass and mediator mass. The mediator mass enters through the denominator of the corresponding propagator (1/(q 2 − m 2 φ ) or 1/(q 2 − m 2 Z )), affecting the range of the q 2 values that enter the integration in Eq. (19), which also multiply terms in the numerator with the effective nucleon mass. For heavy mediators (m 2 Z ,φ |q 2 |), this asymptotically reaches a flat behaviour. This can be seen more clearly in Fig. 3, where we represent the ratio λ k /λ 0 k , computed separately for each contribution (new physics, interference term, and SM) for the innermost shell (a quantity that we define asλ 1 /λ 0 1 ). The non-trivial behaviour observed in Fig. 2 appears when all the contributions are included in the calculation of the mean-free path.
Using the values of the mean free path with medium effects, λ k , we compute the neutrino diffusion time for each point in the the parameter space of our low-mass mediator models through Eq. There are numerous experiments which have probed different aspects of new physics in the neutrino sector (or new light mediators in general), without having found so far any significant deviation with respect to the SM predictions. This has lead to constraints on different combinations of couplings in the new physics model. In the U (1) B−L model, all of these constraints can be translated directly into bounds on g B−L . We take our preexisting constraints on this model from Refs. [4,78]. In the two scalar models, we consider three classes of constraints, divided based on the combination of couplings they apply to. In Figs. 4 and 5, we plot all of the constraints under the assumption of universal couplings to SM fermions (C q = C l = C ν = Y / √ 13.8), though they can be rescaled for other models with different relative couplings. • Neutrino-nucleus interactions.
The constraints we derive from supernova neutrino scattering can only be directly compared with other limits on Y = √ C ν C N . This is the same combination of couplings that appears in the cross section for CEνNS. This elusive process has been measured by two experiments within the COHERENT collaboration, the first time in CsI [6], and again more recently in the CENNS-10 liquid argon experiment [7]. Both results were compatible with the SM prediction, and can therefore be used to place bounds on modifications to the CEνNS scattering rate. Constraints on a light scalar mediator were first derived in Ref. [42,79], and were later improved using an updated quenching factor [10,11]. They are identical for the LNV and LNC scalar models. Constraints can also be obtained on Y by combining limits on the individual couplings C N and C ν . The former have been derived from measurements of neutron scattering [80], while the latter can be obtained from cosmology and, in the case of the LNV scalar, from searches for neutrinoless double-beta decay. Cosmological bounds arise from limits on the effective number of neutrino degrees of freedom during Big-Bang nucleosynthesis (BBN), and differ for the LNV and LNC models. All of these are discussed in more detail in Ref. [42]. We represent these constraints with solid lines in Figs. 4 and 5.
Constraints can be derived on the combination of couplings √ C ν C e from measurements of neutrino-electron scattering. In Figs. 4 and 5, we use dot-dashed lines to represent the equivalent bounds that would apply to Y for the case where C e = C q . These bounds can be adapted for a theory with a different relation between these couplings, by rescaling them by a factor C q /C e . We have obtained constraints from the results of the Borexino experiment [2], following the method used in Ref. [4]. We have also included the limits that can be derived from the GEMMA detector, using the results of Ref [26], and from the recent analysis of electron recoil events at XENON1T [4,26,81,82]. In both cases these constraints are less relevant than those from Borexino.
Constraints have also been placed on models with new light mediators from the decays of mesons containing heavy quarks. In Ref. [83], measurements of branching ratios of B + -meson decays   were used to place limits on a new scalar mixing with the SM Higgs. The constraints in that work were obtained due to the effective coupling generated between a top quark (in a loop) and a µ + µ − pair. This constraint can be rescaled to give a constraint on the combination of couplings C t C µ in our scalar models. The resulting constraints on Y are shown as dotted lines in Figs. 4 and 5 under the assumption of universal couplings to SM fermions.
As with the neutrino-electron scattering bounds, the constraints derived from B + -meson decays only apply to Y in a model-dependent way, and should be rescaled with the relative C t and C µ couplings in a scalar model with non-universal couplings. However, even in the case where there is no tree-level coupling to top quarks, a constraint could still be obtained by replacing the top quark in the loop with either its first-or second-generation counterpart. Such a scenario would require more careful consideration than the simple rescaling used for the neutrino-electron scattering constraints.
In Figs. 4 and 5, and Fig. 6 we show in solid lines and under the labels ∆t = 10 s and ∆t = 2∆t SM the upper bounds based on the diffusion time constraints given by Eqs. (22) and (23) for the LNV (blue and orange) and LNC (blue) scalar mediator and for the U (1) B−L model (blue). In Fig. 4 the orange lines correspond to the case in which E M ν = µ D ν + πT while the blue ones correspond to E M ν = πT . As it has been already shown in Section III.2, the SM diffusion time is smaller for lower values of E ν . This translates into a bigger gap between the ∆t = 10 s and ∆t = 2∆t SM lines when E M ν = πT (as aforementioned, the realistic neutrino energy for Majorana neutrinos should be obtained through numerical simulations). The differences between the diffusion upper bounds for the LNV, Fig. 4, and LNC, Fig. 5, scenarios for the same incoming neutrino energy, stem from the fact that medium effects are different for Majorana and Dirac neutrinos since in the first case µ ν = 0. Note that these constraints have been weakened with respect to the ones of [42], mainly due to the effect of the Pauli blocking which restricts the phase space of nucleons (for LNV and LNC) and neutrinos (only for LNC). As a result of this, SN diffusion constraints are no longer competitive with those of COHERENT.
Regarding the U (1) B−L model, the bounds of Fig. 6 are also less restrictive than other existing limits. The different behaviour of these bounds with respect to the LNC scalar scenario, especially at low mediator masses, is due to the distinct scattering amplitudes, Eqs. (6) and (9), associated to each model. Note that in these cases the nature of the neutrino is the same (Dirac) and medium effects are introduced in an identical way.
For the LNV and the LNC scenarios, we also show in Figs. 4 and 5 the EoS limit given by Eq. (24), and the region (shaded in blue) where the conversion of ν L into ν R for the LNC scenario could distort the neutrino burst flux and time. As already mentioned above, the supernova EoS line is only indicative and shows the region above which LNV interactions could affect the EoS of SN matter. Likewise, the excluded region due to ν R production in the LNC scenario is based on qualitative arguments (the upper limit corresponds to demanding at least ∼ 100 new physics interactions in Eq. (26) to consider that ν R is trapped). This region has shifted upwards with respect to the results of Ref. [42] due to Pauli blocking and other medium effects, disfavouring a narrow band that was allowed in their work for Y ∼ 10 −4 .
Care must be taken when comparing our constraints with previous results in the literature [42,43], as there are substantial differences in the analysis. In particular, none of these previous works have included medium effects, i.e., Pauli blocking for the outgoing states and effective masses, and chemical potentials for nucleons. As we have shown in Fig. 2, these are responsible for a significant reduction in the scattering cross section (an increase in the mean-free path) which, in turn, results in a less stringent bound. The recent analysis of Ref. [43] also considered a radial dependence of the density in the proto-NS, however their input data corresponds to a snapshot at ∆t = 0.25 s after bounce (whereas we take ∆t = 1, 5 s for consistency) and it is based on a simulation that employs a different nuclear matter model (the SFHO [84]) than the one we take, where the neutrinosphere has a larger radius, R = 40 km. Besides, it must be noticed that their computation of the neutrino mean-free-path ignores the central region, where the nuclear density is higher, and, therefore, where neutrinos get to spend more time. This translates into weaker constraints since the region where more scatterings take place is not considered to calculate the time which neutrinos spend streaming out.
Finally, our choice of couplings for the scalar case (Y ) coincides with the notation used in Ref. [42], but the comparison with Ref. [43] must include a rescaling of their couplings by a factor √ Q , with Q ≈ 14 the nucleon coherence factor. For the U (1) B−L model, we can directly compare with Fig. 9 of Ref. [43]. Taking this into account, we notice that our bounds both for the U (1) B−L model and the scalar scenario are approximately a factor 2 more stringent than those of Ref. [43]. This mostly comes from the fact that they neglect the inner part of the star where more scatterings take place.

V. CONCLUSIONS
In this article, we have reevaluated the constraints on particle models with new low-mass scalar and vector mediators in the neutrino sector that can be derived from neutrino diffusion SN. To do this, we have computed the neutrino mean free path for three simplified scenarios, featuring either a new scalar (in a LNV and LNC model) or a new vector field (in a U (1) B−L construction), incorporating medium effects in the determination of the neutrino-nucleon scattering cross-section, and a radial dependence of the density, energy, and temperature inside the proto-NS.
The resulting diffusion time has been compared to the neutrino flux observed from SN1987A, which suggests that neutrinos do not remain trapped for longer than approximately 10 s. Using this as an upper bound, we have derived constraints on the properties of the new low-mass mediators (namely their coupling strength to neutrinos and their mass). We have compared these bounds with those from other experimental techniques.
Our results improve previous estimations, which did not take into account matter effects. In particular, we have found that matter effects lead to an increase of the neutrino mean free path that is more prominent in the central regions of the SN core where the temperature and density are larger. We have also shown that this effect is more important when new physics contributions dominate, for which the neutrino mean free path can increase by more than one order of magnitude with respect to its value in vacuum. This, in turn, leads to shorter diffusion times and relaxes the constraints on the neutrino coupling to nucleons.
The limits derived on the scalar mediator model are less restrictive than current bounds from other experimental sources, and the bounds on coherent elastic neutrino-nucleus scattering are leading for most of the values of the mediator mass. Likewise, in the U (1) B−L scenario, we have shown that the constraint lies comfortably within the area of the parameter space that has already been explored by neutrino experiments.
Our findings motivate dedicated numerical simulations of Majorana neutrinos in SN, which would more precisely incorporate the time and radial dependence of the stellar parameters in the computation of the neutrino diffusion time. This would also permit to study properly the effect of changes in the EoS. In line with previous works, we have indicated when new LNV interactions might alter the EoS, but we have been unable to check whether this has any effect on the SN observables.