Supernova Muons: New Constraints on Z' Bosons, Axions, and ALPs

New light particles produced in supernovae can lead to additional energy loss and a consequent deficit in neutrino production in conflict with the neutrinos observed from Supernova 1987A (SN1987A). Contrary to the majority of previous SN1987A studies, we examine the impact of $Z'$ bosons, axions, and axion-like particles (ALPs) interacting with the muons produced in SN1987A. For the first time, we find constraints on generic $Z'$ bosons coupled to muons, and apply our results to particle models including gauged $L_\mu-L_\tau$ number, $U(1)_{L_\mu-L_\tau}$, and gauged $B-L$ number, $U(1)_{B-L}$. We constrain $Z'$ bosons with masses up to about 250-500 MeV, and down to about $10^{-9}$ in $Z'$-muon coupling. We also extend previous work on axion-muon couplings by examining the importance of loop-level interactions, as well as performing calculations over a wider range of axion masses. We constrain muon-coupled axions from arbitrarily low masses up to about 200-500 MeV, with bounds extending down to axion-muon couplings of approximately $10^{-8}$ GeV$^{-1}$. We conclude that supernovae broadly provide a sensitive probe of new lightly-coupled particles interacting with muons.


Introduction
Deep in the Large Magellanic Cloud, on the outskirts of the Tarantula Nebula, a blue supergiant named Sanduleak once shone with the brightness of over ten thousand Suns. That was, until one day, no more nuclear energy could be released by fusion, hydrostatic burning ended, and the star's iron core began to collapse under its own gravity. Eventually, the core became so massive that even electron degeneracy pressure could not stabilize it. As Sanduleak contracted, photons began to dissociate the iron atoms of the inner core, which decreased the energy of the star even more and caused it to further contract. Electrons were absorbed onto protons and converted into neutrons and electron neutrinos. The escape of the neutrinos further lowered the electron degeneracy pressure, until Sanduleak's core became unstable and collapsed.
When Sanduleak's core reached nuclear densities ρ = 3 × 10 14 g cm −3 , its equation of state stiffened [1] and the collapse was halted [2]. A shock wave formed between the outer and inner core and moved outward from the core through the star; the implosion of the inner core ignited an explosion and a core-collapse supernova was initiated 1 . A mere 0.3 seconds after the collapse, the shock wave ejected the entire contents of Sanduleak's outer layers, leaving only a compact remnant behind. The electron neutrinos from the burst were observed in 1987 by the Kamiokande II [6], IMB [7], and the Baksan [8] collaborations, in an event known as Supernova 1987A (SN1987A).
The observation of SN1987A has provided a unique opportunity to gain insight into the remarkable process that some massive stars undergo at the end of their life: the collapse of their core and the subsequent formation of a neutron star (NS). The neutrinos observed from SN1987A have yielded a wealth of knowledge in the fields of nuclear physics and particle astrophysics. Remarkably, the observation of neutrinos from SN1987A has also allowed the particle physics community to constrain the properties of new Beyond the Standard Model (BSM) particles. If new weakly interacting, long-lived particles were produced within the supernova, their escape would have carried energy away from the star. This new dissipation channel would have contributed to the cooling of the proto-Neutron Star (PNS) remnant of SN1987A [1]. As neutrino observations were consistent with simulations studies with neutrinos as the only cooling mechanism [2,9], a novel cooling channel does not appear to have dominated, implying constraints on the coupling between any new particle and the SM.
Most of the studies of SN1987A cooling have focused on new particles, such as dark photons and axions, that interact with nucleons, electrons, neutrinos, and photons (which are plentiful in supernovae) [1,[9][10][11][12][13][14][15][16][17][18][19][20][21]. In this study, we instead focus on the phenomenology of new particles interacting with muons. While PNS are born with a negligible number of muons, the dense (ρ ∼ 3 × 10 14 g / cm 3 ), hot (T ∼ 40 MeV) environment allows for the production of a sizeable population of muons through weak interaction processes with electrons, thermal photons, and neutrinos [22]. After the muon abundance equilibrates (which happens very shortly after the beginning of the neutrino cooling phase), any new particle sourced by muons can be produced in SN1987A and contribute to its cooling. We will focus on a generic Z -muon coupling, and a generic coupling between muons and axions or axion-like particles (ALPs).
As specific examples of the Z -muon coupling, we will apply our results to some U (1) models, including gauged L µ − L τ , U (1) Lµ−Lτ [23][24][25] (which is a particularly appealing leptophilic model since it is anomaly free), as well as gauged B − L number, U (1) B−L [26][27][28][29]. While limits on the muon-axion coupling were considered in a recent paper [30,31] (and estimated in an earlier paper [32]), we find stronger constraints. We also extend Ref. [30,31] by studying the impact of loop-induced processes, which generically can be expected to arise, and by performing calculations to higher axion masses.
This paper is organized as follows. In Section 2, we describe the processes which give rise to muons in supernovae. In Section 3, we discuss the impact of new particles on the supernova process, and outline a framework for our calculations. We then discuss the Z interactions and U (1) models that we consider in Section 4, detailing the relevant cross sections, processes, rates, complementary constraints, as well as our new results. We discuss the axion model and rates in Section 5, along with relevant constraints, and our new results. We conclude in Section 6. Appendix A contains details of additional particle rates.

Muon Production in Supernovae
The evolution of the core-collapse of SN1987A is well studied: simulations, informed by observations, yield a good understanding of the properties of the PNS (see for instance Ref. [3,4] and Chapter 11 of Ref. [1] for useful reviews). Of key interest to us, the PNS may reach core temperatures as high as 50 MeV, before it is cooled by neutrino diffusion in the first ten seconds after the bounce. The PNS is born with a significant electron to baryon number ratio (from the core of the progenitor star) but no initial muon or tau population. 2 When hydrostatic burning within a massive star ends and the star's iron core grows to a mass close to the Chandrasekhar mass, electron degeneracy pressure can no longer stabilize the core against gravitational collapse, and the core-collapse supernova is initiated. Electrons and protons in the progenitor star combine in a process known as neutronization. The ν e quickly diffuse out of the star, decreasing the net lepton number of the star. So begins the process of forming the neutron-rich core of the NS.
Due to weak magnetism corrections [35], the interaction cross section for matter and neutrinos is slightly larger than for anti neutrinos, i.e.
We make the common assumption that neutrino oscillations may be neglected [22] so that Le, Lµ, and Lτ are separately conserved. Incorporating oscillations is beyond the scope of this work, as it requires a simulation to be run on an order second timescale, whereas existing simulations have timescales of order nanoseconds [33,34].
where N = n, p. This implies that anti-muon-neutrinos diffuse out of the star faster than muon-neutrinos, building up muon number. A net excess of electrons over positrons occurs due to the high initial electron fraction of electrons compensating the positive charge of the protons. These electrons are highly degenerate and have a chemical potential µ e > m µ and can therefore be converted into muons via SM processes [22]. Furthermore, a large population of thermal photons and neutrinos, pair produced in the core, can also lead to a sizable production of muons. Through the combination of these effects, a large population of supernova muons is expected to be produced within the PNS (in contrast, τ leptons are too massive to be produced in substantial numbers at the temperatures reached in the core of the PNS). Importantly for the present study, the muons produced in this way will be significantly less degenerate than the electrons in the supernova [30], such that Pauli-blocking has a much smaller effect on muonic interactions.
Like the electrons, once a significant thermal population exists, the muons will undergo beta processes, where = e, µ. The produced neutrinos and anti-neutrinos diffuse out of the star, and from Eq. (2.1), anti-neutrinos will diffuse out faster than neutrinos. The result is a net muon number in the star-a process referred to as "muonization" -an excess of µ − over µ + . The effect of muons on the evolution of SN1987A was studied in [22]. The muon density depends sensitively on the core temperature and the electron chemical potential. Recent results from gravitational wave studies [36] and NICER [37,38] imply a relatively soft neutron star Equation of State (EoS), which in turn implies relatively high core temperatures [30]. The production of muons effectively softens the EoS additionally, as thermal energy is converted into rest energy of the muons [22]. We will show results for the SFHo EoS [39] applied to progenitors of 18.8 M and 20.0 M ; the temperature and muon number density profiles are determined from simulations in Ref. [30] 3 . In this way we may obtain a conservative and an optimistic constraint respectively.
The simulation snapshots in Ref. [30] consider the muon number density profile at a time stamp of t = 1 second after the bounce. Generically, simulations show that the temperature peaks at roughly 0.5−1 seconds post bounce, attaining a value of roughly 50 MeV at a distance of roughly 10 kilometers away from the center of the core [1,2]. After the first second, the peak temperature stays nearly constant, although the location of the temperature peak moves slowly inward as the PNS discharges neutrinos. Therefore, we expect that muon production does not vary significantly, and we consider n µ at t = 1 second to be representative of the entire 10 second cooling phase.

Impact of New Particles on SN1987A
If new particles are able to transport energy out from behind the neutrinosphere R ν , the radius outside of which most SM neutrinos freestream until arriving at Earth, less energy will be available to produce energetic neutrinos. If sufficiently large, this energy loss would reduce the duration of the neutrino burst below the expectation of 10 seconds [10], which in turn enables constraints on the new particle interaction strength. We now discuss in a general way how we compute constraints on the dark sector particles, following the prescription in Refs. [1,16,17].

Production
We define Γ prod as the rate of production of new long-lived particles (LLP) in the PNS. In the free-streaming region (i.e., when there is no reaborption of the LLP and they all escape the star), the luminosity of the LLP is given by integrating over the differential power

2)
V is the volume inside of the neutrinosphere, m LLP is the mass of the LLP, and ω is the energy of the LLP.

Trapping
Once the LLP interaction strength becomes larger, we move towards the trapping regime. This is when effects such as rescattering and reasborption become important and prevent the novel particle from streaming freely out of the star. To understand these effects, we first consider the absorptive optical depth τ abs for the LLP, which is given by where Γ abs is the absorptive width of the LLP, which is related to the production rate by Γ prod = e −ω/T Γ abs by the principle of detailed balance [40]. R far ∼ O(100) km is a radius outside of which neutrinos are no longer produced efficiently (we use R far = 100 km). The luminosity of the LLP is given by where R ν ∼ O(25km) is the neutrinosphere (at timestamp 1 second post bounce). Following Ref. [16], we define R ν to be the radius where the temperature of the star has fallen to 3 MeV, which is approximately consistent with the condition for neutrino free streaming.
We see that Eq. (3.4) reduces to Eq. (3.1) in the limit τ abs → 0, but we emphasize that Eq. (3.4) is applicable in generality, for all optical depths. In the large-τ abs limit, we expect that Eq. (3.4) should approximately asymptote to the blackbody luminosity for a radius at which the density of source particles (in our case, muons) becomes suppressed.

Cooling
The production of new weakly coupled particles contributes to the cooling of the PNS, which changes the neutrino emission from SN1987A. From observations, it is known that the cooling time of SN1987A is at least ten seconds when driven by neutrino cooling [2,9], but would be shorter if a new energy sink was present [10]. This translates into an upper bound on the luminosity of the new particles referred to as the "Raffelt criterion": the luminosity carried away by the new long lived particle from within the PNS to the outside of the neutrinosphere, L LLP , must not exceed the luminosity that would be carried away by neutrinos, L ν [1]. For SN1987A, this criterion is given by: This criterion will lead to limits on the size of the coupling between the dark particles and the SM, such that the new particle does not transfer energy more efficiently than the neutrinos. We emphasize that the neutrino luminosity on the right-hand side of Eq. (3.5) is a characteristic value at a time 1 second post-bounce. The actual neutrino luminosity is subject to change as a result of different simulation conditions; for example, the neutrino luminosity appears to be slightly enhanced in simulations with muons [22], but the criterion is not significantly affected [30]. Because the supernova muon number density is exponentially suppressed at a radius inside of the neutrinosphere, the large-τ abs limit of Eq. (3.4) can be a blackbody luminosity that is brighter than the neutrino luminosity. Thus, particles sourced by muons and having no other appreciable couplings to the SM may have no "ceiling" to the region in which they are excluded [30]. This approach to the blackbody limit is not exact, however, because of the mismatch between R ν and R far in our Eqs. (3.3) and (3.4) and because of the energy-dependence of τ abs . We explore this phenomenon in greater depth in the following sections.

Z Bosons
In this Section we will consider an LLP identified with a gauge boson, which we refer to in all cases simply as a Z . We discuss several models, calculate relevant cross sections, and present constraints on the strength of the interaction between the Z and the SM.
Throughout this section, we will remain agnostic as to the mass generation mechanism for the Z . Note that an additional Higgs boson, which may be required for gauge invariance for some masses and couplings (see e.g. Refs [41][42][43][44][45][46][47]), can modify the phenomenology; additional couplings to an additional Higgs could produce additional important rates, which may dominate for particular combinations of parameters. In the present paper we only consider vector couplings of the Z , and therefore an additional Higgs mechanism is not strictly required, but could be an interesting extension for future work.

Muon-Coupled Z Models
Here we present the details of the Z models probed by this work. The models we consider have the following the tree-level interaction terms in common, which arise from gauging muon number (see e.g. [48][49][50][51]), where g µ is the Z -muon/muon-neutrino coupling, = (µ L , ν µL ) is the electroweak doublet for the left-handed leptons, and the singlet for the right-handed muon is µ R . In what follows, we will summarize in more detail some examples of models in this form.
We consider extending the SM by the gauge group U (1) Lµ−Lτ . This model is particularly appealing as it is anomaly free. The Lagrangian is given by where g Z is the U (1) Lµ−Lτ gauge coupling, and the electroweak doublets for the left-handed leptons are 1 = (µ L , ν µL ) and 2 = (τ L , ν τ L ), and the singlets for the right-handed leptons are µ R , τ R . At tree-level, the Z in this model only couples to muons, taus, and their neutrinos. However, an effective coupling for the Z to all electromagnetically charged fermions can be induced due to muon and tau loops. The kinetic mixing term above, − 2 Z µν Z µν , arises due to integrating out the muons and taus in the low-energy limit. This leads to irreducible contributions to from the muon and tau loops. We can estimate the amount of mixing likely to arise by [52] = This can be considered a "natural" amount of kinetic mixing. However, the interplay of processes induced by kinetic mixing with those arising from the gauge coupling g Z will depend on the full UV model. To demonstrate variance of results, we will consider three different amounts of kinetic mixing, which are phenomenologically distinct: no mixing ( = 0), one where the couplings are comparable ( = −g Z ), and the natural estimate ( = −g Z /70).

Gauged B − L Model
The last Z model we consider is B − L, which arises from gauging the combination of baryon number B minus lepton number L, as the gauge group U (1) B−L . This model is anomaly free if right-handed neutrinos are included. The popularity of this model has been largely due to providing a simple framework for the implementation of the seesaw mechanism, and its consequent ability to explain the smallness of neutrino masses. We define the Lagrangian for the B − L model as where ν R are additional right handed neutrinos of each flavor (added for anomaly cancellation, although not strictly required), are the left-handed lepton doublets, and q are the left-handed quark doublets.

Cross Sections and Rates
We now consider the rates relevant for Z production using supernova muons and neutrinos. Figure 1 shows the processes in which the Z can be produced from interactions with the muons or neutrinos produced in the supernova. The two main rates of interest that dominate the phenomenology arising from muon interactions are neutrino-pair coalescence (νν → Z ), and the semi-Compton scattering (γµ → Z µ). For estimates of other rates, and why they are subdominant, see Appendix A.
To determine the energy-dependent absorptive width for neutrino-pair coalescence, we have where Z → νν is the decay rate of the Z to neutrinos. By the principle of detailed balance [40], we can substitute the relation f ν (p 1 )f ν (p 2 ) = f Z (p Z ), and use the decay rate Γ Z →νν = αm Z /3, to find the rate of inverse Z decay -which is the neutrino-pair coalescence rate Here α Z = g 2 Z /4π is the new fine structure constant for the gauge group under consideration, with gauge coupling g Z . We assume that m Z 2m ν , such that decay and coalescence are on-shell with unsuppressed phase space. This is violated for masses m Z 1 eV, but, as discussed below, (inverse) decay is negligible in this part of parameter space.
To obtain the rate for semi-Compton scattering, we use the well-known SM Compton cross section of σ (µ) 4 We multiply this by the Z velocity, by the muon number density, and by an additional degeneracy factor F deg to take into account Pauli blocking of the muons. This gives a semi-Compton Z production rate of We use the values of F deg for muons as simulated in Ref. [30] as a function of radius within the supernova. Finally, for m Z ≥ 2m µ , the decay Z → µ + µ − can be on-shell, and pair-coalescence can lead to Z production. This results in a production rate The Bose-enhancement of the production compared to the rate calculated in [52] arises because the muons are in thermal equilibrium in the PNS.

Existing Z Constraints
In the parameter space of interest, we also consider constraints relevant for the Z -muonnumber coupling which arise from other experiments or observations. Some are applicable to the tree-level Z -muon coupling, while others are only applicable when there is some kinetic mixing present. Constraints which are relevant in some portion of parameter space include: • Stellar cooling: The production of new light particles in stellar cores affects energy transport within the star. Therefore, the non-observation of such anomalous energy transport can be used to constrain light degrees of freedom. For kinetically mixed bosons, particular care needs to be paid to plasma effects. The constraints shown below are rescaled from Refs. [53,54].
• Solar Neutrino Scattering: A gauge boson Z , kinetically mixed with the SM photon, may alter the solar neutrino scattering rate, due to its interaction with the muon and taus (as well as muon and tau neutrinos). The Borexino experiment is sensitive to this scenario [55][56][57][58]. Here, we adapt the limits recently calculated in Ref. [58], where we show the conservative case for the low metallicity Sun.
• Muon Magnetic Moment: While a new gauge boson can explain the anomalous contribution to (g − 2) µ [59], a gauge coupling that is too large is excluded by this measurement. The resulting constraints were calculated in Refs. [52,60].
• Effective number of neutrino species: Additional light degrees of freedom may appreciably change the number of relativistic degrees of freedom N eff at recombination. The contribution to N eff from Z gauge bosons depends on the mass and coupling of the Z . For sufficiently large coupling, g Z 4 × 10 −9 , the gauge boson would have been thermalized with the SM at early times [52]. Its decay, in particular if it occurred after neutrino-photon decoupling, implies an effective number of neutrino species potentially in tension with existing constraints. Additionally, in the presence of kinetic mixing, additional interactions may also delay neutrino-photon decoupling itself. Experimentally, N eff is constrained by several effects. Firstly, additional relativistic species modify the expansion rate until just before recombination. This leads to the constraint, with error bars spanning the 95% confidence interval, N eff = 2.97 +0. 58 −0.54 , from CMB polarization and baryon acoustic oscillations (TT, TE, EE + lowE + lensing + BAO) as measured by the Planck satellite [61]. However, this constraint on N eff assumes the H 0 value as measured by Planck, which is in tension with low-redshift measurements. As an additional number of relativistic species modifies the expansion history, the resulting ∆N eff may relieve some of that tension [62,63]. Taking this into account, the central value shifts N eff = 3.27 ± 0.30 [61]. A complementary constraint comes from light element abundances produced during Big Bang Nucleosynthesis (BBN) [64], which would be affected by the annihilation of new light degrees of freedom heating up the photon bath. The 95% confidence interval from light element abundances is given by, N eff = 2.95 +0. 56 −0.52 [61]. In light of these constraints, we show a dashed line for ∆N eff = 0.5 in our plots.
• Neutrino Tridents: A muon neutrino scattering off of the Coulomb field of a nucleus may lead to the production of a µ + µ − -pair; this rare process is referred to as a neutrino trident. New muonphillic gauge bosons alter neutrino tridents through the exchange of a photon and a Z boson. The resulting constraints for L µ − L τ gauge bosons were calculated in Refs. [65,66]. This is probed by the Columbia-Chicago-Fermilab-Rochester (CCFR) Experiment. We label these constraints "CCFR" in the plots.
• Beam Dump Experiments: For the B − L model, there are a number of relevant beam dump constraints in our region of interest, extending down to Z masses of about an MeV. These include electron beam dumps SLAC E137 and SLAC E141 [67][68][69][70], Fermilab E774 [71], and Orsay [72]. Relevant proton beam dumps include LSND [73] and U70/NuCal [74,75] (see also e.g. Refs. [76][77][78] for relevant summaries of these searches for our models).  Figure 2 shows the luminosity L Z as a function of the coupling g Z , using the rates described in the previous subsection and using Eq. (3.4) for all values 5 of g Z . We show the luminosity for the Z model without kinetic mixing (the result is qualitatively similar for the other model variations). On the left-hand side of the plot, it is seen that for m Z 10 MeV, where the production in the core is dominated by semi-Compton scattering (Eq. (4.7)), L Z is m Zindependent. For this mass range, L Z exceeds L ν (defined in Eq. (3.5)) when g Z ∼ 2 × 10 −9 . For m Z 10 MeV, where neutrino-pair coalescence Eq. (4.6) dominates the production, the value of g Z for which L Z = L ν is reached for smaller values of the coupling.

Z Boson Results and Discussion
Many interesting features are evident in Fig. 2: we see the effect of Boltzmann suppression in the reduction of the height of the peak of L Z as m Z MeV; we see the transition between the mass-dependent and mass-independent regimes of L Z ; and we see an interplay between muon-rescattering (the "slanted plateau" shape for low m Z ) and decay contributions to optical depth (the exponential decay at sufficiently large m Z and g Z ). We also demonstrate the effect of reducing the value of R far in Eq. This plot is also applicable to a muon-only coupled Z , though in that case the bounds are weaker by a factor of √ 2 above masses of about 10 MeV (there are also different complementary constraints for that case, see text). New SN1987A limit is shown as magenta shaded region. Inner region is conservative calculation, outer is less conservative. more energy appears to leak out from behind the neutrinosphere. However, this energy will be deposited behind the neutrino gain radius, leading to uncertain impacts on the neutrino luminosity, shockwave formation and revival, and PNS evolution. To be conservative, we do not consider such energy transport to contribute to anomalous cooling, and we always use R far = 100 km in our results. Figure 3 presents our calculated constraint on L µ − L τ Z bosons (when no kinetic mixing is present) arising due to Z production in SN1987A through interactions with muons. The darker/lighter magenta shaded regions are calculated using the SFHo-18.8 and SFHo-20.0 simulations [30], respectively, corresponding to more and less conservative bounds. The range between the two regions can be considered a systematic uncertainty on the muon number and temperature profiles. We include additional relevant constraints in this region of parameter space: "CCFR" which arises from neutrino tridents (dark green), as well as the best-fit region to explain the (g − 2) µ anomaly (bright green). Note that regions above the (g − 2) µ best-fit are excluded for the muon models, as Z bosons with such large couplings would be in conflict with the (g − 2) µ measurement. In this scenario, with no kinetic mixing, our new supernova constraint is the strongest existing bound. The main other competing constraint, which is shown as the dashed line in the figure, is N eff , the effective number of light degrees of freedom in the early Universe [52]. In the standard scenario, we could expect a constraint on the coupling of about 10 −8.5 for masses below about 10 MeV, as such masses and couplings could lead to a change in N eff of about 0.5. Note that additional bounds from Planck may also be relevant for masses up to about 5 × 10 −4 MeV [79]. Previously, the SN1987A bound from earlier work, which had not incorporated charged muons [52], was strongest for Z masses greater than about 10 MeV. Other than these bounds, any remaining competing constraints are substantially weaker.
The shape of our Z supernova bounds in Fig. 3 can be understood due to the interplay of the semi-Compton and pair coalescence processes depicted in Fig. 1 (also see the luminosities in Fig. 2). For Z masses larger than about 10 MeV, the pair coalescence rate dominates, and corresponds to the drop in the limit curve at these masses. This can be compared to the previously estimated SN1987A bound on L µ − L τ gauge bosons [52], which is shown as the dotted line and labelled "Previous SN1987A" in Fig. 3. This previous estimate used only the pair coalescence rate, which is relevant only for the neutrinos in the supernova, and not the muons. As such, the limit found in Ref. [52] roughly traces the limit we would obtain if also only considering muon neutrinos, with the difference arising from different temperature profiles. As we have also included the effects from muons in the supernova, the additional semi-Compton process is relevant, and is what leads to the flat lower end of the bound for Z masses below about 10 MeV. We have shown this constraint to extend down to 10 −8 MeV in the plot; we expect that Z production will begin to experience a suppression due to plasma screening in the core at roughly this mass, since ω µ−τ p 4πα Z n ν {µ,τ } /T ν g Z T ν . The values of this limit would then continue to lower mass as m Z g Z ∝ constant. The upper ends of the bounds are controlled by the rate of trapping of the Z within the supernova. Once considering masses below about 10 −4 MeV, the upper end of the bounds also become flat, reaching a consistent upper limit of about 10 −3 . At higher couplings, we expect there to be anomalous energy transport in the supernova and PNS, but we do not expect excessive energy loss, as shown in Fig. 2. Even with this caveat, our bounds exclude the explanation of the (g − 2) µ anomaly for this model for Z masses below about 10 −5 MeV, in the mass range shown.
As the τ leptons are too heavy to be thermally produced in SN1987A, it may seem that the lepton doublet effectively decouples. However, there will be ν τ neutrinos present; while these do not contribute to the re-scattering of the Z (due to G F suppression), they can coalesce to produce the Z or be produced in its decay. This effect simply contributes an additional factor of 2 in the optical depth τ abs and the differential power Q with respect to the decoupled scenario. Hence, for a model coupled only to muon number, the L µ − L τ results in the coupling-mass plane would be rescaled by √ 2, in the right hand dip. However, we note that any new gauge symmetry involving only muon number needs to involve another fermion doublet with electroweak as well as U (1) charge for anomaly cancellation. Such scenarios are typically very constrained: if the fermions are heavy, they may induce processes like flavor changing meson decay [80,81], if they are light, they are constrained by fourth-generation searches at the LHC [82]. Figure 4 shows our calculated bound for the U (1) Lµ−Lτ model using the muons in the supernova, when kinetic mixing is included. Additional complementary constraints are relevant in the presence of kinetic mixing. We show the bounds from stellar cooling, which are obtained from rescaling the results in Ref. [53]. We also show the limits from Borexino, which are also rescaled simply with the coupling relation. When the amount of kinetic mixing is varied relative to the U (1) Lµ−Lτ gauge couping, both the stellar cooling results and Borexino limits also vary accordingly. The results shown feature two different non-zero two choices of epsilon, = −g Z in the top panel, and = −g Z /70 in the bottom panel. These choices are made to demonstrate how the phenomenology may change in the presence of different amounts of kinetic mixing, as the exact amount expected is highly UV-model dependent. In the case where = −g Z , the main differences are that the complementary bounds become stronger, and most importantly, the supernova bound features an additional peak around a Z mass of about 20 MeV. This is because an additional process becomes relevant for significant kinetic mixing: the bremsstrahlung process from protons, np → npZ , (a process which does not require supernova muons). The peak shown comes from a resonance of the mixed Z , and the whole rate can be compared by simply rescaling the dark photon bremsstrahlung process in Ref. [16]. If the kinetic mixing parameter were larger, the peak at 20 MeV would move down proportionally. On the other hand, if the kinetic mixing parameter is sufficiently small, like in the bottom panel where = −g Z /70, the supernova bounds for the U (1) Lµ−Lτ would not change, and only the complementary constraints would change.
For other leptophilic models, such as U (1) L or a U (1) Le−Lµ model, the bottom of the supernova bounds shown in Fig. 3 and Fig. 4 would be equivalent. This is a consequence of the fact that tau leptons are too massive to be produced in the supernova, so their coupling only affects the constraints by enhancing the neutrino-pair coalescence rate by a factor of two. Moreover, while the electron number density is higher than the muon number density, the electrons are much more degenerate in the core. Therefore the semi-Compton rate involving electrons is highly suppressed by Pauli blocking, which implies that the production rate (which affects the floor of the bounds) would not be affected by the presence of electrons. On the other hand, while electrons are degenerate in the core, further out in radii (i.e. by ∼ 20km), we do not expect this suppression to be large. This means that rescattering on electrons (e.g. inverse semi-Compton absorption) could be important at the ceiling of the SN1987A constraints. As the impact of electron density profiles is outside the focus of this work, we do not plot this scenario. Figure 5 presents the new supernova limits for U (1) B−L gauge bosons. Importantly, compared to our other results, our new B − L bound does not largely stem from the presence of the charged muons, but rather from the neutrinos. The flat region in this case comes from rescaling the U (1) B model limits in Ref. [15], and the peak at around 20 MeV is again from a mediator resonance with the bremsstrahlung process from protons, np → npZ , which was rescaled from Ref. [16]. On the higher mass end (m Z 20 MeV), the dominant process is neutrino-pair coalescence (as before, we show a conservative and less conservative constraint, , which is considered a "natural" amount of mixing for the L µ − L τ scenario. New SN1987A limit is shown as magenta shaded region. Darker region is the conservative calculation, lighter outer region is less conservative (see text).  New SN1987A limit is shown as magenta shaded region. Inner region is conservative calculation, outer is less conservative (see text). Dashed line is previous SN1987A estimate in this region [83].
depending on the profile used). The inclusion of neutrino-pair coalescence explains the main difference with the previous version of this bound, reported based on the rescaled dark photon bounds from Ref. [16], from Ref. [83], which we also show. As the dark photon does not couple to neutrinos (whereas the B − L model does), the high mass region of the B − L constraint was not captured in Ref. [83]. Note that for masses m Z 20 MeV the hadronic interactions dominate over the conservative muon profiles. The darker shaded region for m Z 40 MeV is arising from a proton-related rate. A second difference with Ref. [83] is given by the extent of the resonance at m Z ∼ 20 MeV; we obtain this peak from the same reference (Ref. [16]) as Ref. [83]. This small variation is explained by our use of the conservative bound from Ref. [16], while Ref. [83] reports the fiducial bound. At lower masses, our results overlap with the constraints presented by Ref. [83]. We show complementary constraints from stellar cooling and BBN [83]. The beam dump experiments are shown as one uniform exclusion, though the region plotted contains all the experiments described in the previous section (and is the sum of the regions shown in Ref. [77]).

Axions and Axion-Like Particles
We now consider an alternative class of LLPs: CP-odd scalars known as axions or axion-like particles (henceforth we refer to both axions and ALPs as "axions"). Such particles are wellmotivated and appear in many BSM theories. However, constraints arising from the muon coupling of the axion are not well-studied. A recent paper, Ref. [30], considered the SN1987A limits arising due to the axion-muon coupling in the limit that m a MeV, focusing on tree-level axion-muon couplings. We will extend these results, by investigating the impact of loop-level axion-muon interactions on the parameter space shown in Ref. [30], which extends up to 1 MeV. We will then extend this calculation to higher masses, again including muon-loop interactions.

Muon-Coupled Axion Model
We consider the scenario where axions have muon-dominated interactions, coupled via L ⊃ g aµ (∂ σ a)ψ µ γ σ γ 5 ψ µ = −ig aµ (2m µ )ψ µ γ 5 ψ µ a, (5.1) where m µ is the muon mass, a is the axion field, ψ µ is the muon spinor, and g aµ is the (dimensionful) axion-muon coupling. This is the same model as considered in Ref. [30]. Figure 6 shows the dominant muon processes for the axion in the supernova, semi-Compton scattering (γ +µ → a+µ), and muon-proton bremsstrahlung (µ+p → µ+p+a). As discussed above, we will also include the rate for the axions to be produced/reabsorbed via a muon loop. The semi-Compton production rate is given by [84]

Cross Sections and Rates
where ω is the energy of the emitted axion, and F deg is a degeneracy factor arising from the pauli-blocking of muons, which was discussed above. F deg was simulated in Ref. [30] as a function of radius within the supernova. For bremsstrahlung of an axion in muon-proton collisions, the production rate is [84] Γ brem = α 2 EM (2g aµ m µ ) 2 8π 3 √ 2π Here w ≡ ω/T and y ≡ k S / 2m µ T where k S is the Debye screening scale, which is applicable for our non-degenerate muon scenario. At sufficiently large m a , the decay to two muons a → µ + µ − (or the reverse process: coalescence of two muons into a µ + µ − → a) can also be on-shell. The production rate for this process is The Bose-enhancement of the production compared to the rate calculated in [85] arises because the muons are in thermal equilibrium in the PNS.
As discussed above, we also consider an additional loop process, a → γγ as shown in Fig. 6, which was not considered in Ref. [30]. This process will contribute to both the production and decay rates. We calculate that the (boosted) production rate is where c γγ is a possible UV contribution to the aγγ vertex; we will take this to be zero. The function B 1 (τ ) is which determines how much the heavy lepton contributes to the loop. In the limit m µ m a , we find that B 1 → −m 2 a /12m 2 µ as stated in Ref. [85]. Note that the aγγ vertex vanishes as m a → 0, which differs from the non-decoupling case of CP-odd Higgs, e.g. Fig. 2.12 of Ref. [86], as discussed in Ref. [85].
We may check the importance of the a → γγ loop processes by first considering the m a m µ limit and letting c γγ → 0. In this limit, Eq. where in this mass range, |B 1 | ∼ O(0.1 − 1). Thus, the decay a → γγ can be important in the massive-axion limit.

Existing Axion Constraints
Unlike the Z model, the axion-muon coupling is not currently well probed. Besides the SN1987A constraints, the main other relevant constraints come from the anomalous magnetic moment of the muon (g −2) and from the additional relativistic species N eff , both described in more detail in Section 4.3. Since pseudo-Goldstone bosons contribute negatively to the muon magnetic moment [87], we show constraints for the 5σ deviation from the experimentally measured result [88]. In this case, the N eff constraint depends also on the coupling between axions and SM gluons and photons [89]. As this coupling depends on a UV-scale, and a single thermalized degree of freedom would at most be in very modest tension with the constraints in Section 4.3, we choose not to show it here. In the future, this coupling may be further probed by time-resolved reanalysis of ongoing and upcoming precession experiments [90]. Figure 7 shows the luminosity as a function of the dimensionful axion-muon coupling g aµ , using the rates in the previous subsection and using the SFHo EoS for a 20M progenitor [30]. At small m a m µ , the semi-Compton and bremsstrahlung processes imply a luminosity independent of the mass, exceeding the Raffelt criterion for g aµ 10 −11 MeV −1 . At m a m µ , the diphoton decay dominates the optical depth, and begins to reduce the parameter space in which luminosities in excess of Eq. (3.5) are attained. At m a ≥ 2m µ , the decay a → µ + µ − becomes dominant, leading to a dramatic reduction in the luminosity for a given coupling. The sharpness of this reduction comes from the fact that Γ µµ is tree-level whereas Γ γγ only occurs via a muon loop. Thus, as soon as a → µ + µ − goes on shell, the absorptive width increases very sharply; this is offset partially by µ + µ − → a. The figure also demonstrates the fall-off in luminosity due to Boltzmann suppression at large axion masses, and the asymptote to black body radiation for small axion masses. As in Fig. 2, the dashed line shows the effect of taking R far → R ν for a low-mass axion. The difference between the dashed and solid blue lines represents energy transport from behind the neutrinosphere into the region of shocked matter, and will lead to uncertain impacts on the supernova and PNS evolution. We use R far = 100 km, represented as the solid lines, in all constraints. Figure 8 shows our calculated bound for axions coupling to muons. For axion masses less than about 5 MeV, we find constraints on the axion-muon coupling of ∼ 10 −1.5 − 10 −9 GeV −1 . In this regime, the semi-Compton rate and muon-bremsstrahlung rates dominate the limit. The bound here is flat as these rates are independent of the axion mass. We have confirmed that the upper limit, where trapping occurs, is not impacted by the diphoton decay for small m a . While this loop-induced decay rate is non-zero for even small axion masses 6 , the optical depth over the scale of the proto-neutron star is negligible, as demonstrated in Eq. (5.9). At higher masses, above about 5 MeV, the upper coupling constraints begin to drop, and the lower coupling limits increase, as the loop-induced decay becomes important in this regime. Above m Z = 2m µ , the muon decay comes on-shell, and, except for a small range of couplings in the case of the 18.8M progenitor and a slightly wider in range for the 20M progenitor, this prevents the Z from traveling far enough to drain any energy from the PNS. Compared to the constraint originally obtained in Ref. [30], we find an axion-muon coupling bound that is stronger by about two orders of magnitude in the comparable mass regime. This can be traced to an incorrect factor of 2ω in Eq. 6 of Ref. [30]; we find that by incorporating this factor, we reproduce their results for both their luminosity as a function of coupling and their final results in the low-τ abs limit 7 . Separately, we have also considered the impact of axion-muon loop interactions which were not considered in Ref. [30], and we have shown that neglecting the loop-induced diphoton decay is valid in the m a MeV limit considered in that paper. We extend the constraint by computing the supernova results for the axion limits up to masses of 2m µ where loop-processes are indeed important.

Summary and Conclusions
The transformation of the blue supergiant Sanduleak into Supernova 1987A has gifted us with extensive insight into weakly-coupled particles. Such particles would have streamed out of the star, carrying away energy and cooling the star, leading to a potential tension with the observed neutrino fluxes from 1987A. Previously, limits on particles produced in this way mostly focused on interactions with particles that had a large pre-collapse abundance, such as the nucleons and electrons. We have instead focused on the impact of muonic production of new particles, exploiting recent simulations of the muon density and temperature profiles within the supernova.
For the first time, we considered supernova muons as a probe of muon-coupled Z bosons, setting strong limits on Z models including generic Z -muon couplings, U (1) Lµ−Lτ , and U (1) B−L . We used neutrino-pair coalescence and semi-Compton processes to produce the Z bosons, which can lead to energy loss in the supernova, and conflict with the number of neutrinos observed from SN1987A. We exclude the existence of these gauge bosons up to masses of about 250 − 500 MeV, and couplings down to g Z ∼ 10 −9 − 10 −10 for some masses. Unlike particles produced via kinetic mixing, the muonic production processes imply a flat constraint even at low masses, as they are not sensitive to plasma screening (until reaching below Z masses of about 10 −8 MeV). In the case that there is no kinetic mixing for the Z , these constraints are only rivaled by constraints on N eff in most of the parameter space, which rely on complementary physics. In the case that kinetic mixing is present, we examined the interplay of both proton/neutron and muon processes to produce the Z , finding that the muon-related bounds still dominate for most of the parameter space. Our new bounds exclude the L µ − L τ model as an explanation of (g − 2) µ for Z masses below about 10 −5 MeV, for our mass range shown. Our new constraints on U (1) B−L gauge bosons differ from previous estimates, with the difference largely arising from incorporating the neutrino-pair coalescence process at higher Z masses.
We also substantially extend the previous limits on axion-muon couplings. In particular, for the first time, we performed the calculation outside the m a MeV axion-mass limit, and also examined the impact of including axion-muon loop processes in the supernova calculation across all masses. We find that the loop processes determine the small coupling constraint at larger axion masses by an explicit calculation of the resulting optical depth. When examining the boosted decay rate of the diphoton process over the scale of the proto-neutron star, we find that the diphoton process is negligible for axion masses less than or equal to about 5 MeV. Above 5 MeV, we find that this diphoton process is important, and consequently used it when extending the axion bounds up to m a 2m µ . We find that we constrain the existence of a muon-coupled axion with masses ranging from arbitrarily low masses up to about 5 MeV, down to couplings of about 10 −8 GeV −1 . For axions with masses from 5 MeV up to m a 2m µ , limits gradually become weaker, with the axion becoming trapped and the bounds not extending up to larger couplings.
Across this work, we have pointed out and explicitly demonstrated the broad applicability of supernova muons to provide a sensitive probe of models of new physics. the parameter space, with the exception of regions in which 2 is large compared to α µ and m Z is not small compared to ω p .