Muons in supernovae: implications for the axion-muon coupling

The high temperature and electron degeneracy attained during a supernova allow for the formation of a large muon abundance within the core of the resulting proto-neutron star. If new pseudoscalar degrees of freedom have large couplings to the muon, they can be produced by this muon abundance and contribute to the cooling of the star. By generating the largest collection of supernova simulations with muons to date, we show that observations of the cooling rate of SN 1987A place strong constraints on the coupling of axion-like particles to muons, limiting the coupling to $g_{a\mu}<10^{-7.5}~\text{GeV}^{-1}$.


I. INTRODUCTION
Axions are hypothetical light pseudoscalar degrees of freedom.Initially proposed to solve the strong CP problem [1][2][3], they are a generic prediction of string theory [4] and a good candidate for dark matter [5][6][7][8][9].As a result, a large research effort is dedicated to searching for these "axion-like particles" (henceforth axions), mainly through their coupling to photons [7,[10][11][12][13], electrons [14], or nucleons [15] (see [16] for a review).The coupling to muons, however, has not been well studied, partially by virtue of the fact that the high muon rest mass and short lifetime make performing precision muon experiments over long timescales difficult.Luckily, astrophysics provides us with an alternate means to probe this coupling: the cooling of supernovae. 1s noted recently [17][18][19], though the temperature of a supernova hardly rises above ∼ 50 MeV, there is still a large proportion of thermal photons and neutrinos with energies well above the muon rest mass.Additionally, the electrons acquire a large chemical potential allowing their conversion to muons, with the production of muon antineutrinos compensating to maintain a net zero muon number.Due to the small difference in the neutral current cross-section between ν and ν, muon antineutrinos diffuse out of the PNS at a higher rate than muon neutrinos, leading to the accumulation of a net muon number within the PNS in a process known as "muonization."Recent work has shown this effect may play a nonnegligible role in facilitating supernova explosions [19].
The resulting muon abundance could produce a large axion flux.Here, we exclusively consider axions with mass m a ≪ MeV and muon-dominated interactions, coupled via2 L ⊂ g aµ (∂ σ a) ψµ γ σ γ 5 ψ µ ≡ −ig aµ (2m µ ) ψµ γ 5 ψ µ a , (1) where m µ is the muon mass, a the axion field, ψ µ the muon spinor, and g aµ the axion-muon coupling.Bounds on g aµ arise from axion contributions to the muon g − 2, at the level of g aµ ≲ 10 −2. 4 GeV −1 [20], as well as from cosmology [18,21,22] (see Supplemental Material for a detailed discussion of these limits).The bounds we calculate here extend the experimentally-excluded region by over five orders of magnitude.
Observations of SN 1987A neutrinos indicate that the resulting proto-neutron star (PNS) cooled in roughly ten seconds.A new particle transferring energy more efficiently than the neutrinos would shorten this timescale, leading to the oft-cited "cooling bound" [23][24][25][26][27][28][29] on the axion luminosity of L a ≲ 3 × 10 52 erg/s. 3ote that a rough estimate of this bound was made in Ref. [18], but due to its highly approximative nature, the result was subject to a large degree of uncertainty [22].This is not unexpected, as the high muon rest mass means both the muon density and reaction rates depend sensitively on the core temperature (and muon and electron chemical potentials) of SN 1987A.In this paper, we significantly improve on the previous estimate by running dedicated simulations that make use of recent results [35][36][37][38][39][40][41][42][43] constraining the PNS equation of state and mass, allowing us to reduce the uncertainty and set a robust bound on the coupling over five orders of magnitude.Interestingly, these new constraints on the EoS and mass lead to a generically higher core temperature, which suggests that existing bounds on other axion couplings may be strengthened by including this microphysics as well.

II. PROFILES FROM SIMULATIONS
As with any supernova bound, it is critical to assess the robustness of our results to variations in the choice of model.To do this, we ran four simulations that span a range of allowable final neutron star (NS) masses for SN1987A (see Table I).The simulations were performed in spherical symmetry (1D) with the Prometheus-Vertex code with general-relativistic corrections and six-species neutrino transport, solving iteratively the two-moment equations for neutrino energy and momentum with a Boltzmann closure [44] and using the full set of neutrino processes listed in [45] and [19].PNS convection was taken into account by a mixing-length treatment and explosions were artificially triggered a few 100 ms after bounce at the progenitor's Fe/Si or Si/O composition interface as described in [46].
The main astrophysical uncertainties are connected to the neutron star mass in SN 1987A and the equation of state at supernuclear densities.Model SFHo-18.6 has a canonical neutron star mass well within the range expected for the compact remnant in SN 1987A, while in models SFHo-20.0 and LS220-20.0 the neutron star masses are near the upper edge of the expected range, and in model SFHo-18.8 the neutron star mass is at the lower edge of the allowed range [34].The SFHo equation of state is fully compatible with all current constraints from nuclear theory and experiment [35][36][37] and astrophysics, including pulsar mass measurements [38][39][40] and the radius constraints deduced from gravitational-wave and Neutron Star Interior Composition Explorer measurements [41][42][43].Results for the long-used LS220 equation of state are shown for reference and comparison.Note that these equations of state are considerably softer than those of previous works on axion emission from supernovae [52][53][54], which employed stiff equations of state that are increasingly disfavored by the constraints above.The adoption of softer equations of state generically results in a smaller radius and higher temperature than in PNS models with stiff equations of state.It should be noted that though this change is most relevant for the axion-muon coupling, the generically higher temper- atures will influence many existing supernova limits on new particles, and these results should be revisited with the new EoS constraints in mind.
We have plotted the temperature, density, and muon number density for these simulations at 1 second postbounce in Figures 1, 2, and 3 respectively. 4Note that though the temperature varies by 30%, the ultimate muon density, the quantity to which our bound is most sensitive, does not change by more than an O(1) factor in regions of interest, demonstrating a considerable robustness to large changes in initial parameters.In order to place conservative limits, we ultimately adopt the SFHo-18.8 result as our fiducial profile, which is the coolest and results in the weakest constraint.Additionally, though we place our bound at t = 1 sec, we have also computed the bound at t = 0.5 sec and t = 3 sec, which span the time interval of highest temperature and neutrino luminosity, and find that the ultimate limit on g aµ does not change by more than a factor of two.

III. AXION PRODUCTION BY MUONS
There are two dominant contributions to the axion emissivity due to muons in the supernova: Compton scattering (γ + µ → a + µ) and muon-proton bremsstrahlung (µ + p → µ + p + a).Contributions from muon-muon bremsstrahlung are subdominant as the muon number density is over an order of magnitude below the proton density, muon-electron bremsstrahlung is suppressed by the muon-electron mass ratio, and other channels such as Primakoff or nuclear bremsstrahlung require additional couplings of the axion.
While electrons in the core of the PNS are highly degenerate, suppressing the Compton process [55] and resulting in bremsstrahlung [23] as the dominant axionproduction channel, muons are only mildly degenerate.This reduces the suppression considerably, allowing the Compton process to become the dominant contribution to axion production.This mild degeneracy is displayed in Figure 4, where we plot the ratio µ µ /T as a function of radius for the various profiles considered.Axion production is only effective in regions of high muon density (5 km ≲ r ≲ 15 km), where we observe that the degeneracy parameter remains a small O(1) value.
In order to quantify the effects of the degeneracy on the Compton rate, we compute the approximate suppression factor F deg by averaging over the Pauli blocking factor (see Sec. 3.2.5 of Ref. [23]): where n th µ is the thermal number density and f µ (E) = (e (E−µµ)/T + 1) −1 the energy distribution of muons at the given degeneracy and temperature.Due to the mild degeneracy, this suppression is never in excess of ∼ 15% (see Supplemental Material for a plot of F deg as a function of radius).
We also compute the degeneracy parameter for the protons, as previous work has suggested that proton degeneracy can play a significant role in the core [56], but find that protons are solely degenerate for r ≲ 4 km, which is outside the region of high muon abundance, hence proton degeneracy effects are neglected.
Given the above results and the non-relativistic velocities of the muons, the production rate for the Compton process is [57] 5 with x ≡ ω/m µ and ω the energy of the emitted axion.
For muon-proton bremsstrahlung, we have [57] µ ω e −ω/T F (w, y) where we have defined w ≡ ω/T and y ≡ k S / 2m µ T with k S the Debye screening scale, which is the appropriate scale to use for non-degenerate conditions, and x, t are integration variables.
There are two relevant regimes of the g aµ parameter space: the free-streaming regime (small g aµ ) and the trapping regime (large g aµ ).In the free-streaming regime, the mean-free path for axion absorption is significantly larger than the scale of the PNS, hence axions can escape the inner regions of the star without any further interactions, leading to volume emission.At sufficiently strong couplings, however, we enter the trapping regime.The mean free path becomes small in the interior of the PNS and axions are rapidly produced and reabsorbed out to some radius where the mean free path grows long and they escape.This radius is known as the "axion sphere" [58] and the luminosity can be approximated by blackbody emission from that radius at the local temperature [23].
Let us begin with the free-streaming regime.In order to compute the free-streaming luminosity of axions due to the muon coupling, we use the following [27]: 5 A word of caution: in a previous version of this paper, the following equation was missing a factor of two.This can be traced to Eq. (2.19) in Ref. [57], which does not include the sum over two photon polarizations.Ref. [57] explicitly states below Eq. (2.9) that these rates must be multiplied by the corresponding factor of two in order to yield the total production rate, however this step has been overlooked in much associated literature.Additionally, in a previous version, the term in brackets was taken in the limit ω ≪ mµ, which is not appropriate here as the typical photon energy is ω ∼ 3T ∼ mµ.The inclusion of this factor weakened the initially presented bounds by a factor of 5. A note addressing this appears as an Erratum in the published work.
with τ the optical depth 6 , given by × (e ω/T )(1 − e −ω/T ) , (7) where the factor of e ω/T appears since we are now considering the absorptive widths of the processes, which are related to the production rates by detailed balance: Γ prod = e −ω/T Γ abs , and the factor of (1 − e −ω/T ) appears as a result of using the reduced absorption rate (see Footnote 8).We include the factor of e −τ (ω,r) in the computation of the luminosity to account for the moderate reabsorption of axions by muons even in the "freestreaming" regime. 7 In the trapped regime, we must first identify the axion sphere, which is done by computing the reduced Rosseland mean opacity 8 for axions, given by where w ≡ ω/T .We can then integrate this radially to find the optical depth τ R : Whenever τ R ≫ 1, the axions are efficiently trapped [29].
The axion sphere is defined as the radius r a at which τ R (r a ) = 2 3 .The resulting luminosity is then just approximated as blackbody emission from that radius: As this approximation implicitly assumes that all axions decouple at the same radius, it is only applicable in regions where the coupling is strong enough that absorptive 6 The appropriate optical depth in this expression is, in contrast to what appeared in an earlier version of this work, calculated using the reduced absorptive widths (see Footnote 8).Additionally, the exponential factor in Eq. 6 should be averaged over trajectories, however, in the free-streaming regime, this factor is very near one and hence the differences are negligible. 7A very early version of this paper contained an additional factor of 2ω in the denominator of Eq. 6, which was in error.This caused a factor of one hundred suppression in our bounds.We thank the authors of Ref. [59] for bringing this to our attention. 8There are two different procedures being performed in the computation of the following equation.The first is the usual Rosseland mean of the opacity.The second is an explicit reduction of the absorption rate due to stimulated emission, which is necessary to compute the opacity in the context of radiative transfer.The reduction procedure introduces an additional factor of e w /(e w − 1) in the integrand, as the reduced absorption rate is given by Γ Compton abs × (1 − e −w ).Note as well that in a previous version of this paper, the equation below contained an erroneous additional factor of 1/2.However, due to a compensating factor of two in Eq. 3 (see Footnote 5), the result was unchanged.FIG.5: The free-streaming luminosity in axions as a function of axion-muon coupling.Couplings for which the curves lie above the dotted black line and to the right of the dashed are excluded by the free-streaming luminosity, while couplings to the left of the vertical dashed lines are excluded by the blackbody approximation in the trapped regime.The decrease of the free-streaming luminosity towards weak couplings is due to rapidly declining axion production, whereas the decrease towards strong couplings is due to an increasing absorptive width and corresponding suppression by the e −τ (ω,r) term.Note that the free-streaming luminosity is not the total luminosity, as it neglects the contribution from blackbody surface emission when the optical depth is large.It is shown to the left of the vertical lines purely to demonstrate the expected fall-off in the free-streaming (volume) emission and has hence been plotted as dot-dashed in that region.
processes shut off over a very narrow region.Otherwise, the energy dependence of the absorptive width stretches this single radius of decoupling into an extended decoupling region.At large values of the coupling, the approximation of an axion sphere is applicable.At the lower limit, as we will show below, the trapped regime transitions smoothly into the free-streaming regime while still producing luminosities in excess of the cooling constraint, hence the effect of an extended decoupling region does not affect our overall bound.

IV. RESULTS
Figure 5 shows the results of our luminosity calculations.The curves correspond to the free-streaming luminosity given in Eq. 6 as a function of coupling, with the horizontal axis showing − log 10 (g aµ /GeV −1 ).Note that towards weak couplings, the luminosity falls off as the coupling becomes too weak to allow the efficient production of axions.Towards stronger couplings, the freestreaming luminosity also falls off, now due to the increasing absorptive width (i.e. the factor of e −τ (ω,r) causes a large suppression).
The vertical lines demarcate our free-streaming regime from our trapped regime.We separate these two regimes at the coupling for which the optical depth within the axion sphere rises above 2/3.Note that the transition FIG.6: Our constraints, plotted on axion parameter space.Excluded regions are colored.We show the bound that is placed with the SFHo-18.8 profile as a solid blue line, which we quote as our ultimate bound.Note that this is extremely conservative, as model SFHo-18.8 is on the low edge of those compatible with SN 1987A.To illustrate the model dependence, we also show the result for SFHo-20.0 (light blue), which sits on the high edge of the allowed range.We additionally plot existing bounds from virtual axion contributions to the muon g−2 (purple), which are the most robust experimental constraints [20], and a rough estimate of CMB constraints (dashed line, see Supplemental Material for a thorough discussion).
occurs at couplings where we are still excluded by freestreaming, allowing a smooth transition between regimes and no breaks in our bounds.If, instead, the trapped regime had only begun at much stronger couplings, there could have in principle been a gap in the constraint between the two regimes.But this is far from the considered case, hence the actual details of the transition between regimes do not affect our overall results.Note that even at high couplings, the axion sphere is never further out than ≈ 18 km, as should be evident from the plot of the muon number density (Fig. 3).An axion sphere at 18 km corresponds to a blackbody luminosity of ≳ 10 53 erg/s, so our bounds extend smoothly to high couplings. 9  In summary, we can exclude all couplings to the left of the vertical lines (trapped regime) and can exclude all couplings to the right of the vertical line for which the free-streaming luminosity lies above the horizontal 9 At very high couplings, the axion sphere will be pushed to a radius at which the axion luminosity decreases beneath that of the neutrinos.However, the precise location of this transition is subject to considerable uncertainty as it is highly sensitive to the exponential fall-off of the muon abundance.A naive estimate comparing the blackbody luminosity in axions to the canonical Lν = 3 × 10 52 erg/s places this bound just above the g − 2 limit, however relying on the one-dimensional simulation results presented here is likely inappropriate, as 3D effects may significantly change the decoupling region.Future work is required to better elucidate the behavior in the very high coupling regime.
line corresponding to L a = 3 × 10 52 erg/s (free-streaming regime).These results can be reformulated into a constraint on the axion coupling, the result of which is displayed in Figure 6.Translated into a numerical bound, the constraint from SFHo-18.8 is g aµ < 10 −7.5 GeV −1 .Note that this is a conservative constraint, as it is placed with the coolest profile, a profile that is on the low-mass edge of the allowed range for SN1987A's remnant.As such, we also show the bound from SFHo-20.0 (lighter blue), which sits on the high-mass edge of the allowed range, though we only quote the conservative bound as our final result.
Additionally, we have cut off our bound at m a = 1 MeV as our production calculations implicitly took the axion to be effectively massless, a good approximation when m a ≪ T , as is the case here.More parameter space (up to m a ∼ 100 MeV) could be excluded with a full treatment of a massive axion, though the corresponding Boltzmann suppression will cause the bounds to rise rapidly towards higher masses (see Ref. [59]).
In conclusion, SN cooling constrains the axion-muon coupling to g aµ < 10 −7.5 GeV −1 (11) for axions with masses less than an MeV, severely limiting the parameter space for exotic axion-like particles with muon-dominated interactions.While this analysis was performed in the context of one specific model (an axion with muon-dominated interactions), we wish to note that the general insights presented here open up supernovae as a promising new avenue for robustly probing muonic BSM physics.(For one such extension of this work, see e.g.[59].) Appendix A: Supplemental Material 1.Comparison to bounds from N eff Limits on the axion-muon coupling from CMB measurements have also appeared in the literature [18,21,22].The physics behind these limits is that muons may produce axions in the early universe, increasing the effective number of relativistic species at the time of the CMB, pushing N eff above the measured value.
There is a very rapid turnoff in ∆N eff as a function of axion-muon coupling [21,22].This can be understood from the fact that once the coupling becomes sufficiently strong that the axions reach thermal equilibrium with the Standard Model bath, their contribution approaches a fixed value around ∆N eff ≈ 0.5.This rapid turn-off has a critical impact on computing bounds-it means that while CMB observations can place a strong bound on the axion-muon coupling if the uncertainty on the N eff measurement is ≲ 0.5, the bounds disappear essentially entirely if the uncertainty is ≳ 0.5.This introduces considerable uncertainty on the bound on the axion coupling when the bounds on ∆N eff are around 0.5.
The current precise bound on the axion-muon coupling found in the CMB literature from running a Boltzmann code is around the level of g aµ ≲ 10 −4 GeV −1 [22].However using the results from the same paper (Figure 2 of Ref. [22]) we roughly estimate that a significantly stronger bound could be placed [60].But at the same time we must mention that depending on what is believed about current error bars on ∆N eff , the bound may be significantly weakened, which may be the reason the authors of [22] chose to quote a much weaker bound.Either way, CMB S4 should clear the situation up by providing much smaller uncertainty on ∆N eff [21].Using Figure 2 of Ref. [22] and the value ∆N eff ≈ 0.34 from Ref. [61], we find g aµ ≈ 10 −7 GeV −1 for m a ≪ m µ .This should not be considered a strict limit, however, as a detailed analysis would be required in order to place a robust bound.We have chosen to plot it as a dashed line in Fig. 7 to indicate this.
Depending on the assumptions used when computing ∆N eff , the 2σ uncertainty ranges from 0.34 (Eq.67b of [61]) under the assumption of a pure ΛCDM model and BBN calculations for the helium fraction Y p to 0.58 when Y p floats independently of N eff and stellar measurements are not considered (Eq.82 of [61]) [60].This means that depending on the assumptions made in the calculation of N eff , the bound on the axion-muon coupling could either be strong or non-existent.The favored measurement of ∆N eff , using pure ΛCDM, possesses an uncertainty at 2σ of ∼ 0.34 and a central value slightly below zero, which would make ∆N eff = 0.5 a roughly ∼ 3.3σ deviation from the measured value.However, this deviation could be reduced by any modifications to a cosmological model beyond ΛCDM, e.g. to explain the ongoing Hubble tension, which shift the measurement of N eff upwards.As noted, this would have a significant effect on the bounds.
Additionally, it should be noted that independent of the above considerations, if the reheating temperature of the universe was below m µ , these bounds would be greatly weakened.Of course this is a small part of the currently allowed range for the reheating temperature.But if, for example, a confirmed positive signal was observed in laboratory experiments and CMB analysis indicated small ∆N eff , then we would likely not declare a contradiction but instead take this as interesting evidence for a low reheating scale.In that sense, this CMB bound has some model-dependence.
Taken as a whole, our SN bounds and our estimates of the CMB bounds appear roughly comparable.They are complementary in the sense that the uncertainties and model dependencies in each analysis are completely orthogonal to one another, hence this region of parameter space is strongly disfavored by entirely independent physics.

Plot of F deg
As is explained in the text, we quantify the effects of muon degeneracy on the Compton rate by computing the suppression factor F deg (see Sec. 3.2.5 of Ref. [23]): where n th µ is the thermal number density and f µ (E) = (e (E−µµ)/T + 1) −1 the energy distribution of muons at the given degeneracy and temperature.We have plotted this suppression factor in Figure 7. Due to the mild degeneracy, the suppression is never in excess of ∼ 15%.

3 )FIG. 3 :
FIG.3: Muon number density.Note that despite large differences in peak temperature, the muon number density does not change by more than an order of magnitude in the relevant region of high muon density (5 km ≲ r ≲ 15 km).

FIG. 4 :
FIG.4: Muon degeneracy parameter µµ/T for various profiles.This ratio never exceeds a small O(1) value in the relevant region of high muon density (5 km ≲ r ≲ 15 km).

TABLE I :
Model name Equation of state Progenitor mass (M⊙) NS bary.mass (M⊙) NS grav.mass (M⊙) Supernova model parameters and resulting NS baryonic and gravitational masses.We place our final bound with SFHo-18.8, as it produces the weakest constraint.