Supernova Limits on Muonic Dark Forces

Proto-neutron stars formed during core-collapse supernovae are hot and dense environments that contain a sizable population of muons. If these interact with new long-lived particles with masses up to roughly 100 MeV, the latter can be produced and escape from the stellar plasma, causing an excessive energy loss constrained by observations of SN 1987A. In this article we calculate the emission of light dark fermions that are coupled to leptons via a new massive vector boson, and determine the resulting constraints on the general parameter space. We apply these limits to the gauged $L_\mu-L_\tau$ model with dark fermions, and show that the SN 1987A constraints exclude a significant portion of the parameter space targeted by future experiments. We also extend our analysis to generic effective four-fermion operators that couple dark fermions to muons, electrons, or neutrinos. We find that SN 1987A cooling probes a new-physics scale up to $\sim7$ TeV, which is an order of magnitude larger than current bounds from laboratory experiments.

For definiteness, let us consider a vector mediator Z ′ with mass m Z ′ and couplings to SM leptons and dark fermions χ with mass m χ , where ℓ = e, µ, τ , and g ℓ , g ν ℓ and g χ are generic couplings. Assuming that the dark fermions are sufficiently light (m χ ≲ 150 MeV), we will show in the following that their production from the PNS in SN 1987A leads to stringent constraints on their couplings to SM leptons. For the benchmark U (1) Lµ−Lτ model, this excludes large regions of the parameter space targeted by future experiments [21,[98][99][100][101][102][103][104][105][106][107][108]. While our analysis is valid for any mass of the Z ′ , we can integrate it out and describe its contribution with an effective four-fermion operator, if m Z ′ is much larger than arXiv:2307.03143v3 [hep-ph] 12 Oct 2023 the PNS's temperature and chemical potentials. This leads to significant simplifications in the analysis, and allows us to extend it to derive SN 1987A bounds on completely generic interactions of the dark sector with leptons through heavy portal mediators. In the heavy Z ′ limit, our calculations are in fact analogous to the ones necessary to study the SM production of neutrinos from leptons in stellar plasmas, which have received much attention over the past half century, since the seminal works in the 1960's [109][110][111][112][113][114][115][116][117][118][119][120][121][122][123][124][125] (see also Refs. [115,[126][127][128] for production of light dark fermions from heavy new physics). On the other hand, in case the Z ′ is a light and narrow state, the calculation is similar to the on-shell production of massive vector bosons coupled to leptons in the plasma [61,[129][130][131][132] (see also Refs. [35,59,61,133] for massive axions).
In this article we focus on the portal interactions with muons, but we also study neutrinos, which could be naturally linked to muons via SU (2) L . Electrons however form a highly degenerate and ultra-relativistic plasma in the PNS, which might lead to important medium effects in the electron and photon dispersion relations, requiring the inclusion of other production mechanisms not relevant for muons and neutrinos. With this caveat in mind, our calculations are easily extensible to electrons, generalizing and updating the pioneering work of Ref. [127] and improving the results presented in Ref. [128].
The rest of the paper is organized as follows. In Sec. II we outline the classical SN argument to constrain new exotic cooling agents using the neutrino flux observed from SN 1987A. Besides describing the general theoretical framework, we specify the SN simulations that we employ in our numerical analysis. In Sec. III we describe and compute the rates of the main emission mechanisms induced by the model in Eq. (1). We focus on extracting the main physical features of the rates using different approximations in the various regimes of the Z ′ mass, and on deriving analytical estimates. However, our final results rely on exact numerical computations whose details are deferred to Appendices. In Sec. IV, then, we implement these calculations of the rates in the SN simulations and derive the constraints on the parameter space of the Z ′ model in Eq. (1). We also generalize this analysis in terms of effective operators in the heavy Z ′ limit, and to one particular realization of the model arising from a L µ − L τ gauged symmetry. Finally, in Sec. V we summarize the results of our paper and close with a brief outlook.

II. SUPERNOVA COOLING
In the dense and hot environment within proto-neutron stars [134][135][136][137] neutrinos become trapped and a thermal population of muons is predicted to arise [138,139]. New light dark particles that couple to leptons, e.g. via the interactions in Eq. (1), can be produced efficiently in the stellar plasma, leading to a significant loss of en-ergy if they can escape from the PNS. The corresponding dark luminosity L χ is then subject to the classical bound L χ ≲ L ν at 1 s post-bounce, where L ν is the neutrino luminosity [30,40,41]. This limit is obtained from the observation of a neutrino pulse over ∼ 10 s [140,141] during SN 1987A [142][143][144], which is in accordance with the predictions of the standard theory of core-collapse SN (see Refs. [59,145,146] for a critical reappraisal of this limit 1 ).
Here, we apply this argument to scenarios where light dark fermions couple to leptons with interactions such as those in Eq. (1). One needs to distinguish two regimes based on the mean free path (MFP) of the dark fermions in the plasma or, equivalently, the strength of the portal interactions. If the dark particles are very weakly coupled (or the MFP is much larger than the radius of the PNS) then they free stream out from the SN once produced, whereas for large couplings (or MFP much shorter than the radius of the PNS) they thermalize with the medium and get trapped inside of the PNS.
In the free-streaming regime the general expression for the total energy-loss rate per unit volume, Q, for a given emission process is These are thermal integrals over the phase space of all the initial-and final-state particles weighted by their number density distributions f i and the Pauli blocking or Bose enhancement factors (1 ∓ f j ), respectively. Furthermore, |M| 2 is the the squared matrix element of the given production process and E χ is the total energy carried away by the dark particles. In the calculations of the freestreaming regime, one conventionally uses f χ = 0 for the new particles, because their occupation numbers inside the PNS are very low and not thermalized by assumption.
In the deep trapping regime, on the other hand, the dark-sector particles are in thermal equilibrium with the plasma and they are emitted from a surface with radius r χ (dark sphere) following a law analogous to the one of the black body radiation, where g χ is the number of degrees of freedom of the χ particle (g χ = 2 for massive dark fermions), x m = m χ /T χ and T χ = T (r χ ) is the temperature of the dark sphere. The radius r χ is defined, as is conventional in astrophysics [149][150][151], through the optical depth τ χ (r), by requiring where λ(r) is a suitable spectral average of the dark fermion's MFP at a radius r. In this work, we use a "naive" thermal average for computational simplicity 2 (see Appendix A for our definitions of thermal averages). The energy-dependent MFP λ(r, p χ ) is related to the total rate of interaction of a dark-sector particle in the medium, The contribution to Γ χ of a given process with a bunch of target particles b colliding with χ in the initial state is defined through The quantity C abs is the collision operator describing the absorption rate per unit volume of the medium which uses the same definitions as in Eq. (2), except for |M| 2 which is now the squared matrix element of the given absorption process 3 . For the numerical analyses of this paper we use SN simulations including muons presented in Ref. [59] and whose radial profiles for the relevant quantities are reported in [152]. Our fiducial results are obtained using the simulation labeled as SFHo-18.80, which reaches the 2 We have checked that other averages, such as the conventional Rosseland MFP [149][150][151], give very similar results. In the recent literature also other forms of taking the spectral average have been discussed, for example by including an energy-dependent opacity in the free-streaming spectral luminosity, integrating over the energy only in the very end [61]. For a discussion and comparison of the various approaches see Ref. [35]. 3 Denoting Γours for our definition, the standard approach in the literature [35] is to work in terms of emission rates Γ E = gχfχΓours. However, one then has to use Boltzmann-equation arguments to argue which rate has to be used for absorptive processes, leading to the reduced absorption rate Γ A . The definition Γours is chosen such that Γours = Γ b χ = n b ⟨σv⟩ b , which naturally leads to Γours = Γ A (see Appendix A). lowest temperatures and, therefore, will lead to the most conservative limits on the dark luminosity (at 1 s postbounce). The upper bound is set by the neutrino luminosity calculated within the same simulation, which for SFHo-18.80 is given by 4 L χ ≤ L ν = 5.7 × 10 52 erg s −1 .
For a rough estimate of the systematic uncertainties related to SN modeling, we will also show the more stringent limits obtained from using the hotter SFHo-20.0 simulation, which gives L χ ≤ 1.0 × 10 53 erg s −1 . In the free-streaming regime, the dark luminosity is obtained as a volume integral of Eq. (2), L χ = QdV , while in the trapping regime we use Eq. (3). Finally, it will be useful to estimate the contributions to the dark luminosity of the different processes to understand their relative importance. For this, we define "typical PNS conditions" as those at 1 s post-bounce and at a radius ≈ 10 km. This region dominates the volume emission in L χ and is representative of the bounds in the free-streaming regime. Using the simulation SFHo-18. 80 [152], this approximately corresponds to: Here T denotes the temperature, ρ the density, µ l the chemical potential of the lepton l and Y ℓ is the number density fraction of the charged lepton ℓ relative to the one of baryons. For the Y ℓ we quote the results derived from the rounded temperature and chemical potentials in Eq. (9) and, therefore, they are slightly different to those reported in [152]. Let us stress again that Eq. (9) will be only used for numerical estimates, while our final results and constraints on the models will be obtained using the full radial profiles of all relevant thermodynamical quantities.
Note also that this is not generally true for production from electrons because they are ultra-relativistic and form a highly degenerate system that suppresses photoproduction compared to bremsstrahlung and annihilation [30]. Moreover, there are important plasma effects which, for example, dress the electron with an effective mass m * e ∼ 10 MeV and give rise to pseudo-particle excitations that need to be taken into account in a realistic analysis (see e.g. Ref. [133] for the emission of massive axions from electrons in SN). Nonetheless, the calculations we present in this work can be easily extended to electrons and compared with previous literature where all these effects have been neglected [127,128]. We will estimate some of them below in Sec. III C.
In case of absorption, there are the inverse processes µ − χχ → γµ − and χχ → µ − µ + , whose rates are related by detailed balance to those of the photoproduction and annihilation production, respectively, provided that the χ andχ particles reach thermal equilibrium. In addition, other scattering processes may contribute to the diffusion and energy transport in the trapping regime, such as χµ − → χµ − and χν ℓ → χν ℓ (see middle panel of Fig. 1), or processes in the dark sector such as χχ → χχ.
In Appendices B and C we provide the cross sections for all relevant 2 → 2 and 2 → 3 processes needed to calculate the energy-loss and absorption rates. In the following, we discuss in detail the contributions of the annihilation and photoproduction topologies.

A. Annihilation
The energy-loss rate per unit volume in Eq. (2) for µ − µ + → χχ annihilation can be simplified to In this equation, g µ are the muon's spin degrees of freedom, E −(+) denote the muon (anti-muon) energy in the PNS rest frame, p ± are the absolute values of their 3momenta, p ± ≡ E 2 ± − m 2 µ , and f ± ≡ f ± (E ± ) = 1/(e (E±±µ)/T + 1) are their Fermi-Dirac distributions in the medium. The function I s is a (dimensionless) angular integral over the annihilation cross section σ(s) = σ(µ + µ − → χχ) which depends on the colliding angle θ between the muon and anti-muon through the Mandelstam variable Top-to-bottom panels represent annihilation, scattering and photoproduction processes, respectively. Annihilation and scattering diagrams will also contribute when ℓ is replaced by a neutrino.
s in the PNS frame, s = 2 m 2 µ + E + E − − p + p − cos θ . The cross section is physical only above the 2-particle threshold which imposes the kinematic constraint s ≥ 4 max(m 2 χ , m 2 µ ) in the angular integral. The annihilation cross section in Eq. (11) for the Z ′ model in Eq. (1) is where we have introduced β i (s) = 1 − 4m 2 i /s, κ i (s) = 1 + 2m 2 i /s and the total Z ′ width and where θ(x) is the Heavyside step function. Note that the average energy carried away by the dark fermions in the annihilation process is equal to the thermally averaged center-of-mass (CM) energy of the leptons, For typical conditions in the PNS, Eq. (9), E µ ∼ 280 MeV, which sets the scale for 2m χ above which the production of χ's in the plasma will become exponentially ("Boltzmann") suppressed by the distribution functions f ± . In addition, this energy scale allows one to define three regimes of m Z ′ in Eq. (12) depending on which term dominates the denominator: (i) A "heavy regime" in which m Z ′ ≳ 1 GeV ≫ E µ , so that the Z ′ is too heavy to be produced on-shell; (ii) the "resonant regime" where the Z ′ can be produced on-shell, 2m µ ≤ m Z ′ ≲ E µ ; and (iii) the "light regime" with Z ′ masses below the twomuon threshold, m Z ′ < 2m µ , so that the Z ′ is too light to be produced on-shell. Analogous expressions can be defined for electrons and neutrinos by replacing the couplings and masses accordingly (notice that for neutrinos g ν ℓ = 1). Moreover, the average CM energies in these cases are E e ∼ 160 MeV and E ν ∼ 130 MeV, and analogous regimes to those for the muons can be formulated for neutrino-antineutrino and electron-positron annihilation. The light regime in these cases is restricted to extremely small Z ′ masses, making it irrelevant for the range of vector boson masses we are considering here.
The demarcation of these regimes is useful because one can use approximations to derive analytic results and isolate the main physical factors in control. In the following, we discuss these approximations and describe their contributions to the absorption rate Γ χ .

Heavy regime
In this case, the denominator in the propagator (see Eq. (12)) is dominated by the Z ′ mass. Expanding in powers of s/m 2 Z ′ up to leading order, the cross section can be easily integrated analytically, giving a function I s (E + , E − , m χ , m µ ) proportional to the effective coupling g 2 χ g 2 µ /m 4 Z ′ . We can further approximate this expression by taking the high-energy limit m χ → 0 and m µ → 0, obtaining 5 Also setting m µ → 0 in the integrals in Eq. (10), the integrations can be carried out analytically, giving Here we have re-scaled the chemical potential y = µ/T and introduced the functions where Li n+1 (z) is the polylogarithm of order n + 1. If we also take vanishing chemical potentials, we recover the results in Ref. [127] Q heavy 5 In the following discussion we fix gµ = 2 with the understanding that some intermediate formulas change by factors of 2 for the neutrino case.
in terms of the Riemann ζ-function with F 4 F 3 ≈ 133. We have also included a subindex in Q to indicate that this is a zeroth-order approximation neglecting masses and chemical potentials of the leptons. In order to assess the accuracy of the above approximations, we compare Q in Eq. (10) for massless dark fermions and the cross section in the heavy Z ′ limit with Eq. (18) for different SM particles at the typical conditions of PNS in Eq. (9). For muons one finds Q µ /Q 0 ≈ 0.33 while for neutrinos and electrons one finds Q ν /Q 0 ≈ 0.99 and Q e /Q 0 ≈ 0.54 (using the physical electron mass in vacuum), respectively.
The thermal suppression of the muon population is mild for these conditions in the PNS, Y e ≃ 4Y µ . The positron abundance is also suppressed by the large electronic chemical potential and, hence, for the same couplings to electrons and muons one obtains similar rates.
With these approximations one can estimate the parametric dependence of the energy loss rate per unit mass (i.e. the emissivity) produced by lepton annihilation in the heavy regime as where we have divided Eq. (18) by the density in Eq. (9), and ϵ max = 2.1 × 10 19 erg s −1 g −1 has been estimated dividing L χ in (8) by the total mass of the PNS in this simulation M PNS = 1.351 M ⊙ .

Resonant regime
If the Z ′ can be produced on-shell, then the denominator in Eq. (12) is dominated by the Z ′ width Γ Z ′ , and it can be replaced by π/(m Z ′ Γ Z ′ )δ(s − m 2 Z ′ ) in the narrow width approximation. The δ-function can be used to perform the angular integration in Eq. (11) and for m χ → 0 this gives (neglecting terms of relative size 2m 2 µ /m 2 Z ′ ) The energy integrations of Eq. (10) can be wellapproximated by neglecting the chemical potentials, but keeping a non-zero muon mass in the integrand, giving where BR χ ≡ BR(Z ′ → χχ) denotes the invisible Z ′ branching ratio, and numerically m µ /T e −2mµ/T ≈ 0.004. For muon annihilation this indeed yields a good approximation, with Q µ /Q µ=0 ≈ 0.95, while for electrons one can set m e → 0 in the integrals, giving a result similar to Eq. (16), which can be further approximated by 6 and numerically µ 2 e /4T 2 e −µe/T ≈ 0.06, resulting in Q e /Q m=0 ≈ 1.6. Finally neglecting lepton masses and chemical potential simplifies to where numerically F 1 F 0 ≈ 0.57. This is a good approximation for neutrino annihilation with Q ν /Q 0 ≈ 0.95, while energy-loss rates for electrons (muons) are smaller by a factor 10 (100). Importantly, the contribution to the energy-loss rates in the resonant regime scales perturbatively with the couplings as ∼ O(g 2 ) instead of ∼ O(g 4 ) in the heavy or light regimes. In fact, for BR χ = 1, one should recover the results obtained for the coalescence Z ′ production mechanisms in Ref. [61] 7 . Also, the various rates in Eqs. (22), (23) and (24) all scale quadratically with the Z ′ mass.
From these expressions one can readily obtain the emissivities, which for neutrino annihilation read where we have used the same approximations as in Sec III A 1 which are valid up to m Z ′ ∼ 200 MeV. As discussed above, emissivities for electrons and muons are expected to be smaller by a factor 10 and 100, respectively.

Light regime
In this case the denominator of the propagator is dominated by s, and the cross section can be integrated analytically. For massless χ and muons one obtains which is independent of E + , E − . Similarly as in the resonant regime, the energy integrations of Eq. (10) can be approximated by neglecting the chemical potentials, but keeping a non-zero lepton mass in the integrand where m 3 µ /T 3 e −2mµ/T ≈ 0.04. For muon annihilation this indeed yields a good approximation, with Q µ /Q µ=0 ≈ 0.85, while for electrons one can set m e → 0 in the integrals (but keep non-zero chemical potentials), giving a result similar to Eq. (23) and numerically µ 3 e /6T 3 e −µe/T ≈ 0.2, resulting in Q e /Q m=0 ≈ 1.2. Finally neglecting lepton masses and chemical potential simplifies to where numerically F 2 F 1 ≈ 1.5. This is a good approximation for neutrino annihilation with Q ν /Q 0 ≈ 0.97, while for electrons (muons) the energy-loss rates and emissivities are smaller by a factor 10 (100) than predicted by this formula. Nevertheless, using Eq. (29) we obtain the emissivity

Annihilation and scattering contributions to trapping
Given the scattering of a χ with another particle b, the absorption rate can be approximated by (see Appendix A) where g b are the number of degrees of freedom of the particle b (g b = 2 for b =χ, µ, e and g b = 1 for b = is the Møller velocity and σ(s) is the scattering cross section for χ + b → X 1 + . . . + X n [153]. Moreover, the index i runs over the finalstate particles and we have approximated the effect of the Pauli blocking by its thermal average or degeneracy factors [30] where g i and n i denote their degrees of freedom and number densities, respectively. There are two types of processes related by crossing to the annihilation diagram that are relevant for the trapping regime: Inverse annihilation, χ χ → µ − µ + and scattering, χ µ − → χ µ − andχ µ − →χ µ − . Scattering processes are kinematically more involved as they exchange a Z ′ in the t-channel. For very light Z ′ 's (m Z ′ ≪ T ), they have a differential cross section with a Coulombian enhancement in the forward direction which involves a small momentum transfer and, therefore, little contribution to the thermalization rate between the dark and SM sectors.
The interplay between the contributions of inverse annihilation and scattering to the absorption rate is similar to the case of heavy-lepton neutrinos in SN [154], as recently emphasized in Ref. [155]. In the absence of selfinteractions between the dark fermions, these two processes really define two surfaces that determine different properties of the dark luminosity in the trapping regime. The freeze-out of inverse annihilation first fixes the number flux of χ's. The outgoing flow then thermalizes via scattering processes with the leptons until they decouple at a larger radius. For Z ′ masses in the resonant regime, absorption rates will be dominated by the inverse annihilation and the two surfaces coalesce into a single dark sphere that determines the thermal emission of the χ's. In the heavy regime both mechanisms can be important and this distinction must be kept in mind. This is also reminiscent of models with new neutrino self-interactions and their effect in the dynamical evolution of the SN [156][157][158][159][160][161][162]. Notably, when considering the model in Eq. (1) with neutrino interactions, processes such as ν ℓνℓ → ν ℓνℓ also occur. These effects may lead to a fundamentally different incarnation of the SN 1987A cooling limit, valid in the trapping regime, which is however beyond the scope of the classical SN 1987A cooling bound that we apply in our analysis.
With all this in mind, our fiducial analysis in the trapping regime includes both processes, inverse annihilation and scattering, in the calculation of the MFP to obtain one single dark sphere that determines L χ . However, we repeat the calculations for the case where we do not include the scattering processes. We take the variation of our results as an indicator of the potential systematic uncertainties involved in our treatment of the trapping regime.
In addition, dark elastic scattering, χχ → χχ, may become relevant in the trapping regime, as recently discussed in Refs. [46,163]. In our case, however, this will not play an important role for the calculations of L χ in the trapping regime. The reason is that dark elastic scattering does not directly contribute to maintaining the population of χ's in thermal equilibrium with the SM plasma after they freeze out (after inverse annihilation turns off). Moreover, for low m Z ′ , the rate of dark elastic scattering is resonant and overwhelmingly larger than scattering or inverse photoproduction (discussed below in Sec. III B 3), that do tend to maintain the thermal equilibrium between the dark and SM sectors. This situation is particularly relevant for muons, whose population drops significantly at the outer layers of the PNS where freeze out of the χ's occurs. Therefore, they would effectively decouple immediately after, except for the contributions to the MFP induced by their interactions with the leptons, which are already accounted for in our fiducial analysis 8 . Nevertheless, we have studied the contributions of dark elastic scattering processes for completeness, and discuss in Sec. IV the consequences for our results if these are included at face value in the calculation of Γ χ and the MFP.

B. Photoproduction
For muons we use two approximations: (1) Describe the Pauli blocking of the final muon by ( (32) and, (2) in the phase space integrals of the initial particles we take the extreme non-relativistic limit where the muon is static and recoilless. Thus, the kinematics is evaluated in the muon's rest frame, s = m 2 µ +2m µ ω, with ω the photon's energy, and E χ + Eχ = ω. One arrives at the conventional formula [30] where ω 0 = (s 0 − m 2 µ )/2m µ and s 0 is the kinematic threshold of the process.
While Eq. (33) is approximately valid for muons, electrons are instead ultra-relativistic in the PNS. If we neglect the electron's mass in the phase space integrals 9 and assume that E χ + Eχ ≈ (ω + ω ′ )/2, then where ω ′ is the electron energy and s = 2ωω ′ (1 − cos θ).
Similarly to the annihilation topologies, one can estimate the average energy available for the production of dark particles in the photoproduction process. For this we use the average CM energy at threshold, and one can also define different regimes of m Z ′ for photoproduction: where the Z ′ can be produced on-shell. For the conditions specified in Eq. (9) we find E γ µ ∼ 90 MeV and E γ e ∼ 150 MeV. In Appendix C we provide the results for the photoproduction cross sections along with some details of the calculations. In addition, we have checked the validity of the approximations in Eqs. (33) and (34) by calculating the full 5-body phase space thermal integrals exactly, as described in the seminal work of Ref. [112], where this was done for effective (axial)vector operators. For completeness, we outline this calculation in Appendix D. We find that the approximate formulas provide a very accurate description, within ∼20% (muons) and ∼ 40% (electrons) of the full photoproduction rates in the resonant regime and neglecting degeneracy of the final state lepton. In the heavy regime, Eq. (33) underestimates the full result for muons by a factor ∼ 8. Nevertheless, as shown below in Sec. III B 1, the emission rate for heavy m Z ′ is dominated by annihilation by orders of magnitude. On the other hand, for the highly degenerate electrons in the PNS, we find that approximating the Pauli blocking effects by F deg,e overestimates the full energy loss rate by a factor ∼ 3. In the case of muons, the degeneracy factor F deg,µ instead agree with the exact result within ∼ 5% accuracy.
In the following, we discuss in more detail the two regimes of photoproduction and the contribution of inverse photoproduction to the interaction rate Γ χ .

Heavy regime
The cross section of the process µ − γ → µ − χχ, induced by the exchange of a heavy Z ′ , is × − 55ŝ 6 + 682ŝ 5 + 483ŝ 4 − 968ŝ 3 − 169ŝ 2 + 30ŝ − 3 + 12ŝ 2 2ŝ 4 − 14ŝ 3 − 87ŝ 2 − 52ŝ + 1 logŝ , withŝ = s/m 2 µ , α denotes the fine-structure constant and we have taken for simplicity the massless limit m χ → 0. This cross section is equivalent to the expression obtained by Dicus in [114] for the photoproduction of a neutrino pair using only vectorial couplings 10 . One might attempt to obtain a more simplified expression of the rate by taking the non-relativistic limit in Eq. (36). However, the cross section in this case, grows rapidly with energy and leads to a gross overestimation of the integral in Eq. (33) (see also Ref. [112]). For instance, for the typical PNS conditions used above, one obtains a rate that is larger than the one obtained using the relativistic expression of the cross section by a factor ∼ 25. Nevertheless we use the non-relativistic approximation together with Eq. (33) to get a rough estimate of the emissivity produced by the photoproduction process in this regime, giving Comparing with the equivalent contribution from annihilation in Eq. (20), we conclude that photoproduction gives an emissivity rate that is smaller by a factor ∼ 250 and thus can be neglected in the heavy regime.

Resonant regime
where we have introduced the notation α χ = g 2 µ /4π, and whereŝ = s/m 2 µ , x ′ = m 2 Z ′ /m 2 µ and R(x,ŝ) = (ŝ 2 − 2ŝ(x+1)+(x−1) 2 ) 1/2 . For BR χ → 1 we recover the cross section for semi-Compton production of massive vector bosons, γµ − → Z ′ µ − . In the m Z ′ → 0 limit, one obtains and if we now perform the non-relativistic expansion, we recover the Thomson cross section for α χ → α and BR χ → 1. This is the expression commonly used for the semi-Compton production of vector particles in stellar plasmas [30,61], but it is less appropriate for the PNS where leptons are relativistic. In fact, in case of muons for m Z ′ = 0 and in the typical conditions we have been using for the PNS, we find that the Thomson cross section overestimates the relativistic one by a factor ∼ 2.
On the other hand, the energy-loss rate of the full resonant photoproduction cross section is insensitive to m Z ′ up to m Z ′ ≳ T , at which point it starts dropping due to increased Boltzmann suppression and defines the onset of the heavy regime. Taking as a reference the Thomson cross section, we can estimate the emissivity of the photoproduction in the resonant regime as using our typical PNS conditions in Eq. (9). Comparing this to the emissivity from µ + µ − annihilation, Eq. (30), we observe that the rate of χχ production from muons for light Z ′ will be dominated by photoproduction for many orders of magnitude. For m Z ′ ≲ 10 MeV this process is even more important than resonant neutrinoantineutrino annihilation for the case g ν ℓ = g µ , as demonstrated by comparing to Eq. (25).

Contribution of inverse photoproduction to trapping
Assuming thermal equilibrium, which is adequate in the trapping regime, the contribution of inverse photoproduction ℓ − χχ → ℓ − γ to Γ χ can be related to the production rates by means of detailed balance (see Appendix A). In case of muons, the photoproduction rate of χ per unit volume can be calculated using the same approximations as in Eq. (33), while for electrons we use Then, we estimate the contribution of photoproduction to the MFP by applying detailed balance, using the inverse of ⟨Γ γ χ ⟩ (see Appendix A) Note that this approximation differs from the direct calculation of the contribution of inverse annihilation and scattering in Eq. (31). In order to combine the two contributions in our estimate of the MFP we use the approximate formula where all the thermal averages are understood to be taken with respect to the χ kinematics (see Appendix A).

C. Other processes and neglected plasma effects
In our analysis, we have selected the processes that are dominant for the production and absorption of χ's in muons and have neglected plasma effects which are expected to be small in this case. Electrons in the PNS, on the other hand, are highly degenerate. Moreover, in the plasma the electron and photon dispersion relations are significantly modified. The electron mass effectively increases while the photon acquires a longitudinal mode and an effective mass that could enable the decaỹ γ → χχ, where the "plasmon"γ includes these collective plasma modes [30,111,119,121].
Nonetheless, the production of χ's in the heavy regime and neglecting the χ mass would be analogous to the SM pair-production of heavy-lepton neutrinos from the electrons in the stellar plasma. Plasmon decay is indeed an important process in the conditions of high densities predicted in the PNS. However, at the high temperatures reached in the SN explosions considered in this work, e + e − annihilation becomes the dominant process [121]. Adding mass to the χ's will not affect this conclusion and may, in fact, kinematically close the plasmon decay if m χ ≳ 10 MeV, which is the scale of the plasma frequencies expected in the PNS. Finally, accounting for the increase of the electron mass in the plasma by a similar amount does not affect significantly the annihilation rates as discussed in Sec. III A 1.
On the other hand, for the light Z ′ some of the neglected effects can become important. For instance, resonant bremsstrahlung production e − p → e − p(Z ′ → χχ) could become the dominant process for m Z ′ ∼ 10 MeV, as it is the case for on-shell production of axion-like particles [51]. In addition, for these masses one needs to consider medium-inducedγ − Z ′ mixing which may also have an impact [165]. For all of these reasons, our results for the emission from electrons in the light regime, m Z ′ ≲ 50 MeV, should be considered as an intermediate step towards a more refined and robust calculation, and the results that we report for this case regarded as a rough approximation.

IV. RESULTS
In this section, we apply the upper limit on the dark luminosity in Eq. (8) to obtain the SN 1987A constraint on the parameter space of the different dark sector models we consider. We start by presenting an analysis in terms of effective operators, which corresponds to the heavy regime introduced in Sec. II for the calculation of the relevant processes. This allows us to generalize our analysis to generic interactions of dark fermions coupled to leptons via dimension-6 operators, and extract a SN 1987A limit for any portal mediator with mass much larger than the temperatures and chemical potentials in the PNS. It is interesting to note that this approach was already applied in the context of neutrino emission from stellar plas- mas in the very early days of the SM [115,126]. We then focus on the constraints on the parameter space of the simplified Z ′ model of Eq. (1), which can be regarded as the continuous extension of the EFT bounds to low mediator masses for one particular operator (V V ). Finally, we present the SN 1987A constraints for a phenomenologically relevant and UV-motivated version of the Z ′ model, obtained by gauging the L µ − L τ symmetry.

A. Effective field theory
The analysis for the heavy Z ′ can be generalized in the context of an EFT with four-fermion operators. Focusing on the couplings of dark-sector fermions χ to leptons l = e, µ, ν ℓ , the most general effective Lagrangian at leading order is where X, Y run over V, A, S, P, L, R, T, T ′ , with Γ V = γ µ , and Γ T ′ = σ µν γ 5 , and Lorentz indices properly contracted 11 . Matching the effective operators to the Z ′ model in Eq. (1) coupled to charged leptons yields Λ ℓ = m Z ′ and C ℓ V V = g ℓ g χ . For neutrinos their bilinears in the EFT Lagrangian are constructed with left-handed fields and contribute, instead, to C ν ℓ LV = g ν ℓ g χ . 11 Only two tensor operators are independent and we will take those corresponding to C l T T and C l T ′ T .  47), with mχ = 0. The orange bars show the exclusion range obtained from SN 1987A, using the simulation SFHo-18.8. The upper limit is obtained by the free-streaming regime and the lower limit is set by the trapping regime. The red bars show the lower limits obtained at LEP [166], while the purple and green hatched bars show the projections for Belle-II with 50 ab −1 of integrated luminosity [167] and ILC with √ s = 1 TeV and Lint = 1000 fb −1 [168].

SS V V, AA
In Table I, we show the limits obtained on Λ eff l ≡ (Λ l /C l XY ) 1/2 for these interactions in the limit m χ = 0 for muons, neutrinos and electrons. The upper limits correspond to the free-streaming regime, the lower limits to the trapping regime. Notice that the excluded regions of Λ eff l are in the EFT range of validity as long as Λ l ≳ 1 GeV, which is much larger than all other energy scales relevant in the PNS. This is the case for essentially all operators.
The SN 1987A bounds on the EFT operators are very strong, reaching up to 4 − 7 TeV. This sensitivity to the mediator mass scale is approximately one order of magnitude better than the one achieved by laboratory experiments for similar leptonic interactions [166,167,169,170] 12 . For instance, monophoton searches at LEP have been used to set the following lower bounds on Λ eff e [166], In Fig. 2, we show the SN 1987A limits on these EFT operators compared with those obtained at LEP [166], and with the projections of the sensitivity that could be achieved at Belle II [167] or at a future e + e − linear collider [168]. Remarkably, the SN 1987A bounds are stronger than LEP by roughly one order of magnitude, and will even dominate over future collider limits. Note however that SN 1987A constraints apply only to sufficiently light dark fermions, while collider bounds typically extend to larger χ masses, such as LEP, which provides constraints for m χ ≲ 100 GeV.
As discussed in Sec. II, the dominant effect in the heavy regime of the Z ′ model is annihilation (for typical SN conditions), and photoproduction can be neglected. We have checked that this is indeed true for any EFT operator by calculating their contributions to photoproduction explicitly (see Appendix C). In fact, we note that, for the EFT limit of the Z ′ model, the approximate expression in the high-energy limit (m ℓ → 0) in Eq. (18) leads to a bound on the effective scale that is off only by ∼ 30% with respect to the full calculation. Therefore, since the high-energy limit gives quite accurate results, the numerical bounds for different Lorentz structures and different leptons are of the same order, as in this limit the corresponding cross sections differ at most by O(1) numbers.
In Fig. 3 we show the dependence of the limits on the heavy scale as a function of the dark fermion mass m χ for the effective V V interaction and LV for neutrinos. As expected, the excluded regions shrink with increasing m χ and they are limited to masses m χ ≲ 300 MeV. We also analyze the variations of the constrained region produced by using the simulation SFHo-20.0 or different prescriptions for the processes included in the trapping regime. The upper limits of Λ eff l obtained from this hotter simulation are a factor ∼ 2 stronger in the freestreaming regime, because in this case emission is dominated by the hottest region inside the supernova. On the other hand, the uncertainties of the boundary with the trapping regime for electrons and neutrinos are relatively small and our fiducial calculation is, again, on the conservative side. The reason for this behavior is that the dark sphere is not located in the region of highest temperature. This makes trapping sensitive to the shape of the temperature profile, which is similar for both simulations.
Moreover, the scale marking the onset of the trapping regime in these cases is Λ eff l ≈ 100 GeV, which is of the order of the electroweak scale. This is consistent with the fact that the boundary of the trapping regime is set by the luminosity of the trapped neutrinos via Eq. (8), which interact with SM leptons precisely through dimension-6 operators suppressed by the Fermi scale.
For muons, however, the bound extends to scales that are a few orders of magnitude lower than for electrons and neutrinos. In addition, there is a large variation in the location of the boundary depending on the selection of processes contributing to trapping. The reason is that muons are relatively heavy and there is a maximal radius where they can be produced by thermal fluctuations in the plasma. Putting it differently, there are no muons to scatter with in case of inverse photoproduction and inverse annihilation becomes ineffective because of the strong phase-space suppression. In case that only µ − χ interactions are included in the calculation of the MFP, the radius of the dark sphere is typically smaller than the one of the neutrino sphere and the bounds on Λ eff µ extend down to ∼ 1 GeV. However, if one were to also include χχ elastic scattering in the calculation of the MFP, and the χ is light, m χ ≲ T , then they would produce contributions analogous in size to those given by the neutrino and electron interactions and the boundary of the trapping regime would be increased again to ∼ 100 GeV,  cf. the discussion in Section III A 4.

B. Simplified Z ′ models
We now present the SN 1987A constraints in the parameter space of the simplified Z ′ model in Eq. (1). In the derivation of our results for the annihilation contributions to Q we numerically solve the integral in Eq. (10) (and the equivalent ones for neutrinos and electrons), as described in Appendix E. This procedure allows us to track the contribution of the annihilation rates across the whole m Z ′ range, including the transition regions between the three regimes defined in Sec. III A. Instead, for photoproduction we only use the expressions in the resonant regime, because photoproduction is only relevant for the low-m Z ′ regime of the charged leptons, see Eq. (38).
In Fig. 4 we show the SN 1987A limits for the model in Eq. (1), for the case when only one of the leptonic couplings is present and with g χ = g l . In the lower (upper) part of the plots we identify the dominant process for production (absorption) in the indicated mass range. We also show the variation of the constrained region obtained by using the SFHo-20.0 simulation (dotted curve) and, independently, by omitting the scattering contribution to the MFP in the trapping regime (hatched region). In all cases we observe the onset at high masses of the power-law behavior of the constraints ∝ 1/m 4 Z ′ , characteristic of the EFT. Interestingly, this occurs for m Z ′ ≳ 1 GeV in the free-streaming domain but already for m Z ′ ≳ 100 MeV in the trapping regime. This is because in the former case emission is governed by the conditions in the hottest region of the PNS, while in the latter it corresponds to the cooler outmost layer of the dark sphere, where all relevant energy scales are smaller. The boundary of the trapping regime changes very little with respect to adding or not scattering contributions in the absorption rate. Adding dark elastic scattering (χχ → χχ) in the muonic case (for m χ ≲ T ) would instead make the trapping region similar to the neutrino case.
The shape of the constrained region below m Z ′ ≈ 1 GeV depends on the lepton considered. For muons the low m Z ′ region is dominated by (resonant) photoproduction, which gives a flat bound up to m Z ′ ≳ T , where the on-shell production of the Z ′ starts decreasing due to Boltzmann suppression. Nevertheless, it remains more important than µ + µ − -annihilation (in the light regime), until it becomes resonant at the µ + µ − threshold, m Z ′ ≥ 2m µ . However, this occurs already at large energies and the production quickly suffers from the Boltzmannian suppression, converging to the EFT scaling at higher m Z ′ .
For electrons, photoproduction dominates again for low m Z ′ , even though the e + e − threshold is much lower than for muons. Resonant annihilation of electrons quickly replaces photoproduction as the dominant process above m Z ′ ≳ T in the free-streaming domain, until the resonance starts suffering from Boltzmann suppression and the EFT takes over. For the trapping regime, there is no range of m Z ′ where resonant annihilation dominates and the EFT directly replaces inversephotoproduction at m Z ′ ≳ 100 MeV.
For neutrinos there is no photoproduction and the production and absorption rates are given by annihilation in the resonant and heavy regimes. In the free-streaming regime we observe a strengthening of the bound up to T ≲ 100 MeV. This is due to the m 2 Z ′ scaling of the emission rate in the resonant regime, see Eq. (24), which is quickly overcome by Boltzmann suppression until the heavy regime takes over.

C. Gauged Lµ − Lτ models
Finally, we study SN 1987A constraints on a UVmotivated realization of the simplified Z ′ model in Eq. (1). This is the gauged L µ − L τ model coupled to DM fermions, described by the interaction Lagrangian where χ is the dark fermion and j µ SM is the SM part of the L µ − L τ current. These interactions induce an irreducible contribution to the kinetic mixing of the Z ′ with the photon, through µ and τ loops, giving a mixing parameter of the order ϵ ∼ g µ−τ /70 [172]. This would presumably give only small corrections to our analysis (see e.g. [61]) and thus will be neglected.
It is well known that such models can accommodate the present (g − 2) µ anomaly with couplings of order g µ−τ ∼ 10 −4 for a light Z ′ gauge boson, m Z ′ ≪ m µ [173][174][175]. They also allow to reproduce the DM relic abundance through resonant s-channel annihilation, when muon and DM fermion have similar couplings to the Z ′ gauge boson and the latter is heavier than the DM fermion by a factor 2 − 3 [77,78]. For such light masses, Z ′ decays and DM annihilation can heat the SM bath after neutrino decoupling, thereby increasing the effective number of relativistic degrees of freedom, usually expressed as an effective number of neutrino species N eff . Such a contribution could help to reduce the longstanding tension between local and cosmological determinations of the Hubble constant [95,97], if the new contributions is of order ∆N eff ∼ 0.1 − 0.4 [1].
At present, the most relevant laboratory constraints (see Ref. [70] for an overview) stem from BaBar searches for Z ′ -bosons above 212 MeV decaying to muons [176], neutrino trident production [177] at CCFR [178], bounds on coherent elastic neutrino nucleus scattering from the COHERENT collaboration [179] and constraints on neutrino-electron scattering [180,181] at BOREX-INO [182]. A variety of accelerator-based searches have been proposed to explore the unconstrained parameter space, such as NA62 [98], which looks for final state radiation of Z ′ -bosons in K + → µ + ν µ , and dedicated searches using muon beam facilities such as the NA64µ experiment [183] at CERN and M 3 [99] at Fermilab.
Astrophysical limits from white dwarfs have been studied already in Refs. [70,184]. Here we show that also constraints from SN 1987A disfavor a large region of parameter space with significant overlap with the expected reach of planned experiments and with the region that could address the H 0 tension. In Fig. 5, we show the SN 1987A limits on the muon coupling as a function of the Z ′ mass, in the scenario where muons and dark matter couple equally to the Z ′ , g µ−τ = g χ , and m Z ′ = 3m χ . This is compared to the current bounds discussed above, shown in gray, and the regions preferred at 95% C.L. by (g − 2) µ and H 0 .
In the free-streaming regime we use the same method as for the simplified model described in Appendix E, including the contributions from µ + µ − , ν µ ν µ and ν τ ν τ an-nihilations. We also add the resonant photoproduction off muons which dominates the rate up to m Z ′ ≲ 10 MeV, where neutrino annihilation starts to give the largest contribution to the rate, producing the characteristic strengthening of the bound with m Z ′ . At m Z ′ = 2 m µ a small feature signalizes the onset of resonant µ + µ − annihilation. Instead, the boundary of the constrained region in the deep trapping regime is dominated by inverse resonant annihilation into neutrinos. Note that the χ mass scales with the Z ′ mass, which has the effect of slightly suppressing the rate as compared to the masslessχ case, cf. Fig. 4. A rough estimate of the uncertainty of the excluded region is indicated with a dashed orange line, which shows the limits obtained from employing the hottest SN simulation SFHo-20.0, as described in Sec. II.

V. SUMMARY AND CONCLUSIONS
In this paper, we have studied the SN 1987A cooling constraints on dark-sector models induced by the emission of new light dark fermions χ coupled to leptons. To provide a concrete framework, we consider general vector portal interactions arising from the exchange of a massive Z ′ vector mediator. We focus primarily on the couplings to muons, which are predicted to have sizable number densities within the hot and dense environments of the proto-neutron star formed during core-collapse supernovae. However, we also extend our analysis to couplings to neutrinos and electrons.
We have considered various mechanisms for the production and absorption of the χ's and different regimes that depend on the ranges of the parameters of the model. Firstly, the constraints depend on the mass of χ, as their pair production becomes Boltzmann suppressed for m χ ≫ T . Secondly, different regions of m Z ′ can be identified based on whether the dark fermions are resonantly produced or generated from the tail of the Z ′ resonance. This distinction arises, for instance, when the Z ′ is heavy and cannot be produced on-shell by the thermal fluctuations in the medium. Finally, there exist two distinct regimes of coupling values, depending on whether the dark fermions free-stream or become trapped within the PNS. Consequently, by analyzing the χχ production within these two regimes, for given masses of the Z ′ and the χ, we can determine the range of couplings that is excluded by the observations from SN 1987A.
For Z ′ particles with masses m Z ′ ≲ T ∼ 10 MeV and massless χ, the observations from SN 1987A place constraints on the couplings between ∼ 10 −1 and ∼ 10 −9 , for equal couplings of the vector mediator to leptons and dark fermions, cf. Fig. 4. However, the range of these bounds strongly depends on the Z ′ mass and the specific lepton to which the Z ′ couples.
These calculations can be readily extended to explicit Z ′ -models, for example motivated by a gauged lepton flavor symmetry with a dark sector charged under it. We specifically investigate the case of L µ − L τ , which has been proposed in the literature as a combined solution to the (g − 2) µ anomaly, the Hubble tension and the dark matter puzzle. The SN 1987A limit covers a large region of parameter space that overlaps with the forecasts of future experiments and with part of the region that could address some of the tensions, see Fig. 5.
On the other hand, when the Z ′ mass is larger than the temperature and chemical potentials in the PNS, the interactions mediated by the Z ′ can be accurately described by effective operators. This allows us to generalize the analysis to completely generic heavy portal interactions between dark fermions and SM leptons, summarized in Fig. 3. We find that SN 1987A cooling can probe new-physics scales up to 4 − 7 TeV (cf. Table I), which surpasses current bounds from laboratory experiments by an order of magnitude (see Fig. 2).
We emphasize that our analysis is not complete when a light Z ′ is coupled to electrons. In this case, bremsstrahlung processes are expected to provide the dominant contributions to the emission of dark fermions, and plasma effects can have a significant impact on the analysis. Nevertheless, we consider our results as an important step into this direction, which significantly extends previous studies in the literature.
Finally, other aspects of SN physics could lead to constraints complementary to the ones obtained in this work (see e.g. Refs. [185,186]). In particular, there have been recent efforts towards a better understanding of the effect of neutrino self-interactions [159,161,162] in the trapping regime. Moreover, it would be interesting to check whether bounds obtained from energy transport by dark particles exceed the bounds derived from energy loss alone, following the analysis in Ref. [185].

CODE AVAILABILITY
We provide a minimal code example to test different parameter points of the L µ − L τ model at https: //github.com/spinjo/SNforMuTau.git.

ACKNOWLEDGMENTS
We thank Andrea Caputo, Pilar Coloma, Miguel Escudero, Patrick Foldenauer, Felix Kahlhöfer, Enrico Nardi, Uli Nierste, Filippo Sala and Stefan Vogl for useful discussions. We want to thank Robert Bollig and Hans Thomas Janka for sharing data of the simulations with us. This work is partially supported by project C3b of the DFG-funded Collaborative Research Center TRR257, "Particle Physics Phenomenology after the Higgs Discovery" and has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska -Curie grant agreement No 860881-HIDDeN. The work of C.A.M. is supported by the Office of High Energy Physics of the U.S. Department of Energy under contract DE-AC02-05CH11231. Work by JMC is supported by PGC2018-102016-A-I00, and the "Ramón y Cajal" program RYC-2016-20672.

Appendix A: Generalities of absorption rates
In the case of generic 2 → n processes, where a dark particle χ interacts with a particle b in the initial state, the resulting contributions Γ b χ to the absorption rate of χ can be approximated by where ⟨·⟩ b denotes the thermal average taken over the particle b, defined by Furthermore, σ denotes the cross section of the process, the Møller velocity and we have approximated the Pauli-blocking effects by introducing the degeneracy factors defined in Eq. (32).
By performing the d 3 p χ integral in Eq. (6), one arrives at where ⟨σv⟩ χ b denotes the thermal average over the complete initial state kinematics, i.e.
The inverse process defines an analogous collision operator for production of χ's, C b prod , which in conditions of thermal and chemical equilibrium reads This detailed balance relation can then be used to estimate the MFP in the PNS Appendix B: Cross sections for annihilation (2 → 2) In this appendix, we list cross sections for the 2 → 2 processes ℓℓ → χχ, χχ → ℓℓ (s-channel) and ℓχ → ℓχ  13 We disagree with Ref. [128] on the t-channel recovering their results only in the limit m ℓ , mχ → 0.
We can also generalize the cross sections for the Z ′ model in Eq. (1) by including generic vector (V i ) and axial (A i ) couplings to the leptons (such that V ℓ = 1, A ℓ = 0 gives back the model in Eq. (1)): The t-channel for small Z ′ masses is more involved because the propagator depends on the Mandelstam variable t, over which one integrates to obtain the total cross section. Neglecting the Z ′ decay width, the resulting ex-pression reads The cross sections for the inverse process χχ → ℓℓ can be obtained using In the limit m χ → 0, the final integral can be evaluated analytically, where g S = − 79ŝ 6 − 14ŝ 5 − 189ŝ 4 − 296ŝ 3 + 527ŝ 2 + 54ŝ − 3 + 12ŝ 2 2ŝ 4 + 10ŝ 3 + 9ŝ 2 + 44ŝ + 25 logŝ , g P = − 79ŝ 6 + 338ŝ 5 + 675ŝ 4 − 1160ŝ 3 + 175ŝ 2 + 54ŝ − 3 + 12ŝ 2 2ŝ 4 + 2ŝ 3 − 63ŝ 2 − 28ŝ + 17 logŝ . (C12) The matrix element for vector-like interactions reads After squaring and spin-averaging, one again obtains two traces for the SM and dark part of the amplitude, respectively. Due to the vector-like nature of the interaction, these traces are contracted with two Lorentz indices, one arising from M and one from M † |M| 2 = X µν (p 1 , p 2 )L µν (p a , p b , p 3 ) .