# Muonic Boson Limits: Supernova Redux

Andrea Caputo  School of Physics and Astronomy, Tel-Aviv University, Tel-Aviv
69978, Israel  Department of Particle Physics and Astrophysics, Weizmann
Institute of Science, Rehovot 7610001, Israel  Georg Raffelt  Max-Planck-
Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, 80805
München, Germany  Edoardo Vitagliano  Department of Physics and Astronomy,
University of California, Los Angeles, California 90095-1547, USA

(September 7, 2021)

###### Abstract

We derive supernova (SN) bounds on muon-philic bosons, taking advantage of the
recent emergence of muonic SN models. Our main innovations are to consider
scalars $\phi$ in addition to pseudoscalars $a$ and to include systematically
the generic two-photon coupling $G_{\gamma\gamma}$ implied by a muon triangle
loop. This interaction allows for Primakoff scattering and radiative boson
decays. The globular-cluster bound $G_{\gamma\gamma}<0.67\times
10^{-10}~{}{\rm GeV}^{-1}$ carries over to the muonic Yukawa couplings as
$g_{a}<3.1\times 10^{-9}$ and $g_{\phi}<4.6\times 10^{-9}$ for
$m_{a,\phi}\lesssim 100$ keV, so SN arguments become interesting mainly for
larger masses. If bosons escape freely from the SN core the main constraints
originate from SN 1987A $\gamma$ rays and the diffuse cosmic $\gamma$-ray
background. The latter allows at most $10^{-4}$ of a typical total SN energy
of $E_{\rm SN}\simeq 3\times 10^{53}$ erg to show up as $\gamma$ rays, for
$m_{a,\phi}\gtrsim 100$ keV implying $g_{a}\lesssim 0.9\times 10^{-10}$ and
$g_{\phi}\lesssim 0.4\times 10^{-10}$. In the trapping regime the bosons
emerge as quasi-thermal radiation from a region near the neutrino sphere and
match $L_{\nu}$ for $g_{a,\phi}\simeq 10^{-4}$. However, the $2\gamma$ decay
is so fast that all the energy is dumped into the surrounding progenitor-star
matter, whereas at most $10^{-2}E_{\rm SN}$ may show up in the explosion. To
suppress boson emission below this level we need yet larger couplings,
$g_{a}\gtrsim 2\times 10^{-3}$ and $g_{\phi}\gtrsim 4\times 10^{-3}$. Muonic
scalars can explain the muon magnetic-moment anomaly for $g_{\phi}\simeq
0.4\times 10^{-3}$, a value hard to reconcile with SN physics despite the
uncertainty of the explosion-energy bound. For generic axion-like particles,
this argument covers the “cosmological triangle” in the
$G_{a\gamma\gamma}$–$m_{a}$ parameter space.

††preprint: MPP-2021-154

###### Contents

  1. I Introduction
  2. II Muonic boson interactions
    1. II.1 Tree-level interaction structure
    2. II.2 Photo production
    3. II.3 Bremsstrahlung
    4. II.4 Two-Photon Coupling
    5. II.5 Primakoff process
    6. II.6 Bounds from Primakoff conversion
    7. II.7 Two-photon decay
    8. II.8 Two-photon coalescence
  3. III SN 1987A energy loss
    1. III.1 Free-streaming case
      1. III.1.1 Introductory remarks
      2. III.1.2 Application to muonic bosons
      3. III.1.3 Using the Garching SN models
      4. III.1.4 Emission rate
      5. III.1.5 Relativistic corrections
    2. III.2 Trapping limit
      1. III.2.1 Introductory remarks
      2. III.2.2 Simple bound from energy transfer
      3. III.2.3 Energy loss from the boson sphere
      4. III.2.4 Reduced absorption rate
      5. III.2.5 Spectral average
      6. III.2.6 Application to muonic bosons
      7. III.2.7 Volume emission vs. SB approximation
  4. IV Explosion energy
  5. V Decay photons from SN 1987A
    1. V.1 SMM observations
    2. V.2 Gamma signal from massive-particle decays
    3. V.3 ALP decays
    4. V.4 Muonic bosons
  6. VI Diffuse gamma-ray background
    1. VI.1 Redshift integral
    2. VI.2 Cosmic core-collapse rate
    3. VI.3 Generic limit
      1. VI.3.1 Short-lived bosons
      2. VI.3.2 Long-lived bosons
    4. VI.4 ALPs
    5. VI.5 Muonic bosons
      1. VI.5.1 Simplified estimates
      2. VI.5.2 Full numerical distributions
  7. VII From Colliders to Cosmology
    1. VII.1 Cosmology
    2. VII.2 Colliders
  8. VIII Axion-like particles
  9. IX Discussion and summary

##  I Introduction

Traditionally muons have been ignored in core-collapse supernova (SN)
simulations, although it is well known that neutron stars contain lots of
muons. Moreover, comparing the muon mass of $m_{\mu}=105.66~{}{\rm MeV}$ with
temperatures of 30–60 MeV in the hottest regions of collapsed SN cores reveals
that muons are not strongly suppressed. However, it is only recently that
muons and concomitant six-species neutrino transport was implemented in the
Garching group’s Prometheus Vertex code in a unique effort [1, 2]. While
directly after collapse the trapped electron lepton number provides for large
$e$ and $\nu_{e}$ chemical potentials, the core soon begins to deleptonize by
$\nu_{e}$ emission and to muonize by $\bar{\nu}_{\mu}$ losses. For
illustration we show in Fig. 1 the Garching model SFHo-18.8 at 1 s postbounce
(pb), the coldest of the models of Ref. [3]. In the critical region around 10
km with a typical temperature of 30 MeV, the muon density is around a quarter
that of electrons.

Figure 1: Profile of the Garching muonic SN model SFHo-18.8 at $t_{\rm pb}=1$
s [3] that we will use as our “cold reference model.” Top: Chemical
potentials, temperature, plasma frequency $\omega_{\rm p}$, and Debye
screening scale $k_{\rm s}$. Bottom: Number densities $n_{i}$, normalized to
$n_{0}=0.181~{}{\rm fm}^{-3}$, corresponding to nuclear density of $3\times
10^{14}~{}{\rm g}~{}{\rm cm}^{-3}$.

This large muon density invites one to extend the traditional SN 1987A
particle bounds [4, 5, 6] to bosons that couple specifically to muons in
violation of flavor universality. Otherwise more information would be expected
from interactions with first-generation fermions. One case in point that has
been studied in the SN context is a new $Z$ boson from a gauged
$L_{\mu}{-}L_{\tau}$ number [7]. For very small masses, such bosons can also
engender long-range muonic forces between neutron stars [8] or carry away
energy from binary pulsars [9].

We here focus on a new muon-philic scalar $\phi$ (muonic scalar for short)
that is motivated as one explanation for the observed discrepancy between the
measured and predicted muon magnetic moment [10] that now stands at a
$4.2\sigma$ significance [11, 12, 13, 14]. The required muonic Yukawa coupling
is [15]

$g_{\phi}\simeq 0.4\times 10^{-3},$ (1a) whereas a much larger value makes the
discrepancy worse. A pseudoscalar always makes it worse, implying an upper
bound [16] $g_{a}\lesssim 0.95\times 10^{-3}.$ (1b)

In both cases Yukawa couplings larger than the $10^{-3}$ level are excluded by
$g_{\mu}{-}2$ and thus marks the largest coupling worth studying with
astrophysical arguments.

Our work is inspired by a recent analysis concerning a muonic axion-like
pseudoscalar [3], but scalars may be more interesting in that they are
actually ruled in rather than ruled out by $g_{\mu}{-}2$ alone. In the SN
context, differences arise from the cross section of the main production
process $\gamma+\mu\to\mu+a$ or $\phi$ that we show explicitly in Fig. 3
below. The gist is that for the same Yukawa coupling, the pseudoscalar cross
section is always smaller. Deep in the SN core where photon energies are of
the same order as the muon mass, the reduction is around a factor of 3. In the
trapping limit where decoupling is in the neighborhood of the neutrino-sphere,
the relative factor is $(\omega/m_{\mu})^{2}$ and thus the pseudoscalar cross
section is relatively suppressed by a factor of around 10. In both cases we
find a gap between the $g_{\mu}{-}2$ inspired values and the SN 1987A cooling
argument.

However, these specific results depend on the exponentially decreasing muon
abundance at the proto-neutron star (PNS) surface in our 1D reference models.
A very different picture may arise in more realistic 3D models where the
material in this region can be subject to strong motions, redistributing the
muons.111We thank Thomas Janka for stressing this point in a private
communication. In our case this question becomes moot once we include
Primakoff scattering as the main opacity source.

Figure 2: Constraints on the muonic Yukawa coupling $g_{a}$ of pseudoscalars
(left panel) and $g_{\phi}$ of scalars (right panel) as a function of boson
mass. We also show the constraint where the muon $g_{\mu}{-}2$ discrepancy
would get worse (pseudoscalars) or could be explained (scalars). The shading
or double arrows indicate the excluded range, except for the scalar range
ruled in by $g_{\mu}{-}2$. The SN bounds are derived in this paper and
summarized in Table 4 using both a hot and a cold SN reference model. We here
show conservative limits (the larger number in the free-streaming case and the
smaller one in the trapping regime). The HB-star bounds are taken from Ref.
[17], the cosmological $N_{\rm eff}$ constraint from Ref. [18].

This crucial effect comes about because actually our main innovation is to
include systematically the generic effective two-photon interaction. A
(pseudo)scalar, even if it couples exclusively to muons, at one loop
inevitably obtains a two-photon vertex that allows for $\phi$–$\gamma$
Primakoff conversion and for the decay $\phi\to 2\gamma$. Unless this
$2\gamma$ coupling is fine-tuned to disappear in a UV complete theory, it
dominates our arguments and especially the $2\gamma$ decays provide the most
restrictive SN limits. So we do not consider (pseudo)vectors because their
$2\gamma$ decay is forbidden by the Landau-Yang Theorem.

In accordance with the previous literature we normalize the two-photon
couplings $G_{a\gamma\gamma}$ and $G_{\phi\gamma\gamma}$ such that
pseudo(scalars) have the same Primakoff cross section and the same decay rate.
Conversely this means that for a given limit on
$G_{\gamma\gamma}=G_{a\gamma\gamma}=G_{\phi\gamma\gamma}$, notably from the
CAST experiment [19] and the helium-burning lifetime of horizontal-branch (HB)
stars [20, 21, 17], translates to a bound on the underlying $g_{\phi}$ that is
a factor of 2/3 less restrictive than on $g_{a}$ as can be seen in our summary
plot Fig. 2. Here and in the following, the differences between scalars and
pseudoscalars are never dramatic, yet always warrant separate treatments.

The loop-induced two-photon vertex of muonic pseudoscalars was also mentioned
by Croon et al. in their discussion of SN limits [7]. However, they assumed an
axial-vector derivative interaction structure of the form
$(g_{a}/2m_{\mu})\partial^{\nu}a\,\bar{\mu}\gamma_{\nu}\gamma_{5}\mu$ instead
of the pseudoscalar form $-ig_{a}a\,\bar{\mu}\gamma_{5}\mu$. These two
possibilities are often equivalent [22], but provide different loop factors in
the present situation of virtual muons. For massless pseudoscalars, the loop
contribution vanishes for the derivative structure, while it is unsuppressed
for the pseudoscalar case. Which form is appropriate depends on the UV
completion of the theory and the possible Nambu-Goldstone nature of the
pseudoscalar. We will always use the pseudoscalar coupling because in this
case the radiative decay is not suppressed and dominates our arguments.
However, if one assumes the opposite one can always fall back on our “tree-
level only” constraints that are explicitly listed, for example, in our
summary Table 4 below and that broadly agree with the earlier literature. For
scalars there is no such ambiguity.

For large couplings, on the SN trapping side, the pseudoscalar opacity in the
decoupling region is dominated by Primakoff scattering on charged particles
and comparable for scalars. As a consequence, the SN 1987A energy-loss limits
are essentially the same for scalars and pseudoscalars shown in Fig. 2 and not
even close to the $g_{\mu}{-}2$ inspired values.

However, the most dramatic consequence of the two-photon vertex is the
radiative decay $\phi,a\to 2\gamma$ even though it is strongly phase-space
suppressed for low-mass bosons. However, the HB-star bounds on
$G_{\gamma\gamma}$ pertain up to masses of around 100 keV, so the SN arguments
are anyway interesting mainly for larger masses. In this case the $2\gamma$
decay becomes so fast that it is a dominating effect and far more important
than Primakoff scattering for most of our arguments.

The effect of decays is particularly dramatic on the trapping side. If the
boson luminosity is comparable to that of neutrinos, essentially corresponding
to the SN 1987A energy-loss limit, the entire boson energy is dumped into the
surrounding matter of the progenitor star and makes SN explosions far too
energetic.222We call this the Falk-Schramm argument because in 1978 these
authors used similar reasoning to constrain radiative decays of neutrinos
[23]. The visible energy of a core-collapse SN explosion of 1–$2\times
10^{51}\,{\rm erg}$ is less than 1% of the binding energy of the final neutron
star of $E_{\rm SN}=2$–$4\times 10^{53}~{}{\rm erg}$, most of which is
normally carried away by neutrinos. Therefore, the energy carried away by
bosons must be much less than what is allowed by SN 1987A neutrinos. We find
that this argument actually closes the gap to the $g_{\mu}{-}2$ level as shown
in Fig. 2 if we demand that bosons carry less than 1% of the total, but even
10% would be enough to close the gap. This argument also covers what has been
called the “cosmological triangle” for generic axion-like particles (ALPs) as
can be seen in Fig. 12.

These and the following discussions heavily use the Garching muonic SN models
[3], using their SFHo-18.8 model with an inner $T$ of around 30 MeV as our
cold reference case and LS220-20.0 as our hot one with around twice the inner
$T$. We derive nominal bounds by post-processing these models and show them
separately for the hot and cold models, for example in our summary Table 4.
This approach ignores the feedback of the new particles on SN physics, but
probably gives us a good sense of the relevant parameter range.

These remarks do not necessarily pertain to the Falk-Schramm argument because
to suppress the emission to $0.01\,L_{\nu}$ means that the bosons derive from
a region at significantly larger radii than the neutrino sphere. Post-
processing an unperturbed model could yield an unreliable answer in a region
where the atmospheric structure must depend on boson energy transfer. In this
sense we are least confident of the excluded region “Explosion energy” in Fig.
2. On the other hand, the class of subluminous Type-II Plateau SNe has much
smaller explosion energies, perhaps as low as $10^{50}~{}{\rm erg}$ or less
[24, 25, 26, 27, 28, 29]. Our reference models do not necessarily provide good
proxies for this class. Either way, the explosion-energy argument may deserve
a dedicated investigation of energy transfer by ALPs in this SN region.

Our arguments become much more straightforward on the feeble-interaction side
of the exclusion range where the new particles stream freely from the SN core
once produced. The boson luminosity turns on slowly after collapse in sync
with the inner core heating up and a significant muon population building up.
So the feedback on SN explosion physics is minimal.

For SN 1987A, the boson decay photons, not neutrinos, provide the most
sensitive probe. The 100-MeV-range $\gamma$ rays from bosons decaying between
SN 1987A and Earth would have been picked up by the Gamma-Ray Spectrometer
(GRS) on board the Solar Maximum Mission (SMM) satellite that was operational
at the time. A long time ago, these data were used to constrain neutrino
radiative decays [30, 31, 32] and very recently ALPs, particles that interact
only by their two-photon vertex [33]. Our case is analogous, except that here
photo production on muons dominates, not Primakoff production.

This argument relies on the GRS $\gamma$-ray signal integrated over a time
window of 223.2 s after the neutrino burst and as such does not depend on the
detailed time structure of the putative $\gamma$-ray signal, only on the total
SN energy emitted in the form of bosons and their typical energies. With this
information, the constraint follows from a simple analytic expression derived
e.g. in Ref. [32]. The constraint on the boson coupling strength here requires
taking a fourth root relative to the limiting $\gamma$ fluence because the
coupling strength enters quadratically both at production and decay, so these
constraints are particularly forgiving of uncertainties of the assumed SN
1987A model. On the other hand, boson emission is here indeed a perturbative
effect, so in principle the unperturbed neutrino signal might discriminate
between possible SN 1987A models and make the boson emission characteristics
more specific.

For most masses, however, the boson decay photons from all past SNe to the
diffuse cosmic $\gamma$-ray background are yet more constraining. Recently
this approach was used for ALPs that are emitted from a SN core by different
processes [34]. We investigate the dependence of this argument on the cosmic
core-collapse rate and redshift distribution as well as the boson emission
characteristics of an average SN. We find that the redshift distribution
impacts the constraint only on the 10% level, so the only relevant information
is the total number of past SNe per comoving volume of $n_{\rm
cc}=0.6$–$1\times 10^{7}~{}{\rm Mpc}^{-3}$. We thus find that an average SN is
allowed to emit at most $10^{-4}$ of the neutron-star binding energy in the
form of radiatively unstable bosons, practically independent of their exact
spectral distribution. In this form the limit applies if most bosons decay
within a Hubble time, which in our context applies for $m_{a,\phi}\gtrsim 100$
keV. For smaller masses the constraint deteriorates as can be seen in Fig. 2,
yet continues to dominate the SN 1987A gamma-ray limit.

The diffuse $\gamma$-ray limit is far more constraining than the one from the
explosion energy, so it can be avoided only if the bosons decay within the
progenitor-star radius. By the same token, the explosion-energy argument
cannot be avoided by decays outside of the progenitor star, so both arguments
can be avoided only by producing fewer bosons.

The main part of our paper is devoted to working out these arguments in detail
and comparing with recent results of other authors. In Sec. II we discuss the
interaction structure, the two-photon coupling, and production processes for
(pseudo)scalar muonic bosons. In Sec. III we turn to the traditional SN 1987A
energy-loss argument, paying particular attention to the trapping regime to
clear up some confusion that has crept into the recent literature. This is
followed in Sec. IV by a discussion of the Falk-Schramm argument of excessive
energy deposition by decaying particles that would render SN explosions too
energetic. We then turn to the feeble-interaction side of our exclusion plot
and first consider in Sec. V decay photons from SN 1987A that would have been
picked up by the Gamma-Ray Spectrometer (GRS) on board the Solar Maximum
Mission (SMM) satellite. This discussion is followed in Sec. VI by our most
constraining argument, the contribution of decay photons to the cosmic diffuse
$\gamma$-ray background in the sub-100-MeV range from all past SNe. In Sec.
VII we briefly summarize some pertinent limits from cosmology and colliders.
Our arguments are often very similar to those for generic ALPs, so we
summarize our results in Sec. VIII specifically for this case, facilitating a
direct comparison with the earlier literature. The final Sec. IX is given over
to a discussion and summary of our findings. Our constraints are summarized in
Figs. 2 and 12 and Table 4.

Figure 3: Cross section for the muonic Compton process with a final-state
vector (solid), scalar (dashed) or pseudoscalar (dotted). The vector case
requires a factor of 2 for the final-state polarizations. The energy $\omega$
is considered in the muon rest frame. For $\omega\gtrsim m_{\mu}$ the scalar
and pseudoscalar cross sections quickly approach each other, while asymptotic
agreement with the vector case requires very large energies.

##  II Muonic boson interactions

###  II.1 Tree-level interaction structure

The starting point for our discussion is the assumption of scalars $\phi$ or
pseudoscalars $a$ that couple to muons through the Yukawa operators

$\mathcal{L}_{\phi}\supset-
g_{\phi}\phi\bar{\mu}\mu\quad\hbox{and}\quad\mathcal{L}_{a}\supset-
ig_{a}a\bar{\mu}\gamma_{5}\mu.$ (2)

For completeness sometimes we will also consider the vector coupling
$\mathcal{L}_{Z}\supset g_{Z}Z_{\mu}\bar{\mu}\gamma^{\mu}\mu$, with $Z$ being
a new massive vector boson coupled to muons only.

No other tree-level interactions are assumed to exist. For tree-level
processes, the pseudoscalar interaction is often equivalent to the axial-
vector derivative structure [22]

$\mathcal{L}_{a}\supset\frac{g_{a}}{2m_{\mu}}\,\partial^{\nu}a\,\bar{\mu}\gamma_{\nu}\gamma_{5}\mu.$
(3)

We will see, however, that the effective two-photon vertex through a muon
triangle loop is strongly suppressed in the derivative case. Notice also that
in the scalar case there is no equivalent derivative structure of the form
$\partial^{\nu}\phi\bar{\mu}\gamma_{\nu}\mu$ because after a partial
integration it is equivalent to a total derivative of the conserved muon
vector current, i.e., such a structure yields a vanishing matrix element in a
scattering process.

In the following we always ignore both the boson mass and the photon plasma
mass, which for the parameters and conditions of interest give negligible
modifications.

###  II.2 Photo production

The main boson production process in the inner SN core where muons abound is
photo production of the form $\gamma+\mu\to\mu+\phi$. The nonrelativistic
cross section for this “semi-Compton process” is the Thomson cross section

$\sigma_{\rm T}=\frac{\alpha g_{\phi}^{2}}{3m_{\mu}^{2}}.$ (4)

It is understood as the cross section of one unpolarized photon to produce
$\phi$. The reverse process of $\phi$ absorption sports a factor of 2 to count
two possible final-state photon polarizations.333This photo-production cross
section is the polarization-averaged cross section of a single photon. The
scalar or axion emission rate from a medium thus involves the number density
of photons summed over both polarizations. The often-cited photo production
rate (for the case of pseudoscalars) was given in Eq. (2.19) of Ref. [35] with
the additional explanation after Eq. (2.9) that the total rate requires a
factor of 2 to count both photon polarizations. Apparently this proviso was
frequently overlooked, leading to a missing factor of 2, for example, in Eq.
(5.2) of Ref. [7], Eq. (6) of the original version of Ref. [3], and Eq. (9)
and the subsequent inline equations of Ref. [36].

The usual Thomson cross section for photon scattering, $\sigma_{\rm
T}=8\pi\alpha^{2}/3m_{\mu}^{2}$, arises with $g_{\phi}\to e$,
$\alpha=e^{2}/4\pi$, and a factor of 2 for the two final-state polarizations,
i.e., the initial-state polarizations are averaged, the final-state ones
summed as usual.

For pseudoscalars we substitute $g_{\phi}\to g_{a}$ and need an additional
factor $(\omega/m_{\mu})^{2}$ for photon energy $\omega$ that arises from the
spin-dependent nature of the low-energy pseudoscalar coupling, so the cross
section is $\sigma_{\rm T}(\omega/m_{\mu})^{2}$ at low energies.444Notice that
this expression applies only for $\omega/m_{\mu}\ll 1$, whereas in a SN core
with $T\gtrsim 30~{}{\rm MeV}$, a typical $\omega\simeq 3T$ is similar to
$m_{\mu}$. The nonrelativistic expansion at such large energies [3, 7] vastly
overestimates the cross section as acknowledged in the updated version of Ref.
[3]. In the trapping limit the relevant conditions are those where the bosons
decouple with much smaller $T$, so in this case the low-energy expansion is
appropriate.

For general kinematics when $\omega$ is not small relative to $m_{\mu}$, the
semi-Compton cross section for the production of (pseudo)scalars and photons
is [4, 37, 38]

$\displaystyle\kern-20.00003pt\sigma_{\phi}$ $\displaystyle=$
$\displaystyle\sigma_{\rm
T}\frac{3}{4}\left[\frac{1-3\hat{s}}{2\hat{s}^{2}}-\frac{16}{(\hat{s}-1)^{2}}+\frac{(\hat{s}+3)^{2}}{(\hat{s}-1)^{3}}\ln\hat{s}\right],$
(5a) $\displaystyle\kern-20.00003pt\sigma_{a}$ $\displaystyle=$
$\displaystyle\sigma_{\rm
T}\frac{3}{4}\left[\frac{1-3\hat{s}}{2\hat{s}^{2}}+\frac{\ln\hat{s}}{\hat{s}-1}\right],$
(5b) $\displaystyle\kern-20.00003pt\sigma_{\gamma}$ $\displaystyle=$
$\displaystyle\sigma_{\rm
T}\frac{3}{4}\left[\frac{16}{(\hat{s}-1)^{2}}+\frac{\hat{s}+1}{\hat{s}^{2}}+\frac{2(\hat{s}^{2}-6\hat{s}-3)}{(\hat{s}-1)^{3}}\ln\hat{s}\right],$
(5c)

with $\hat{s}\equiv s/m_{\mu}^{2}$ and $\sqrt{s}$ the center of mass energy.
In Fig. 3 we show these expressions as a function of photon energy $\omega$ in
the rest frame of the muon where $s=m_{\mu}^{2}+2m_{\mu}\omega$. At large and
small energies the scalar cross section is half the cross section of the
photon case, but for large energies the asymptotic convergence is very slow.
We also notice that for $\omega\gg m_{\mu}$ the pseudoscalar cross section is
identical with the scalar case, whereas for $\omega\ll m_{\mu}$ it is
suppressed according to $\sigma_{a}=\sigma_{\rm T}(\omega/m_{\mu})^{2}$ as
mentioned earlier.

###  II.3 Bremsstrahlung

The bremsstrahlung process $\mu+p\to p+\mu+\phi$ is another potential particle
source. The corresponding axion emission rate for nonrelativistic electrons
was calculated in Ref. [39]. We follow the steps in that paper and note that
in the squared matrix element we need to substitute $\omega_{a}^{2}\to
2m_{e}^{2}$ to go from pseudoscalar to scalar emission. This factor follows in
analogy to the transition between pseudoscalar and photon emission in free-
free, free-bound and bound-bound transitions that was discussed in Refs. [40,
35]. Ignoring screening effects, we thus find for the energy-loss rate from
$\mu+p\to p+\mu+\phi$

$\displaystyle\epsilon_{\phi}$ $\displaystyle=$ $\displaystyle
g_{\phi}^{2}\,Y_{\mu}Y_{p}\frac{n_{B}}{m_{\mu}m_{p}}\frac{4\alpha^{2}}{3\pi}\sqrt{\frac{2}{\pi}}\,\left(\frac{T}{m_{\mu}}\right)^{1/2}$
(6) $\displaystyle=$ $\displaystyle
g_{\phi}^{2}Y_{\mu}Y_{p}\,\frac{n_{B}}{n_{0}}\,\left(\frac{T}{30~{}{\rm
MeV}}\right)^{1/2}1.86\times 10^{38}~{}\frac{{\rm erg}}{{\rm g}\,{\rm s}},$

where $n_{0}=0.181~{}{\rm fm}^{-3}$ is the baryon density corresponding to
$3\times 10^{14}~{}{\rm g}/{\rm cm}^{3}$ (nuclear density).555The scalar
bremsstrahlung emission rate was also calculated in Ref. [41] and the result
was given in terms of a multi-dimensional phase-space integral over angular
coordinates. Evaluating it numerically for our conditions we find that our
result is larger by a factor that is $\sqrt{2}$ within numerical accuracy, but
as Ref. [41] does not document enough details we cannot pin-point the origin
of the discrepancy.

For our reference conditions, bremsstrahlung is about an order of magnitude
smaller than the Compton rate. In addition, bremsstrahlung is somewhat
suppressed by screening effects. On the other hand, the nonrelativistic
approximation is not very good in either case, so the exact ratio is somewhat
uncertain. We neglect bremsstrahlung, but if it were to contribute some amount
of additional emission, our final bounds err in the conservative direction by
neglecting it. In the trapping limit, what matters are interactions in the
neutrino-sphere region where bremsstrahlung for sure is negligible.

###  II.4 Two-Photon Coupling

A (pseudo)scalar that couples to muons according to Eq. (2) inevitably also
has a effective two-photon coupling through a triangle loop. The effective
coupling can be written in the form

$\displaystyle\mathcal{L}_{\phi\gamma\gamma}$ $\displaystyle=$ $\displaystyle
G_{\phi\gamma\gamma}\phi\,({\bf E}^{2}-{\bf B}^{2}),$ (7a)
$\displaystyle\mathcal{L}_{a\gamma\gamma}$ $\displaystyle=$ $\displaystyle
G_{a\gamma\gamma}a\,{\bf E}\cdot{\bf B},$ (7b)

with the effective couplings

$\displaystyle G_{\phi\gamma\gamma}$ $\displaystyle=$
$\displaystyle\frac{2\alpha}{3\pi}\,\frac{g_{\phi}}{m_{\mu}}\,B_{a}(x),$ (8a)
$\displaystyle G_{a\gamma\gamma}$ $\displaystyle=$
$\displaystyle\frac{\alpha}{\pi}\,\frac{g_{a}}{m_{\mu}}\,B_{\phi}(x).$ (8b)

The loop factors $B$ depend on $x=m_{a,\phi}/2m_{\mu}$. Assuming $x<1$ so that
the bosons cannot decay into a muon pair, the loop factors are [42, 43]

pseudoscalar: $\displaystyle B_{a}(x)$ $\displaystyle=$
$\displaystyle\frac{{\rm ArcSin}^{2}(x)}{x^{2}}=1+\frac{x^{2}}{3}+{\cal
O}(x^{4}),$ (9a) derivative: $\displaystyle B_{a}(x)$ $\displaystyle=$
$\displaystyle\frac{{\rm ArcSin}^{2}(x)}{x^{2}}-1=\frac{x^{2}}{3}+{\cal
O}(x^{4}),$ (9b) scalar: $\displaystyle B_{\phi}(x)$ $\displaystyle=$
$\displaystyle\frac{3}{2x^{4}}\left[x^{2}-(1-x^{2}){\rm ArcSin}^{2}(x)\right]$
(9c) $\displaystyle=$ $\displaystyle 1+\frac{7x^{2}}{30}+{\cal O}(x^{4}).$

We consider bosons with $m\lesssim 10$ MeV so that $x^{2}=(m/2m_{\mu})\lesssim
2\times 10^{-3}$. Therefore, the loop factors can be safely neglected except
in the derivative case where the coupling is strongly suppressed by
$x^{2}/3\lesssim 0.7\times 10^{-3}$. So, for pseudoscalars, the two-photon
vertex depends on the assumed interaction structure, i.e., the axial-vector
derivative vs. the pseudoscalar one which are neither identical nor always
equivalent as explained earlier.

For boson decay, in the rest-frame of the decaying particles and with
$m_{a,\phi}\ll m_{\mu}$ the muon mass is largest scale. In Primakoff
scattering, the external photon energy can be around $m_{\mu}$ in the deep
interior where a typical $\omega=3T$ and $T$ can be as large as 30–50 MeV.
While a possible modification of the low-energy expansion will not be large,
it occurs in the deep interior where the semi-Compton process dominates. Of
course, for the latter one also needs to avoid the low-energy expansion. In
our context, the Primakoff effect dominates in the decoupling region where the
approximation $\omega\ll\operatorname{m}_{\mu}$ holds.

###  II.5 Primakoff process

The two-photon coupling allows, for example, for the Primakoff conversion
between (pseudo)scalars and photons in external electric or magnetic fields.
The way we have written the interaction structure, the Primakoff amplitude on
a charged particle comes either from ${\bf E}^{2}$ (scalar) or ${\bf
E}\cdot{\bf B}$ (pseudoscalar) and thus leads to the same cross section for
the same value of $G_{\gamma\gamma}$ which now stands for either
$G_{a\gamma\gamma}$ or $G_{\phi\gamma\gamma}$.

The Primakoff cross section for $\gamma+Ze\to Ze+a$ or $\phi$ on a
nonrelativistic charged particle with charge $Ze$, averaged over photon
polarizations, is

$\sigma_{Z}=\frac{Z^{2}\alpha G_{\gamma\gamma}^{2}}{2}\,f_{\rm s},$ (10)

where the screening factor is [39]

$f_{\rm s}=\frac{1}{4}\left[\left(1+\frac{k_{\rm
s}^{2}}{4\omega^{2}}\right)\log\left(1+\frac{4\omega^{2}}{k_{\rm
s}^{2}}\right)-1\right].$ (11)

In the Debye approximation this is

$k_{\rm
s}^{2}=\frac{4\pi\alpha\hat{n}}{T}\quad\hbox{where}\quad\hat{n}=\sum_{j}Z_{j}^{2}n_{j}\equiv\hat{Y}n_{B},$
(12)

which defines the effective charge $\hat{Y}$ per baryon. As we neglect the
boson mass and photon plasma mass, the cross section diverges logarithmically
in the forward direction, an effect that we control with Debye screening.

In this form the rate was derived for a nonrelativistic and nondegenerate
stellar plasma. For relativistic electrons, the Primakoff cross section should
not be very different, but electrons are degenerate in all regions of interest
in the SN core. So their contribution both as scattering targets and for
screening is significantly smaller and we neglect them. We also neglect muons
because in those regions where they abound the Compton process is much more
important.

One often thinks of the SN medium to consist primarily of neutrons and
protons, so the main Primakoff targets would be protons. However, in the
neutrino decoupling region of the muonic SN models there also appear small
nuclear clusters or small nuclei. In Fig. 4 we show the profiles of charged-
particle abundances in our reference model, where $Y_{j}$ is the number
fraction per baryon. The difference between $Y_{e}+Y_{\mu}$ and $Y_{p}$
highlights the appearance of these structures.

Figure 4: Top: Charged-particle abundances in our reference model shown in
Fig. 1. The difference between $Y_{e}+Y_{\mu}$ and $Y_{p}$ highlights the
appearance of light nuclei or nuclear clusters in the decoupling region.
Bottom: Average screening factor $f_{\rm s}$ for Primakoff scattering defined
in Eq. (11). The red line is the Rosseland average, the blue line a thermal
average for boson emission as explained in the text.

If we had abundance profiles for the individual clusters $j$ of charge $Z_{j}$
we could simply use their abundances, but what is explicitly listed are
profiles for light and heavy clusters as well as alpha particles. Notice that
what is listed are mass fractions $X_{j}$ which are constrained by
$\sum_{j}X_{j}=1$. If all nuclei were alpha particles, we would have
$1=X_{p}+X_{n}+X_{\alpha}$ and $\hat{Y}=Y_{p}+4X_{\alpha}/4=Y_{p}+X_{\alpha}$.
For protons and neutrons, $Y$ and $X$ is the same, so
$X_{\alpha}=1-Y_{n}-Y_{p}$ and finally

$\hat{Y}=1-Y_{n}.$ (13)

In reality, some of the clusters have smaller or larger $Z/A$ ratios than
$\alpha$ particles, so this is an approximation, but perhaps not worse than
neglecting electrons. Another approximation is to neglect degeneracy effects
of the charged particles that should be small wherever the Primakoff process
is important, but in any case is yet another small modification of the
effective screening scale $k_{\rm s}$. The benefits of a more refined
treatment are limited because of the overall uncertainties of the SN
arguments.

For our reference model, the screening scale is shown in the bottom panel of
Fig. 1. For energy losses in the deep interior, relevant for the energy-loss
argument in the free-streaming limit, we need the average Primakoff scattering
rate for a thermal distribution of initial-state photons, and an additional
factor of $\omega$ because we actually need the energy-loss rate. In this case
what we call thermal average is computed as

$\langle f_{\rm s}\rangle_{\rm
T}=\frac{15}{\pi^{4}}\int_{0}^{\infty}\\!\\!dx\,f_{\rm s}(x,k_{\rm
s}/T)\,\frac{x^{3}}{e^{x}-1},$ (14)

where we use $\omega=xT$. We show $\langle fs\rangle_{\rm T}$ as a blue line
in the bottom panel of Fig. 4.

For our purposes, Primakoff scattering is mostly important in the trapping
case where bosons decouple in the 16–18 km region where $k_{\rm s}\simeq T$.
To identify the “axiosphere” we use the Rosseland average of the mean free
path (MFP) as explained in Sec. III.2.5 with the weight function defined in
Eq. (37). Notice that $\langle f_{\rm s}\rangle_{\rm R}$ is calculated by
integrating $f_{\rm s}^{-1}$ over the normalized weight function and taking
the inverse afterward because we need the average MFP, not the average
interaction rate. We show $\langle f_{\rm s}\rangle_{\rm R}$ as a function of
radius as a red line in the bottom panel of Fig. 4.

Notice that $f_{\rm s}(\omega,k_{\rm s})=f_{\rm s}(x,k_{\rm s}/T)$ depends
only on $x$ and the dimensionless ratio $k_{\rm s}/T$. The two averages are
fortuitously very similar in the region where $k_{\rm s}\simeq T$ because of
the weak $\omega$ dependence and do not change much as a function of radius.
For simple estimates on the 20% precision level one could use something like
$f_{\rm s}=0.8$ everywhere, in particular as the overall precision of the
Primakoff cross section is probably not better than this because of the
discussed uncertainties of the chemical composition, electron contribution,
degeneracy effects, and concomitant screening prescription.

###  II.6 Bounds from Primakoff conversion

The two-photon vertex is the main interaction channel to search for axions,
leading to many constraints and ongoing and future projects [44, 22]. One
intriguing approach, first proposed by Sikivie [45], is to consider axion-
photon conversion in large-scale external magnetic fields in analogy to
neutrino flavor oscillations [46]. This method is particularly powerful for
smaller axion masses than we consider here, so we mention explicitly only the
conversion of axion-like particles emitted by SN 1987A in the galactic $B$
field, leading to $G_{a\gamma\gamma}<5.3\times 10^{-12}~{}{\rm GeV}^{-1}$ for
$m_{a}<4.4\times 10^{-10}~{}{\rm eV}$ [47, 48].

The solar axion search by the CAST experiment has established the constraint
[19]

$G_{a\gamma\gamma}<0.66\times 10^{-10}~{}{\rm GeV}^{-1}\quad\hbox{(95\%
C.L.)}$ (15)

that has become a reference value for this coupling, but only applies for
$m_{a}\lesssim 0.2~{}{\rm eV}$. Another long-standing constraint derives from
the energy loss of horizontal-branch (HB) stars that would reduce their
lifetime. As a result, fewer HB stars would be observable in globular clusters
relative to other phases of evolution [20]. A recent update happens to be
identical with Eq. (15) at two significant digits [21], but extends to masses
up to some 10 keV, corresponding to the internal HB-star temperature. For
larger masses, the $e^{-m_{a}/T}$ suppression kicks in only slowly with
increasing mass, but then drops sharply at $m_{a}\gtrsim 200~{}{\rm keV}$
[17].666We use the curve from Fig. 4 of Ref. [17]. However notice that the SN
bound (green region) in that figure which seems to cover the entire globular-
cluster excluded region is incorrect and actually in contradiction with e.g.
Fig. 16 of Ref. [49] by some of the same authors. So the globular-cluster
argument indeed covers a range of parameters not excluded by SN arguments.

As reference limits for the Yukawa couplings of our muonic (pseudo)scalar
bosons we translate Eq. (15) according to Eq. (8) and find the 95% C.L. limits

$\displaystyle g_{\phi}$ $\displaystyle<$ $\displaystyle 4.6\times 10^{-9},$
(16a) $\displaystyle g_{a}$ $\displaystyle<$ $\displaystyle 3.1\times
10^{-9}.$ (16b)

We show them in Fig. 2, where the mass dependence follows from Ref. [17].

###  II.7 Two-photon decay

In addition, for a nonvanishing (pseudo)scalar mass, the two-photon decay
$\phi$ or $a\to 2\gamma$ is possible and occurs with the rate

$\Gamma_{\gamma\gamma}=\frac{G_{\gamma\gamma}^{2}m_{\phi}^{3}}{64\pi}$ (17)

for both the scalar and pseudoscalar case. In terms of the Yukawa couplings,
the decay rates are

$\Gamma_{\phi\gamma\gamma}=\left(\frac{2}{3}\right)^{2}\frac{g_{\phi}^{2}\alpha^{2}}{64\pi^{3}}\frac{m_{\phi}^{3}}{m_{\mu}^{2}},\quad\hbox{and}\quad\Gamma_{a\gamma\gamma}=\frac{g_{a}^{2}\alpha^{2}}{64\pi^{3}}\frac{m_{a}^{3}}{m_{\mu}^{2}}.$
(18)

Numerically this is in the scalar case

$\Gamma_{\phi\gamma\gamma}=260~{}{\rm s}^{-1}\,\left(\frac{g_{\phi}}{0.4\times
10^{-3}}\right)^{2}\,\left(\frac{m_{\phi}}{{\rm MeV}}\right)^{3},$ (19)

where the reference coupling strength corresponds to the approximate value
required to explain the muon magnetic moment anomaly. It is allowed by the SN
1987A energy-loss argument.

The lab-frame decay rate involves a Lorentz factor $m_{\phi}/E_{\phi}$, so the
MFP against decay is

$\lambda^{\rm lab}_{\phi\gamma\gamma}=1.2\times 10^{9}\,{\rm
cm}\,\left(\frac{0.4\times 10^{-3}}{g_{\phi}}\right)^{2}\,\left(\frac{{\rm
MeV}}{m_{\phi}}\right)^{4}\frac{E_{\phi}}{10\,{\rm MeV}}.$ (20)

The reference energy of $E_{\phi}=10$ MeV is characteristic of bosons emitted
from the neutrino-sphere region. The smallest mass before entering the HB-star
exclusion range is perhaps 0.2 MeV, implying $\lambda^{\rm
lab}_{a\gamma\gamma}\simeq 0.7\times 10^{12}~{}{\rm cm}$, within the envelope
of the progenitor of SN 1987A.

###  II.8 Two-photon coalescence

The reverse process of photon coalescence $2\gamma\to\phi$ is also possible.
Of course the in-medium effective photon mass must be small enough compared
with the boson mass. Lucente et al. [49] studied SN limits on ALPs and found
that photon coalescence as an emission process is negligible compared to
Primakoff scattering except for boson masses beyond a few 10 MeV. As we will
restrict our discussion to boson masses up to 10 MeV we may ignore this
process.

##  III SN 1987A energy loss

###  III.1 Free-streaming case

####  III.1.1 Introductory remarks

Shortly after the observation of the SN 1987A neutrino signal it became clear
that the duration of several seconds and the observed energy is incompatible
with excessive energy loss in hypothetical new forms of radiation such as
axions [50, 51, 52, 53]. After the explosion, probably some hundreds of ms
after collapse, the evolution of the remaining proto neutron star (PNS) is
deleptonization and cooling on a diffusion time scale of a few seconds because
neutrinos are trapped [54, 55, 56]. In this way one probes the interaction
strength of very feebly interacting particles that escape freely from the SN
core.

From a modern perspective, to derive such constraints in earnest one should
implement the new energy sink in a state-of-the-art numerical simulation of SN
1987A for a plausible range of input assumptions (equation of state,
progenitor properties), derive the expected neutrino signal, and compare it
with the data. One could thus derive a quantitative confidence range for the
allowed particle coupling strength. In practice this has never been fully
done. Shortly after SN 1987A, besides many simple estimates [57], numerical
simulations with axion losses were performed [58, 59, 60], sometimes using the
time $\Delta t_{90\%}$ when 90% of the neutrino signal would have been
registered as a measure of signal duration [60].

One of us later developed a simple criterion to put the previous work on a
common footing: the new energy loss should not exceed $10^{19}\,{\rm
erg}\,{\rm g}^{-1}\,{\rm s}^{-1}$, or an overall luminosity of around $3\times
10^{52}\,{\rm erg}\,{\rm s}^{-1}$, to be calculated at nuclear density
$\rho=3\times 10^{14}\,{\rm g}\,{\rm cm}^{-3}$ and $T=30~{}{\rm MeV}$ [4, 20].
These conditions are meant to represent the PNS somewhat after the explosion
when the cold interior has heated up. The $T$ profile has a maximum (cf. Fig.
1) that moves inward as the core deleptonizes, so volume particle emission
will reach full strength only some time after core bounce. This behavior is
also manifest from Fig. 5, where we show a contour plot for the radial and
time evolution of the temperature, baryon density, and muon abundance for the
same SN model.

In contrast, the standard neutrino luminosity is largest at the beginning when
it is partly powered by accretion and may carry away up to half the full
amount during the first second. Particle emission removes energy from deep
inside the PNS that would otherwise power the late neutrino signal. So the SN
1987A signal duration is considered the main observable to constrain energy
losses from the PNS interior.

Figure 5: Time and radial evolution of the temperature (top panel), baryon
number density normalized to nuclear density $n_{0}=0.181\,{\rm fm}^{-3}$
(central panel), and muon abundance $Y_{\mu}$ for the SN model SFHo-18.8 [3]
that was also used in Fig. 1.

A modern analysis might be more constraining because contemporary models,
especially with muons, tend to be hotter. Moreover, PNS convection speeds up
cooling, leaving less room for a yet shorter signal.

Ideally, of course, the neutrino signal of the next galactic SN would become
available to obtain a high-statistics result. This would also overcome the
doubts that have been cast on the SN 1987A particle constraints because the SN
1987A neutron star has not yet clearly shown up [61], although ALMA radio
observations [62] are best explained as first evidence for a non-pulsar
compact remnant heating the overlying dust layer [63]. Improving particle
bounds would be one of the many benefits of observing the neutrino signal of
the next nearby SN.

####  III.1.2 Application to muonic bosons

As a first estimate we use this simple back-of-the-envelope recipe to
constrain bosons $b$ that couple to muons with Yukawa strength $g_{b}$. The
main emission process is photo production on muons that at first we treat as
nonrelativistic and nondegenerate. The scale for the production cross section
is set by the Thomson cross section of Eq. (4). So the energy-loss rate per
unit volume is simply the thermal energy density of photons,
$(\pi^{2}/15)\,T^{4}$, times the cross section times the muon number density
$n_{\mu}=Y_{\mu}n_{B}$. To obtain the energy loss per unit mass we divide by
the mass density777In SN physics, what is called the mass density $\rho$ is
usually defined as the number density of baryons times the atomic mass unit
$m_{u}=931.494\,{\rm MeV}$ that is based on the 12C atom. The local
gravitating energy density depends on the local composition and temperature
together with the equation of state. $\rho=n_{B}m_{u}$ so that

$\displaystyle\epsilon_{b}$ $\displaystyle=$
$\displaystyle\xi\,Y_{\mu}\,\sigma_{\rm
T}\,\frac{\pi^{2}}{15}\,\frac{T^{4}}{m_{u}}$ (21) $\displaystyle=$
$\displaystyle 1.70\times 10^{38}\,\frac{\rm erg}{\rm
g\,s}\,\xi\,g_{b}^{2}\,Y_{\mu}\left(\frac{T}{30\,{\rm MeV}}\right)^{4}$

where $\xi$ is a fudge factor accounting for corrections to the cross section.

To estimate the muon density we consider a PNS somewhat after SN explosion
when the outer core has deleptonized and heated. If we neglect the neutrino
chemical potentials, those of electrons and muons are the same. For nuclear
density and an assumed proton fraction of $Y_{p}=0.1$ and that the negative
$e$ and $\mu$ charges balance the protons implies $\mu_{e,\mu}\simeq
120~{}{\rm MeV}$ and $Y_{\mu}\simeq Y_{p}/3$, i.e. a muon fraction of
$Y_{\mu}\simeq 0.03$. The numerical model of Fig. 1 confirms this estimate.
The criterion $\epsilon_{b}<10^{19}\,{\rm erg}\,{\rm g}^{-1}\,{\rm s}^{-1}$
then implies

$g_{b}<1.4\times 10^{-9}\big{/}\sqrt{\xi}$ (22)

as a first limit.

However, for $T=30\,{\rm MeV}$ a typical photon energy of $3T$ is about the
same as the muon mass, so from Fig. 3 we glean that for scalars the semi-
Compton cross section is around 3 times smaller than the Thomson value, so
$\xi_{\phi}\simeq 1/3$. For vectors, the two final-state polarizations
introduce a factor of 2, so $\xi_{Z}\simeq 2/3$. For pseudoscalars finally
$\xi_{a}\simeq 1/10$, implying the bounds shown in Table 1.

Table 1: Muonic boson upper coupling limit from SN 1987A energy-loss argument. The simple estimate is from the often-used argument of Sec. III.1.2, whereas the Garching muonic SN models were used in Sec. III.1.3. We show two significant digits to avoid blurring the differences by rounding errors. Particle | Coupling | Simple | Numerical SN Models
---|---|---|---
| $[10^{-9}]$ | Estimate | Cold |  Hot
Scalar | $g_{\phi}$ | $2.4$ | $4.2$ |  $1.9$
Vector | $g_{Z}$ | $1.7$ | $2.7$ |  $1.2$
Pseudoscalar | $g_{a}$ | $4.4$ | $9.1$ |  $3.5$

####  III.1.3 Using the Garching SN models

To go beyond a back-of-the-envelope estimate we next calculate the emission
rates numerically for the Garching Group’s muonic SN models [3]. In their
Core-Collapse Supernova Archive [64] they provide radial hydrodynamical and
neutrino profiles for different equations of state and different progenitors
at many time shots between core bounce and 10 s afterwards. We use these
models to calculate the muonic-boson luminosity for each time shot and in Fig.
6 compare it with the instantaneous neutrino luminosity. Details about the
emission-rate calculation are given in Sec. III.1.4 below as well as redshift
corrections to the emission rate in Sec. III.1.5. What is shown in Fig. 6 are
luminosities in the reference frame of a distant observer.

Figure 6: Muonic boson luminosity of the numerical Garching models [3]
compared with their instantaneous $L_{\nu}$ as a function of postbounce time.
Luminosities are given for a distant observer. As in Fig. 1 we use the coldest
muonic model SFHo-18.8 (top) and in addition the hottest one LS220-20.0
(bottom). The shown boson couplings were taken such that in each case the
boson luminosity equals $L_{\nu}$ at 1 s.

Unsurprisingly, the boson luminosity is initially very small and only gets
larger as the core heats up and muons become abundant—the emission rate is
essentially proportional to $Y_{\mu}T^{4}$. On the other hand, $L_{\nu}$ is
largest at the beginning, mostly powered by accretion and energy from the
outer SN core layers. Therefore, as in the generic argument of Sec. III.1.2,
comparing $L_{\phi}$ with $L_{\nu}$ at around 1 s looks reasonable and we
adopt this criteria as a nominal bound. For the case of the coldest model that
was used in the upper panel of Fig. 6, the boson luminosity, calculated on the
basis of the unperturbed model, takes over and would carry away much of the
energy that otherwise would power late-time neutrino emission.

For pseudoscalars, the revised version of Ref. [3] found $g_{a}<8.4\times
10^{-9}$ based on the same numerical model, very similar to our bound,
although this agreement is somewhat fortuitous. There should have been a
factor of 2 in their emission rate (see our footnote 3). On the other hand,
they compared with $L_{\nu}=3\times 10^{52}~{}{\rm erg}~{}{\rm s}^{-1}$ from
the generic argument, whereas the native $L_{\nu}$ of the actual model,
ignoring redshift effects, is around $5.7\times 10^{52}~{}{\rm erg}~{}{\rm
s}^{-1}$, almost a factor of 2 larger. These factors approximately cancel and
the remaining differences may be due to different treatments of the cross
section.

Next we consider one of the hottest models of Ref. [3] and specifically use
LS220-20.0 which is based on a different equation of state. The inner
temperature is about a factor of 1.7 larger, so the $\phi$ emission rate about
a factor of 10 larger. The neutrino luminosity, on the other hand, is around a
factor of 2 larger. Notice that the neutron-star binding energy liberated in
the coldest model is $1.98\times 10^{53}\,{\rm erg}$ and $3.94\times
10^{53}\,{\rm erg}$ for the hottest one (see Table I of Ref. [3]). So the
nominal limit should be about a factor of 5 more restrictive on the luminosity
and a factor of 2–3 on the coupling constant.

However, the time profile of the boson luminosity, shown in the bottom panel
of Fig. 6, is quite different in that it drops much more quickly than in the
earlier case. Actually we have checked that for the hot model SFHo-20.0 with
the same equation of state as the cold model, the time profile looks very
similar to the latter, except that the boson luminosity overall is roughly 8
times larger. Therefore, the exact impact of boson emission on the neutrino
signal of SN 1987A probably would differ significantly between the two hot
models for the same boson coupling constant.

In Table 1 we compare the nominal limits from the coldest and hottest Garching
models thus obtained with the ones from our earlier back-of-the-envelope
generic argument. The estimates are somewhere between the two Garching
extremes and thus provide a reasonable magnitude. Without a detailed analysis
of the neutrino signal in the historical detectors for different cases, it is
not obvious if the data would clearly distinguish between these models or if
one of them rules new bosons out, the other might rule them in.

To be specific and conservative, we use the constraints derived from the
coldest numerical SN model in the summary plot of Fig. 2.

####  III.1.4 Emission rate

To complete this discussion we finally provide some details about our
numerical integration. For the emission rate we note that at $T\gtrsim 30$ MeV
the muons ($m_{\mu}=105.66$ MeV) are not strictly nonrelativistic, so one
cannot use the low-energy expansion of the cross section, and recoil effects
are important. No simple approximation is very good for these conditions, so
one should evaluate the boson production rate from the thermal environment,
including the appropriate photon and muon occupation numbers, by a numerical
evaluation.

However, this leads to a multi-dimensional numerical integral that we have
found too cumbersome to deal with in view of the large other uncertainties,
e.g. the appropriate SN 1987A model. Moreover, we have checked that the main
impact comes from the reduced cross section seen in Fig. 3. Therefore, we have
opted to calculate the emission rate in analogy to Eq. (21), determining the
fudge factor by averaging the cross sections given in Eq. (5) over a thermal
Bose-Einstein distribution of $\omega$.

####  III.1.5 Relativistic corrections

To integrate over the SN core we include redshift corrections and show the
neutrino and boson luminosities for a distant observer. The Garching SN models
include general-relativistic effects in an approximate way as described by
Rampp and Janka [65] and Case A of Marek et al. [66]. In practice this means
for us that the spatial coordinates should be interpreted as in flat space, so
the volume integral is performed with flat-space coordinates. Moreover, the
time coordinate is the one of a distant observer and does not require any
transformations.

However, the emitted particles suffer a gravitational redshift before reaching
infinity. In the tabulated models, this effect is encoded in the
“gravitational lapse” that is listed for every radius and is to be understood
as $(1+z)^{-1}$, where $z$ is the redshift. So the particle energy at infinity
is “local energy $\times$ lapse.” Moreover, the rate of emission suffers
another redshift factor, so the contribution to the luminosity at infinity by
a given spherical shell in the PNS is the naively calculated one times
(lapse)2 or times $(1+z)^{-2}$. Of course, this is only a 20–30% correction,
but it is trivial to include when post-processing the numerical models.

Moreover, the physical properties of the medium are given in Lagrange
coordinates and thus co-moving with the medium, also causing redshift effects.
Within the PNS the radial velocity $v_{r}$ is small, whereas at larger
distances it is not negligible and suffers a discontinuity at the shock-wave
radius, in turn causing a discontinuity in the tabulated $L_{\nu}(r)$. The
Doppler effect causes a redshift or blueshift of
$1+z=\sqrt{(1-v_{r})/(1+v_{r})}=1-v_{r}+{\cal O}(v_{r}^{2})$. Again $L_{\nu}$
is affected by the square of this factor, so in the appropriate limit of
$v_{r}\ll 1$ we multiply the tabulated $L_{\nu}(r)$ with $1+2v_{r}$ to
interpret $L_{\nu}$ in the rest frame of a distant observer.

After both corrections have been applied, $L_{\nu}(r)$ is constant beyond the
decoupling radius of around 18 km, but somewhat increases at very large
distances. The required time to reach a large distance means that $L_{\nu}(r)$
for a given time shot reflects earlier neutrino emission. To be specific and
to avoid this time-of-flight effect, we extract $L_{\nu}$ at 400 km, but the
exact radius is not important.

The reference model shown in Fig. 1 at 1 s pb has a neutrino luminosity at
infinity of $4.4\times 10^{52}\,{\rm erg}\,{\rm s}^{-1}$, whereas the maximum
value in local coordinates at around 16 km is 5.7 in these units, some 30%
larger than the value seen by a distant observer. In the upper panel of Fig. 6
we show $L_{\nu}(t)$ for this model, including the redshift corrections. The
integral up to the largest available time is $L_{\nu,{\rm tot}}=1.95\times
10^{53}\,{\rm erg}$, in good agreement with the final mass deficit of
$1.98\times 10^{53}\,{\rm erg}$ of this model [3]. Therefore, while redshift
corrections are not a huge effect, one needs to include them to obtain
consistent results and cross-checks.

###  III.2 Trapping limit

####  III.2.1 Introductory remarks

The main attraction of the SN 1987A energy-loss argument is that it probes
particles that are more feebly interacting than neutrinos, often providing
unique information. However, sometimes one may wish to consider new particles
that interact so strongly that they are trapped and can escape only from their
own decoupling region, in the axion case called the axiosphere in analogy to
the neutrino sphere [50, 51, 52, 53]. Typically such particles will be
excluded by other arguments, ranging from laboratory to cosmological evidence,
although we will see that for muonic bosons there is a sliver of parameter
space on the trapping side where SN arguments are unique.

Such particles can compete with neutrinos for several different tasks in SN
physics. They can radiatively transfer energy from deep inside the SN core to
the PNS surface in competition with neutrinos and convection and in this way
speed up PNS cooling. They can carry away some of the liberated neutron-star
binding energy, leaving less for detectable neutrinos. Even after beginning to
stream freely in their decoupling region, they can deposit some energy behind
the shock wave and contribute to reviving the shock wave in the delayed
explosion scenario. But also the opposite can be the case: a self-consistent
SN model may no longer have a gain radius within the shock radius beyond which
there is net energy deposition. The new particles could also show up in
neutrino detectors and could have caused some of the SN 1987A events [67].
Some of these effects may be more important than others, so it is hard to
formulate a generic argument.

A self-consistent simulation with axions in the trapping limit revealed that
the SN 1987A signal was shortened in terms of $\Delta t_{90\%}$, whereas the
number of events in the KamII and IMB detectors remained the same [53]. The
axion energy transfer heats the PNS surface, leading to larger emitted
neutrino energies, and thus to a larger detection rate.

Another self-consistent study, motivated by nuclear-physics issues, decreased
the effective neutrino cross section on nucleons, speeding up PNS cooling and
thus decreasing $\Delta t_{90\%}$, yet increasing the SN 1987A event numbers
by the same heating effect of the PNS surface [68].

Notice that the energy flux from the PNS surface is largely fixed deeply
inside, driven by the gradients of $T$ and lepton number. Because the surface
area is essentially fixed, the need to radiate a larger flux requires a larger
effective surface $T$. So the correlation between a reduced PNS cooling speed
and increased observable neutrino event energies is quite generic.

####  III.2.2 Simple bound from energy transfer

The energy flux carried by a trapped (pseudo)scalar boson at radius $r$ is
$L_{r}=-(\lambda/3)\,[d(aT_{r}^{4})/dr]\,4\pi r^{2}$ with $a=\pi^{2}/30$ is
the thermal energy content of one bosonic degree of freedom and $\lambda$ is
the MFP. As a simple estimate we use $dT/dr=-T/r$ with $T=30~{}{\rm MeV}$,
$r=10~{}{\rm km}$ and $L=3\times 10^{52}$ erg/s, all meant to mimic a
situation similar to the earlier energy-loss argument. (Notice that the model
in Fig. 1 shows a nearly constant $dT/dr\simeq-4~{}{\rm MeV}/{\rm km}$ for
$r=8$–15 km, which is stratified by convection and numerically very similar to
our estimate.) So the MFP needed to compete with standard energy transfer is
$\lambda=45L/8\pi^{3}rT^{4}=11$ m.

For our muonic scalar we use $\lambda^{-1}=(0.3\sigma_{\rm T})Y_{\mu}n_{0}$,
so $\lambda=g_{\phi}^{-2}~{}6.5\times 10^{-9}~{}{\rm cm}$, implying
$g_{\phi}\gtrsim 2\times 10^{-6}$ to avoid excessive energy transfer. We will
see, however, that avoiding excessive energy loss here provides a far more
restrictive limit because the density of muons quickly drops toward the PNS
surface. Therefore, the bosons may have rather large couplings to compete with
neutrinos in the decoupling region and then will be irrelevant for energy
transfer deeper inside.

####  III.2.3 Energy loss from the boson sphere

In the trapping limit, our bosons emerge from a region near the PNS surface
whence they escape without being reabsorbed on their way out, in analogy to
the neutrino sphere. If we picture them being emitted as blackbody radiation
by a spherical surface with radius $R_{\phi}$, the luminosity is given by the
Stefan-Boltzmann (SB) law

$L_{\phi}=4\pi R_{\phi}^{2}\,\frac{\pi^{2}}{120}\,T^{4}(R_{\phi}),$ (23)

where $T(R_{\phi})$ is the SN temperature at radius $R_{\phi}$.

To find $L_{\phi}$ as a function of coupling strength we determine $R_{\phi}$
by the requirement that the optical depth at that radius is
$\tau(R_{\phi})=2/3$ beyond which the bosons essentially stream freely. The
optical depth is defined as

$\tau(r)=\int_{r}^{\infty}dr^{\prime}\,\Gamma(r^{\prime}),$ (24)

where in natural units the interaction rate is the same as the inverse MFP,
$\Gamma=\lambda^{-1}$.

For our constraint we finally seek the coupling strength such that
$L_{\phi}=L_{\nu}$. In this context we define the Stefan-Boltzmann radius
$R_{\rm SB}$ of a given SN model by

$L_{\nu}=4\pi R_{\rm SB}^{2}\,\frac{\pi^{2}}{120}\,T^{4}(R_{\rm SB}).$ (25)

It is defined as the radius where the SB law for one boson degree of freedom
matches this model’s $L_{\nu}$. The SB radius is a property of a given SN
model without reference to any particular particle-physics assumption. Of
course, $R_{\rm SB}$ is always close to the neutrino-sphere radius $R_{\nu}$.

Our constraint then derives from the requirement that $\tau(R_{\rm SB})=2/3$,
i.e., the coupling strength is such that the boson-sphere radius $R_{\phi}$
coincides with the SB radius.

####  III.2.4 Reduced absorption rate

Before applying this method a number of clarifications are in order. To
compute the optical depth we need to distinguish carefully between different
variants of interaction rates that apply to our bosons where the MFP is
dominated by absorption (not scattering). The absorption rate, e.g. by inverse
bremsstrahlung, is termed $\Gamma_{\rm A}$, whereas the corresponding
spontaneous emission rate is $\Gamma_{\rm E}$. In the Boltzmann collision
equation for a given mode of the boson field with occupation number $f$ that
travels along some ray with spatial coordinate $s$, these quantities appear as

$(\partial_{t}+\partial_{s})\,f=\Gamma_{\rm E}(1+f)-\Gamma_{\rm
A}f=\Gamma_{\rm E}-\underbrace{(\Gamma_{\rm A}-\Gamma_{\rm
E})}_{\hbox{$\Gamma_{\rm A}^{*}\equiv\Gamma$}}f.$ (26)

Here $\Gamma_{\rm A}^{*}$ is called the reduced absorption rate, including
stimulated emission in the form of a negative absorption rate. If the medium
is in thermal equilibrium, detailed balance implies $\Gamma_{\rm
E}=e^{-\omega/T}\Gamma_{\rm A}$ so that the reduced absorption rate is

$\Gamma\equiv\Gamma_{\rm A}^{*}=\Gamma_{\rm A}(1-e^{-\omega/T}),$ (27)

which we use as the absorption rate. The spontaneous emission rate is then
expressed as

$\Gamma_{\rm E}=\frac{\Gamma}{e^{\omega/T}-1}.$ (28)

It is $\Gamma$, the reduced absorption rate, that appears in the optical-depth
integral of Eq. (24).

In a stationary and homogeneous situation, the lhs of Eq. (26) vanishes and
the equation is solved by a thermal Bose-Einstein distribution $f_{\rm
eq}=(e^{\omega/T}-1)^{-1}$. So we may write this equation instead for the
deviation from equilibrium $\Delta f=f-f_{\rm eq}$ and then reads

$(\partial_{t}+\partial_{s})\,\Delta f=-\Gamma\,\Delta f.$ (29)

So it is the reduced absorption rate $\Gamma$ which damps the deviation of $f$
from equilibrium, explaining its central importance for radiative transfer.

For our case of muonic bosons, the situation simplifies because we only
consider the absorption on either muons or the Primakoff conversion on charged
particles which are both of the type $\phi+X\to X+\gamma$. Therefore,
$\Gamma_{\rm A}=\Gamma_{\rm S}(1+f_{\gamma})$, where $\Gamma_{\rm S}=\sigma
n_{X}$ and $f_{\gamma}$ a boson stimulation factor for the final-state photon
in the thermal environment.888Here the cross section $\sigma$ is for the boson
in the initial state and includes a factor of 2 for the final-state photon
polarizations. The cross sections listed in Eq. (5), on the other hand, are
for photon scattering and thus averaged over the initial-state photon
polarizations. If the targets do not recoil much, the photon energy is nearly
the same as that of the boson, so $f_{\gamma}=1/(e^{\omega/T}-1)$ and
$1+f_{\gamma}=1/(1-e^{-\omega/T})$. So overall one finds

$\Gamma=\Gamma_{\rm A}^{*}=\Gamma_{\rm S},$ (30)

i.e., in our case the reduced absorption rate fortuitously is simply the naive
“cross section $\times$ target density.”999In some of the recent literature on
muonic bosons [69, 7] the optical depth was based on the un-reduced absorption
rate. On the other hand, in the updated version of Ref. [3] and in Ref. [49]
dedicated to ALPs the reduced opacity was used.

####  III.2.5 Spectral average

In our cases of interest the reduced opacity barely depends on boson energy
$\omega$, but in general this spectral dependence could be pronounced. In
differential form, the SB boson luminosity is

$\frac{dL_{\phi}}{d\omega}=4\pi
R_{\omega}^{2}\,\frac{1}{4}\,\frac{4\pi\omega^{3}}{(2\pi)^{3}}\,\frac{1}{e^{\omega/T(R_{\omega})}-1},$
(31)

where the factor $1/4$ includes one factor $1/2$ to count only the outgoing
modes of a blackbody distribution, another factor $1/2$ for the average
velocity of the outgoing modes because we are determining a flux. The final
factors include the boson phase space and Bose-Einstein occupation number.
Here $R_{\omega}$ is the radius where the $\omega$-dependent optical depth is
2/3. The critical coupling strength would be found by choosing the coupling
strength such that the spectral integral matches $L_{\nu}$.

Instead it may be more practical to use an average opacity and apply the SB
argument in integral form. One approach is to use the Rosseland mean opacity.
It is based on the diffusion limit where one boson degree of freedom at some
radius $r$ carries an energy flux given by

$\frac{1}{4\pi
r^{2}}\frac{dL_{\phi}}{d\omega}=-\frac{\lambda_{\omega}}{3}\,\nabla
B_{\omega},$ (32)

where for a massless boson the thermal energy density at energy $\omega$ is

$B_{\omega}=\frac{4\pi}{(2\pi)^{3}}\,\frac{\omega^{3}}{e^{\omega/T}-1},$ (33)

implying

$\nabla
B_{\omega}=\frac{1}{2\pi^{2}}\,\frac{e^{\omega/T}\omega^{4}/T}{(e^{\omega/T}-1)^{2}}\,\frac{\nabla
T}{T}.$ (34)

Here $(\nabla T)/T=\nabla\log T$ is the inverse length scale of temperature
decrease; the diffusion approximation is justified when this length scale is
large compared with $\lambda$.

For purely absorptive boson interactions as in our case, the MFP
$\lambda_{\omega}=\Gamma_{\omega}^{-1}$ is based on the reduced absorption
rate. If it does not depend on $\omega$, explicit integration yields

$\frac{1}{4\pi
r^{2}}\,L_{\phi}=-\frac{\lambda}{3}\,\frac{2\pi^{2}T^{4}}{15}\,\frac{\nabla
T}{T}=-\frac{\lambda}{3}\,\nabla(aT^{4})\,$ (35)

where $a=\pi^{2}/30$ is the radiation constant for a single boson degree of
freedom.

If the reduced absorption rate does depend on $\omega$, comparing Eq. (32)
with Eq. (35), the integrated flux can be written in the form

$\frac{1}{4\pi
r^{2}}\,L_{\phi}=-\frac{\langle\lambda\rangle}{3}\,\nabla(aT^{4})$ (36)

with the Rosseland average of the MFP

$\displaystyle\kern-20.00003pt\langle\lambda\rangle=\left\langle\Gamma^{-1}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{\int_{0}^{\infty}\\!d\omega\,\Gamma_{\omega}^{-1}\,dB_{\omega}/dT}{\int_{0}^{\infty}\\!d\omega\,dB_{\omega}/dT}$
(37) $\displaystyle=$
$\displaystyle\frac{15}{4\pi^{4}}\,\int_{0}^{\infty}\\!\\!dx\,\frac{1}{\Gamma({xT})}\,\frac{x^{4}\,e^{x}}{(e^{x}-1)^{2}},$

where we have used $\omega=xT$. The key point is that it is the MFP
$\lambda=\Gamma^{-1}$, not the interaction rate $\Gamma$, that is averaged
with a weight function derived from the black-body energy distribution.

While the Rosseland average is the appropriate quantity to compute the boson
radiative energy flux in the diffusion limit deep inside the PNS, it is less
obvious how good it quantifies the energy-dependent decoupling process in the
spirit of the integral SB approach.

####  III.2.6 Application to muonic bosons

We will see that in our context the spectral dependence of the decoupling
process is weak so that is most practical to apply the SB argument in integral
form, using the Rosseland average for the residual spectral opacity
dependence.

For scalars, we write the tree-level scattering rate on muons,
$\phi+\mu\to\mu+\gamma$, in the form

$\Gamma^{\phi,{\rm tree}}=\hat{Y}_{\mu}^{\phi,{\rm tree}}\,n_{B}\,2\sigma_{\rm
T},$ (38)

where $\sigma_{\rm T}$ is the Thomson cross section given in Eq. (4), the
factor of 2 counts two final-state photon polarizations, and $n_{B}$ is the
baryon density. For nonrelativistic and nondegenerate muons,
$\hat{Y}_{\mu}^{\phi,{\rm tree}}$ is the same as the muon abundance $Y_{\mu}$.
Otherwise it includes all corrections coming from muon degeneracy and the
deviation of the true cross section $\sigma_{\phi}$ given in Eq. (5) from the
Thomson value. Because $\sigma_{\phi}$ depends on $\omega$, the Rosseland
average needs to be taken. So finally

$\hat{Y}_{\mu}^{\phi,{\rm tree}}=\left\langle\frac{\sigma_{\phi}}{\sigma_{\rm
T}}\right\rangle Y_{\mu}$ (39)

is what we call the effective tree-level muon density for scalars.

The loop-induced two-photon coupling provides an additional source of opacity
that is irrelevant in the deep interior, but becomes dominant in the
decoupling region. The corresponding Primakoff scattering rate is

$\displaystyle\Gamma^{\phi,{\rm loop}}$ $\displaystyle=$ $\displaystyle\alpha
G_{\phi\gamma\gamma}^{2}\langle f_{\rm s}\rangle(1-Y_{n})n_{B}$ (40a)
$\displaystyle=$ $\displaystyle\frac{2\alpha^{2}}{3\pi^{2}}\,\langle f_{\rm
s}\rangle(1-Y_{n})n_{B}\,2\sigma_{\rm T}$ (40b)

where we have also included a factor of 2 for final-state photon polarizations
relative to Eq. (10), $\langle f_{\rm s}\rangle$ is the Rosseland average of
the screening factor of Eq. (11), and we have used Eq. (8a) for the scalar-
photon coupling. Comparing this expression with the tree-level contribution of
Eq. (38) we may define an effective muon density contributed by Primakoff
scattering of

$\hat{Y}_{\mu}^{\phi,{\rm loop}}=\frac{2\alpha^{2}}{3\pi^{2}}\,\langle f_{\rm
s}\rangle(1-Y_{n}).$ (41)

While the numerical factor $2\alpha^{2}/3\pi^{2}=3.60\times 10^{-6}$ is very
small, the loop contribution still dominates in the decoupling region. For our
reference SN model we show the radial profile of $\hat{Y}_{\mu}^{\phi,{\rm
tree}}$ and $\hat{Y}_{\mu}^{\phi,{\rm loop}}$ in Fig. 7; they cross at
$r=17.71~{}{\rm km}$.

Figure 7: Effective tree-level and loop-induced muon abundances in our cold
reference model for scalars and pseudoscalars as indicated and defined in the
text. The dashed red line is the physical muon abundance $Y_{\mu}$. The
Stefan-Boltzmann radius of this model is $R_{\rm SB}=16.97~{}{\rm km}$, where
the temperature is $T=7.43$ MeV.

The Stefan-Boltzmann radius of this model, i.e., the radius where the scalar
optical depth should be 2/3 so that the scalar luminosity matches $L_{\nu}$,
is $R_{\rm SB}=16.97~{}{\rm km}$, and the corresponding temperature $T=7.43$
MeV. We conclude that tree-level scattering somewhat dominates, but both
sources of opacity are important.

For pseudoscalars, we normalize the rates in the same way to $\sigma_{\rm T}$,
so the effective tree-level muon density is

$\hat{Y}_{\mu}^{a,{\rm tree}}=\left\langle\frac{\sigma_{a}}{\sigma_{\rm
T}}\right\rangle Y_{\mu}.$ (42)

We recall that for $\omega\ll m_{\mu}$ we have $\sigma_{a}/\sigma_{\rm
T}=\omega^{2}/m_{\mu}^{2}$, so here the effective muon density never
corresponds to the naive one. The relationship between the two-photon and
Yukawa coupling is now given by Eq. (8b), implying a larger loop-induced
contribution of

$\hat{Y}_{\mu}^{a,{\rm
loop}}=\left(\frac{3}{2}\right)^{2}\hat{Y}_{\mu}^{\phi,{\rm loop}}.$ (43)

So for pseudoscalars, $\hat{Y}_{\mu}^{a,{\rm tree}}$ is smaller,
$\hat{Y}_{\mu}^{a,{\rm loop}}$ larger relative to scalars and the radius of
equality is now deeper inside at $r=16.50~{}{\rm km}$ as seen in Fig. 7. As
this radius is smaller than $R_{\rm SB}$, the dominant opacity source is
Primakoff scattering so that the case of pseudoscalars is practically
identical to ALPs that interact only by their two-photon coupling.

With these ingredients we finally determine the Yukawa couplings such that the
optical depth at the SB radius is 2/3 and find

$\displaystyle g_{\phi}$ $\displaystyle>$ $\displaystyle 0.84\times 10^{-4}$
(44a) $\displaystyle g_{a}$ $\displaystyle>$ $\displaystyle 0.96\times
10^{-4}$ (44b)

as our nominal lower bounds from SN 1987A energy loss that we show in our
summary plot Fig. 2. If we use the hottest SN model LS220-20.0 instead, these
limits are $g_{\phi}>0.56\times 10^{-4}$ and $g_{a}>1.2\times 10^{-4}$.

We now use the limiting coupling strength thus identified at 1 s postbounce
and compute the SB luminosity at $\tau=2/3$ for all time shots and compare the
SB luminosity with $L_{\nu}$ in Fig. 8 for the cold and hot SN models. In
contrast to the free-streaming case, the boson luminosity here essentially
tracks $L_{\nu}$ because like $L_{\nu}$ it is generated in the PNS surface
region and thus governed by similar physical conditions. The idea that one
should derive a limit by a comparison at 1 s postbounce that was advanced as a
“modified luminosity constraint” is not particularly motivated, but on the
other hand, it makes little difference at which exact time one compares the
luminosities.

Figure 8: Muonic boson luminosity in the trapping limit of the numerical
Garching models [3] compared with their instantaneous neutrino luminosity as a
function of postbounce time. As in Fig. 6 we use the coldest muonic model
SFHo-18.8 (top panel) and in addition the hottest one LS220-20.0 (bottom). The
boson couplings were chosen such that the boson luminosity equals $L_{\nu}$ at
1 s.

One can perform the same exercise without the loop-induced Primakoff
contribution and instead use only the tree-level muon interaction. In this
case we find $g_{\phi}>1.1\times 10^{-4}$ and $g_{a}>6.2\times 10^{-4}$,
significantly more restrictive especially for pseudoscalars where the
scattering cross section on muons is reduced by the factor
$\omega^{2}/m_{\mu}^{2}$ as discussed earlier. However, these bounds do not
quite reach the $10^{-3}$ level and do not exclude the $g_{\mu}{-}2$ motivated
value. Moreover, one would be hard-pressed to take these nominal bounds very
seriously because they depend on the exponentially decreasing muon abundance
in the decoupling region. The muon distribution in this region would be
strongly affected by 3D effects in a non-spherically symmetric SN simulation.
Therefore, a 1D muonic model may be a poor proxy for this situation.

Including the loop effect, the limiting coupling strength is essentially set
by the Primakoff interaction channel, so the trapping case of our bosons is
practically identical with that of generic ALPs that interact only with
photons. So we now ignore the tree-level effect and again use the 1 s time
shot of the cold model. We thus find

$G_{\gamma\gamma}>2.1\times 10^{-6}~{}{\rm GeV}^{-1},$ (45)

where $G_{\gamma\gamma}$ stands for either $G_{\phi\gamma\gamma}$ or
$G_{a\gamma\gamma}$. For the hot model we find $G_{\gamma\gamma}>3.0\times
10^{-6}~{}{\rm GeV}^{-1}$.

Based on the SB argument, Lucente et al. [49] found the corresponding bound
$G_{\gamma\gamma}>7.7\times 10^{-6}~{}{\rm GeV}^{-1}$, a factor of 3.7 more
restrictive. This is a significant difference because it involves taking a
square root of the limiting flux. Besides a different SN model, sources of
difference include: (i) The Primakoff opacity that in Ref. [49] was based on
the proton abundance, ignoring small nuclear clusters. (ii) Integrating to
$\tau=2/3$ at a radius that was identified as the neutrino sphere on the basis
of assumed neutrino opacities.101010We thank the authors for explaining this
procedure in a private communication. It is somewhat different from what is
described in their Ref. [49].

It is plausible to expect that the SN 1987A neutrino signal will be strongly
affected when the boson luminosity is comparable to $L_{\nu}$, but of course
the exact modification has not been shown by a detailed analysis. The latter
would require implementing boson energy transfer and losses self-consistently
in a SN simulation, a formidable task probably not much less demanding than
neutrino transport itself.

####  III.2.7 Volume emission vs. SB approximation

In the recent literature on SN particle bounds, the SB procedure was critiqued
and instead a “modified luminosity constraint” was formulated for the trapping
limit [69, 6] and followed in subsequent papers [49]. It was correctly noted
that physically the boson emission was not blackbody emission from a
hypothetical axiosphere (or the equivalent for other bosons), but from an
extended volume in the outer regions of the star and it was asserted that the
luminosity thus determined exceeded the SB estimate.

For an unperturbed SN background model, the differential boson luminosity for
any degree of trapping is

$\frac{dL_{\phi}}{d\omega}=\int_{0}^{\infty}\\!\\!dr\,4\pi
r^{2}\,\frac{4\pi\omega^{3}}{(2\pi)^{3}}\,\Gamma_{\rm
E}(\omega,r)\,\bigl{\langle}e^{-\tau(\omega,r)}\bigr{\rangle},$ (46)

where the optical depth $\tau(\omega,r)$ is based on the reduced absorption
rate $\Gamma$ and then $\Gamma_{\rm E}$ is given by Eq. (28). The directional
average of the absorption factor is

$\bigl{\langle}e^{-\tau(\omega,r)}\bigr{\rangle}=\frac{1}{2}\int_{-1}^{+1}d\mu\,e^{-\int_{0}^{\infty}ds\,\Gamma\left(\omega,\sqrt{r^{2}+s^{2}+2rs\mu}\right)},$
(47)

where $\mu=\cos\beta$ and $\beta$ is the angle between the outward radial
direction and a given ray of propagation along which $ds$ is integrated.

So $\beta=0$ corresponds to the radial direction out of the star and thus to
the smallest optical depth $\tau_{0}(r)$ for a boson emitted in any direction
at radius $r$. The quantity $\tau_{0}(r)$ is called $\tau_{\rm radial,out}(r)$
in Eq. (B.1) of Ref. [69] or simply $\tau(\omega,r)$ in Eq. (7) of Ref. [3].
So if one uses $e^{-\tau_{0}}$ instead of $\langle e^{-\tau}\rangle$ as in
Ref. [3] one overestimates the luminosity in the trapping regime by a factor
of a few. In Eq. (B.4) of Ref. [69], and propagating into Eq. (4.3) of Ref.
[49], a geometrical correction to $e^{-\tau_{0}}$ was applied where typically
$\tau$ would actually become smaller than $\tau_{0}$, thus implying less
damping than the minimal possible amount. So the boson luminosity must have
been overestimated even more.

It is true, however, that $L_{\phi}$ for any degree of trapping is given by
the volume integral Eq. (46), so when postprocessing a numerical SN model, one
can evaluate this integral without reference to the SB argument and compare
$L_{\nu}$ with $L_{\phi}$ as a function of coupling strength.

We have investigated the comparison between the boson luminosity thus obtained
with the one provided by the SB approach and find them to agree very well if
the reduced absorption rate does not depend on $\omega$. One key element is
that in the decoupling region of the SN core the density falls much more
quickly as a function of radius than the temperature. If the absorption rate
is proportional to the density, as in the Primakoff case, and if we express
the temperature profile as a function of optical depth, then typically
$T\propto\tau^{p}$ with the power-law index $p\simeq 1/4$. (Notice that the
optical depth as a radial coordinate increases from outside in.)

The SB picture does not imply that the emitted bosons derive from a sharp
geometric radius. They are emitted from an extended region as shown, for
example, in Fig. 19 of Ref. [49], but this is not in contradiction with the SB
approach which does not imply that the emission region is a quasi-delta
function of geometric radius.

Making these arguments more precise is a somewhat extended exercise in the
physics of radiative transfer, especially when spherical geometry is included
as in Eq. (46). Actually it is somewhat magical how the volume integration of
Eq. (46) in the trapping regime turns itself into a quasi-surface emission
with all the right flux factors, although physically this has to be the case,
of course. When the reduced absorption rate depends strongly on $\omega$ it is
not clear if the Rosseland mean gives a very good approximate description, a
question that we have not yet investigated. We defer these discussion to a
dedicated future paper.

##  IV Explosion energy

Using the SN 1987A neutrino signal as a constraining argument is motivated by
thinking of the boson losses as an invisible energy-loss channel. However, the
loop-induced two-photon decay allows our muonic bosons to strongly show up in
electromagnetic radiation and, depending on parameters, can light up the
material of the SN progenitor, in the case of SN 1987A show up in the Gamma-
Ray Spectrometer on board the Solar Maximum Mission, or contribute to the
cosmic diffuse gamma-ray background from all past SNe.

We have seen that in the trapping limit, muonic boson decoupling is dominated
by Primakoff scattering because of the exponential decline of the muon density
near the PNS surface. Therefore, we now deal with generic scalar or
pseudoscalar ALPs in the $G_{\gamma\gamma}$–$m$ parameter space. The SN 1987A
energy-loss argument implies the lower limit $G_{\gamma\gamma}\gtrsim 2\times
10^{-6}~{}{\rm GeV}^{-1}$ given in Eq. (45).

From the two-photon decay rate of Eq. (17) follows, after including the
Lorentz factor, a MFP against radiative decay of

$\displaystyle\kern-20.00003pt\lambda_{a\to 2\gamma}$ $\displaystyle=$
$\displaystyle\frac{64\pi}{G_{\gamma\gamma}^{2}}\,\frac{E}{m^{4}}=4.0\times
10^{10}~{}{\rm cm}$ (48) $\displaystyle\times$
$\displaystyle\left(\frac{10^{-6}~{}{\rm
GeV}^{-1}}{G_{\gamma\gamma}}\right)^{2}\left(\frac{1\,{\rm
MeV}}{m}\right)^{4}\left(\frac{E}{10\,{\rm MeV}}\right).$

This is very much less than the Hubble scale, so the decay photons would show
up in the cosmic diffuse $\gamma$-ray background. We will see in Sec. VI that
at most about $10^{-4}$ of an average SN energy may appear in this form, so
one would conclude that $G_{\gamma\gamma}$ must be so large, the ALPs so
strongly trapped, that the overall energy they carry remains below this limit.

This argument, however, does not necessarily apply because the decay can be so
fast that it happens within the surrounding matter of the progenitor star. The
smallest relevant mass beyond the HB-star limit is around 0.2 MeV, extending
$\lambda_{a\to 2\gamma}$ to $2.5\times 10^{13}~{}{\rm cm}$. Typical core-
collapse progenitors are red supergiants and as such may have radii up to some
$1000\,R_{\odot}=6\times 10^{13}~{}{\rm cm}$, corresponding e.g. to the radius
of the red supergiant Betelgeuze. The progenitor of SN 1987A was the blue
supergiant Sanduleak $-69^{\circ}202$ with a somewhat smaller radius of $(3\pm
1)\times 10^{12}~{}{\rm cm}$ [70].

So for boson masses so small that they are covered by the HB-star argument,
the decays can happen beyond the progenitor-star radius and are thus also
constrained by diffuse cosmic $\gamma$ rays. Therefore, there is little
benefit in making the SN constraint more precise in this mass range. For
larger masses, where the HB-star argument no longer applies, the cosmic
diffuse $\gamma$ rays also cease to apply and we need to consider the energy
deposition in the progenitor star as the only relevant argument.

The outer layers of a red supergiant have a density of the order of
$10^{-8}~{}{\rm g}~{}{\rm cm}^{-3}$, corresponding to an electron density of
around $6\times 10^{15}~{}{\rm cm}^{-3}$. For 10 MeV $\gamma$ rays, the
Compton cross section is around $5\times 10^{-26}~{}{\rm cm}^{2}$, leading to
a MFP against Compton scattering of around $3\times 10^{9}~{}{\rm cm}$ which
is much smaller than the radius. Therefore, the entire energy emitted in ALPs
is dumped into the progenitor star and thus becomes visible in the form of SN
explosion energy and luminosity.

The idea of using the progenitor matter surrounding the collapsing core as a
calorimeter for decay photons was first advocated by Falk and Schramm [23] in
1978 in the context of putative radiative neutrino decays, years before the
now-standard delayed explosion paradigm was developed. The neutron-star
binding energy released in a core collapse is 2–$4\times 10^{53}~{}{\rm erg}$,
whereas a typical explosion energy is some $10^{51}~{}{\rm erg}$, so less than
1% of the total energy release shows up in the explosion (see also Ref. [71]
for a recent application of this criterion to the dark photon case).

The ALP energy deposition can be reduced to this limit only by a smaller flux
through a larger $G_{\gamma\gamma}$, causing the decays to be yet faster and
even more conservatively within the progenitor star. To reduce the boson
energy release to a level below 1% of $E_{\rm SN}$ we require

$G_{\gamma\gamma}>5.3\,(4.8)\times 10^{-5}~{}{\rm GeV}^{-1}$ (49)

based on the SB flux calculated on the unperturbed reference cold (hot) model.
The corresponding Yukawa couplings are

$\displaystyle g_{\phi}$ $\displaystyle>$ $\displaystyle 3.6\,(3.3)\times
10^{-3}$ (50a) $\displaystyle g_{a}$ $\displaystyle>$ $\displaystyle
2.4\,(2.2)\times 10^{-3}$ (50b)

shown in our summary plot Fig. 2. For the scalar case, this bound excludes the
$g_{\mu}{-}2$ explanation in that it is an order of magnitude more restrictive
than what is the required value of Eq. (1a).

These values were actually determined calculating the radius at which the SB
emission matches 1% of $L_{\nu}$ at the reference time of one second, and then
imposing the optical depth to be 2/3. The SB radius in the cold (hot) model is
42.1 (34.2) km, far beyond the neutrino sphere. This is also the reason why
the bounds from the cold and the hot models are basically the same, as the two
profiles are significantly different only in the inner cores. We have verified
that this procedure is not sensitive to the time at which the luminosity is
calculated. This matters because the Falk-Schramm argument refers to the time-
integrated energy deposition, not the instantaneous luminosity.

If we repeat this exercise for a less stringent requirement of a 10% energy
deposition, the SB radius for the cold (hot) model is at 21.4 (19.6) km,
closer to the neutrino sphere. The corresponding nominal bounds on the Yukawa
couplings are $g_{\phi}>0.81\,(0.75)\times 10^{-3}$ and
$g_{a}>0.54\,(0.51)\times 10^{-3}$, still excluding the scalar $g_{\mu}{-2}$
explanation.

However, our nominal “1% criterion” of a typical neutron star binding energy
is very conservative because SN explosion energies can be much smaller than
the canonical $10^{51}~{}{\rm erg}$. The class of subluminous type II plateau
SNe, besides having small Ni masses, also have small explosion energies even
below $10^{50}~{}{\rm erg}$ [24, 25, 26]. Reconstruction of the explosion
energy of SN 1054 that has led to the Crab Nebula and Pulsar reveals a value
around $10^{50}~{}{\rm erg}$ or less [27, 28]. Some or all of such low-energy
explosions could correspond to the lowest-mass progenitors that evolve as
electron-capture SNe [29]. So the most restrictive Falk-Schramm constraints
will arise from core collapses with the lowest-energy explosions.

On the other hand, our treatment is based on calculating the boson losses from
specific unperturbed numerical models. To take advantage of the lowest
observed explosion energies one should consider appropriate SN models. More
importantly, the bosons themselves will be the dominant agents of energy
transfer in the region between the neutrino and boson sphere, so one may not
necessarily assume one can compute the boson luminosity reliably from post
processing an unperturbed model. Conceivably one could develop an approximate
model of this region without solving the entire SN evolution self-
consistently, but this is a project for future research.

However, given the small amount of energy transfer to the progenitor-star
matter that is enough to violate the explosion-energy constraint, it looks
very hard to hide a scalar with the required coupling strength to explain the
muon magnetic-moment anomaly.

##  V Decay photons from SN 1987A

###  V.1 SMM observations

Photons from putative radiative decays of neutrinos or new particles emitted
by SN 1987A would have been picked up by the Gamma-Ray Spectrometer (GRS) on
board the Solar Maximum Mission (SMM) satellite that operated 02/1980–12/1989.
The GRS consisted of seven NaI detectors surrounded on the sides by a CsI
annulus and at the back by a CsI detector plate [72]. Because the GRS was
observing the Sun, $\gamma$-rays associated with the SN 1987A neutrino burst
would have hit almost exactly from the side and first had to traverse about
$2.5~{}{\rm g~{}cm^{2}}$ of spacecraft aluminum and $11.45~{}{\rm g~{}cm^{2}}$
of CsI shielding, effects that are small but were included to calculate the
effective detector areas [31, 32].

Table 2: GRS $3\sigma$ upper fluence limits for the two indicated time intervals after the arrival of the first SN 1987A neutrino. Channel | Energy band | Gamma fluence limits [cm-2]
---|---|---
| [MeV] | 10 s [31] | 223.2 s [32]
1 | 4.1–6.4 | 0.9 | 6.11
2 | 10–25 | 0.4 | 1.48
3 | 25–100 | 0.6 | 1.84

Neutrinos with up to 10 eV-range masses would not have strongly dispersed the
roughly 10 s SN 1987A burst and the same applies to the hypothetical burst of
decay photons. To constrain radiative decays, Chupp et al. [31] have analysed
the three energy channels shown in Table 2 for 10 s after the arrival of the
first IMB neutrino at UT 7:35:41.37 on 23 February 1987 where no excess counts
were found, leading to the fluence limits shown in Table 2. One key element of
this analysis was to construct and verify a time-dependent background model
because for a satellite in orbit the background rate is not constant. For the
spectrum of decay photons it was assumed that it is flat in the lowest channel
and that it follows $E_{\gamma}^{-2}$ for higher-energies.

If SN 1987A (distance 51.4 kpc) emitted $1\times 10^{53}~{}{\rm erg}$ in one
species of massive neutrinos ($\nu$ plus $\bar{\nu}$) with an average energy
of 15 MeV, the fluence at Earth was around $1.3\times 10^{10}~{}{\rm
cm}^{-2}$, so the SMM $\gamma$ limits imply that less than some $10^{-10}$ of
them should have decayed on their way to Earth.

While such limits may look impressive at first, what one is really
constraining is the underlying effective transition dipole moment $\mu_{\nu}$
and so the decay scales as $\mu_{\nu}^{2}m_{\nu}^{3}$ while the relativistic
Lorentz factor provides another factor $m_{\nu}/E_{\nu}$, so we are punished
with an $m_{\nu}^{4}$ phase-space factor. A similar remark applies to ALPs
where the decay rate scales as $G_{a\gamma\gamma}^{2}m_{a}^{3}$ and in both
cases theoretical models for the effective photon coupling often introduce yet
more powers of the mass. Therefore, typically much more information is gained
from looking at processes such as the plasmon decay $\gamma_{\rm
pl}\to\nu\bar{\nu}$ in stars [73] or the coherent Primakoff conversion of very
low-mass ALPs in astrophysical $B$-fields (see Refs. [48, 74] for recent
discussions and references to the earlier literature).

Motivated by the option of MeV-mass $\tau$ neutrinos, still a distinct
possibility 30 years ago, Oberauer et al. [32] in 1993 extended this
discussion to higher masses. In this case time-of-flight dispersion extends
the hypothetical $\gamma$ signal to a much longer time than the original burst
of low-mass neutrinos. Therefore, these authors considered the signal up to
232.2 s, after which the GRS went into a calibration mode for 10 min. More
data are available later until the satellite passed through the South Atlantic
radiation anomaly and the instruments were switched off. The later data was
not used because it would have required a dedicated background study. The
fluence limits for the 232.2 s interval are shown in Table 2.

###  V.2 Gamma signal from massive-particle decays

To predict the $\gamma$ burst from decaying particles with masses above a few
tens of eV, several new effects need to be included. The first decays happen
directly at the source, so the first photons arrive contemporaneously with the
first (massless) neutrinos, but the signal is stretched because, strictly
speaking, the decays never stop, even when the massive neutrinos have passed
the Earth. Of course, the first photons really come from decays outside of the
progenitor star, not immediately from the PNS surface, so in this sense there
is a brief delay for the onset of the signal. Moreover, the laboratory-frame
photon energies are the boosted rest-frame energies and thus depend on the
rest-frame emission angle. The received photons come from a small angle away
from the line of sight, depending on emission angle, so that the trajectory of
the parent neutrino and decay photon trace out a triangle, implying a larger
time of flight for a larger emission angle. So for a fixed rest-frame energy
of emission, the length of the photon burst is larger for smaller received
energies. Taking all of these effects into account, Oberauer et al. [32]
derived a general expression for the time-energy structure of the expected
signal. Then several simplifications were made: (i) The neutrino burst of a
few seconds is treated as instantaneous emission and so we only need the
fluence spectrum of parent particles $\Phi_{a}(E_{a})$ at Earth (units ${\rm
cm}^{-2}~{}{\rm MeV}^{-1}$), i.e., the time-integrated flux passing the Earth
if there were no decays. (ii) The short period between leaving the PNS and
passing the progenitor surface is ignored, so the signal onset is a step
function at the time of the first massless neutrino. (iii) Photons arriving
within the 232.2 s interval come from a spatial region around the source which
is much smaller than our distance to SN 1987A. (iv) We use isotropic rest-
frame emission (in Ref. [32] the neutrino Majorana case), but also appropriate
for (pseudo)scalar bosons. (v) We include a factor of 2 for two photons per
decay. (vi) The parent particle is taken to be very relativistic so that the
range of lab-frame photon energies is a box spectrum on the interval 0–$E_{a}$
if $E_{a}$ is the lab-frame parent energy. In this case the expected
$\gamma$-ray flux spectrum according to Eq. (18) of [32] is, later also
confirmed e.g. in Eq. (5) of Ref. [75],

$\frac{dF_{\gamma}}{dE_{\gamma}dt}=2\,\frac{2E_{\gamma}}{m_{a}\tau_{a}}\,e^{-2E_{\gamma}t/m_{a}\tau_{a}}\int_{E_{\gamma}}^{\infty}dE_{a}\,\frac{\Phi_{a}(E_{a})}{E_{a}},$
(51)

where $\tau_{a}$ is the rest-frame boson lifetime. In practice it will be so
long that we can neglect the exponential. In this case, the time structure is
flat, so the hypothetical signal should show up as a sudden upward jump of the
GRS counting rate at the time of the neutrino burst.

To check the key assumptions explicitly, we recall that the time-of-flight
difference to travel a distance $L$ between a massless particle and one with
mass $m_{a}$ is

$\Delta t=\frac{m_{a}^{2}}{2E_{a}^{2}}\,L.$ (52)

With $L=3\times 10^{12}~{}{\rm cm}$ for the radius of the SN 1987A progenitor
and with our most extreme mass $m_{a}=10~{}{\rm MeV}$ and $E_{a}=100~{}{\rm
MeV}$ as a typical boson energy, one finds $\Delta t=0.5~{}{\rm s}$. So the
envelope absorption at the source cuts off only a negligible period at the
signal onset.

Conversely we may ask where those photons originate that we see at the end of
the observation period of 223.2 s. Taking now our smallest mass of interest,
$m_{a}=0.1~{}{\rm MeV}$, with $\Delta t=223.2~{}{\rm s}$ we find $L=14~{}{\rm
yr}$ to be compared with the distance to SN 1987A of 160,000 yr. So indeed the
detectable decays would happen within a very small angular range around the
direction of the source.

###  V.3 ALP decays

Recently the SMM data were used to constrain radiative ALP decays, i.e.,
pseudoscalars that are emitted from the SN core by their photon coupling alone
[33]. These authors considered all of the effects leading to Eq. (51), but
being unaware of Ref. [32] they derived the detection spectrum by a Monte
Carlo simulation instead of an analytical integration. Therefore, we here
briefly reconsider this case as a mutual test of consistency.

From their Fig. 8, the lower edge of the excluded region corresponds to

$G_{a\gamma\gamma}<1.75\times 10^{-11}~{}{\rm GeV}^{-1}\,\sqrt{\frac{\rm
MeV}{m_{a}}}\,,$ (53)

where the scaling with mass is stated in their Eq. (19) but otherwise follows
from their numerical analysis.

The fluence spectrum $\Phi_{a}(E_{a})$ was obtained from the time integration
of a $18.0\,M_{\odot}$ SN model of the Wroclaw group. The same model was
previously used to constrain $G_{a\gamma\gamma}$ from the galactic $B$-field
conversion of very low-mass ALPs [48], where many details of the SN model are
documented. The analytic representation for $\Phi_{a}(E_{a})$ provided in Ref.
[33] is an excellent match to the numerical result as shown in their Fig. 3
and corresponds to a total number of emitted ALPs of $N_{a}=5.26\times
10^{53}\,G_{10}^{2}$ where $G_{10}=G_{a\gamma\gamma}/10^{-10}~{}{\rm
GeV}^{-1}$. The average axion energy is $\left\langle
E_{a}\right\rangle=102.0\,{\rm MeV}$, in agreement with typical interior
temperatures of around 30–35 MeV.

Inserting their analytic $\Phi_{a}(E_{a})$ into Eq. (51), ignoring the
exponential, and demanding that multiplied with 223.2 s the $\gamma$-fluence
in Channel 3 obeys the $3\sigma$ limit given in Table 2 of $1.84~{}{\rm
cm}^{-2}$, we find

$G_{a\gamma\gamma}<1.53\times 10^{-11}~{}{\rm GeV}^{-1}\,\sqrt{\frac{\rm
MeV}{m_{a}}}\,,$ (54)

similar to Eq. (53). However, the limit on $G_{a\gamma\gamma}$ involves taking
the fourth root. So in terms of the predicted photon fluence, our bound is a
factor 1.7 more restrictive, i.e., we have a larger photon fluence by this
factor. Actually they used a slightly more restrictive fluence limit of
$1.78~{}{\rm cm}^{-2}$ instead of our $1.84~{}{\rm cm}^{-2}$, making the
discrepancy slightly worse. We have carefully checked the derivation of Eq.
(51) but could not find any extraneous factor of 2 that might have crept in.
Alternatively, a small error may have sneaked into the Monte Carlo
implementation of Ref. [33].111111We thank J. Jaeckel for a private
communication, explaining that the probable origin is insufficient resolution
of the original MC.

For comparison we may perform the same analysis for our cold reference model
that has similar interior temperatures of around 30 MeV. We find
$N_{a}=1.73\times 10^{53}\,G_{10}^{2}$ with $\left\langle
E_{a}\right\rangle=89.9\,{\rm MeV}$, providing

$G_{a\gamma\gamma}<2.0\times 10^{-11}~{}{\rm GeV}^{-1}\,\sqrt{\frac{\rm
MeV}{m_{a}}}\,.$ (55)

For a similar interior $T$ our model emits a factor of 3 fewer axions, a
difference that follows from the time period of emission. The average time of
ALP emission in our model is $\left\langle t_{\rm e}\right\rangle\simeq
2\,{\rm s}$, whereas in their model it is $6\,{\rm s}$. The long time it takes
for their model to cool can also be seen, e.g., in Fig. 1 of Ref. [48].
Conversely, the short time scale of our model is probably explained by the
role of PNS convection in the muonic Garching models.

Finally the same exercise for our hot model provides $N_{a}=1.28\times
10^{54}\,G_{10}^{2}$ with $\left\langle E_{a}\right\rangle=136.8\,{\rm MeV}$
and almost the same average time of emission of $\left\langle t_{\rm
e}\right\rangle\simeq 2\,{\rm s}$. The corresponding constraint is

$G_{a\gamma\gamma}<1.24\times 10^{-11}~{}{\rm GeV}^{-1}\,\sqrt{\frac{\rm
MeV}{m_{a}}}.$ (56)

This result is most restrictive because of the large interior $T$ of the hot
model.

The constraint on $G_{a\gamma\gamma}$ involves taking a fourth root because
the coupling strength enters both at production and decay. Even though the
total number of emitted axions varies by a large factor between the models,
the spread of the limiting $G_{a\gamma\gamma}$ is only a factor of 1.6.

Of course, it may not be completely arbitrary which model best represents SN
1987A. In principle one could derive the expected neutrino signal and compare
it with the historical data. Conceivably one could discriminate between the
models. Notice that here ALP emission is but a small perturbation because the
constraint comes from ALP decays, not from the backreaction on the PNS cooling
speed, so this comparison would be between the unperturbed numerical models.

Table 3: Limits on the Yukawa couplings of (pseudo)scalar muonic bosons from the Garching muonic SN models and the SMM $\gamma$-ray limits. (Notice that here we do not consider the vector case because vectors do not decay into photons.) Model | Boson | $N_{a,\phi}$ | $\langle E_{a,\phi}\rangle$ | $t_{\rm e}$ | Limit on $g_{a,\phi}$
---|---|---|---|---|---
|  | $\times g_{a,\phi}^{2}$ | [MeV] | [s] | $\times\sqrt{{\rm MeV}/m_{a,\phi}}$
Cold | $a$ | $1.48\times 10^{73}$ | 105.8 | 2.3 | $1.74\times 10^{-10}$
| $\phi$ | $11.5\times 10^{73}$ | 67.8 | 2.5 | $1.11\times 10^{-10}$
Hot | $a$ | $13.4\times 10^{73}$ | 142.3 | 2.4 | $1.01\times 10^{-10}$
| $\phi$ | $78.8\times 10^{73}$ | 92.5 | 2.9 | $0.68\times 10^{-10}$

###  V.4 Muonic bosons

We finally turn to our main case of interest and calculate the fluence
$\Phi_{a}$ and $\Phi_{\phi}$ for (pseudo)scalars based on the photo production
process $\gamma+\mu\to\mu+a$ or $\phi$ as a function of the Yukawa couplings,
both for our cold and hot reference models. For future reference we provide a
simple fit function for the boson fluence in terms of a Gamma distribution

$\frac{dN_{a,\phi}}{dE_{a,\phi}}=C\left(\frac{E_{a,\phi}}{\langle
E_{a,\phi}\rangle}\right)^{\alpha}e^{-(\alpha+1)\frac{E_{a,\phi}}{\langle
E_{a,\phi}\rangle}}$ (57)

where $\langle E_{a,\phi}\rangle$ is the average energy of the boson fluence
as reported in Tab. 3, $\alpha$ represents the amount of spectral pinching,
and $C$ is an overall normalization. For $\alpha=2$ we recover a Maxwell-
Boltzmann distribution with $\langle E_{a,\phi}\rangle=3T$. For the cold model
we find these parameters to be $\alpha=0.937$ (2.49) and $C=6.42\,(3.38)\times
10^{72}\,\rm MeV^{-1}$ for the scalar (pseudoscalar) case. For the hot model
we find $\alpha=0.77\,(2.08)$ and $C=2.54\,(1.41)\times 10^{73}\,\rm
MeV^{-1}$.

Using the numerical fluence results to predict the expected $\gamma$ fluence
we finally find the constraints shown in Table 3. The cold model provides the
most conservative constraints which we use in our summary plot of Fig. 2. For
masses larger than $0.70~{}(0.47)$ keV for the scalar (pseudoscalar) case
these results are more restrictive than those from the SN 1987A energy-loss
argument.

##  VI Diffuse gamma-ray background

###  VI.1 Redshift integral

The (pseudo)scalar boson emission from all past SNe creates a cosmic
background density in analogy to the diffuse SN neutrino background (DSNB)
[76, 77, 78, 56, 79]. Radiative decays of these particles contribute to the
diffuse cosmic $\gamma$-ray background and thus can be constrained, an idea
that probably goes back to an early paper by Cowsik [80]. Because of the
phase-space issue discussed earlier, such arguments are not especially
powerful for low-mass particles such as ordinary neutrinos, but very useful
for MeV-range bosons that we consider. Recent discussions include ALPs that
are emitted by their photon coupling [34] or by processes involving nucleons
or electrons [81].

Closely following Ref. [82] we recall that the SN boson density spectrum
accumulated from all cosmic epochs today is given by the redshift integral

$\frac{dn_{a}}{dE}=\int_{0}^{\infty}\\!\\!dz\,(z+1)\,F_{a}(E_{z})\,n^{\prime}_{\rm
cc}(z),$ (58)

where $F_{a}(E)=dN_{a}/dE$ is the boson number spectrum emitted by an average
SN; the integral is the total number $N_{a}$ of emitted bosons. Moreover,
$E_{z}=(1+z)E$ is the blue-shifted energy at emission of the arrival energy
$E$. The first factor $(1+z)$ is a Jacobean $dE_{z}/dE=(1+z)$ between emitted
and detected energy interval.

Finally $n^{\prime}_{\rm cc}(z)=dn_{\rm cc}/dz$ is the core-collapse number
per comoving volume per redshift interval. It is usually expressed as

$n^{\prime}_{\rm cc}(z)=R_{\rm cc}(z)\,t^{\prime}(z)$ (59)

where the rate-of-change of cosmic time with regard to redshift is

$t^{\prime}(z)=\frac{dt}{dz}=\frac{1}{H_{0}\,(1+z)\sqrt{\Omega_{\rm
M}(1+z)^{3}+\Omega_{\Lambda}}}.$ (60)

Here $H_{0}$ is the Hubble expansion parameter, while $\Omega_{\rm M}$ and
$\Omega_{\Lambda}$ are the present-day cosmic matter and dark-energy
fractions. In the usual flat $\Lambda$CDM cosmology
$\Omega_{\Lambda}=1-\Omega_{\rm M}$. In the literature one usually finds
$R_{\rm cc}(z)$, the number of core collapses per comoving volume per unit
time (units ${\rm Mpc}^{-3}~{}{\rm yr}^{-1}$). However, $R_{\rm cc}(z)$ is
derived in terms of an assumed cosmological model because observations for a
given redshift interval need to be translated to intervals of cosmic time, so
only $n^{\prime}_{\rm cc}(z)$ has direct meaning. We will use $H_{0}=70~{}{\rm
km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}=(13.9~{}{\rm Gyr})^{-1}$, $\Omega_{\rm
M}=0.3$, and $\Omega_{\Lambda}=0.7$. These are not the latest best-fit
parameters, but consistent with $R_{\rm cc}(z)$ that we will use.

To derive the present-day number density of decay photons we assume that the
boson mass is sufficiently small to treat them as ultra-relativistic. In this
case the spectrum of decay photons is box-shaped on the range 0–$E_{a}$ where
$E_{a}$ is the particle energy at the instant of decay. Moreover, in the
ultra-relativistic limit the parent energy and that of a decay photon redshift
the same way, so this condition does not depend on when the decay takes place.
However, the rate of decay at redshift ${z_{\rm D}}$ does depend on the cosmic
epoch through the Lorentz factor

$\Gamma_{{z_{\rm D}}}=\frac{1}{\tau_{a}}\,\frac{m_{a}}{E_{z_{\rm D}}},$ (61)

where $E_{z_{\rm D}}$ is the parent energy at redshift ${z_{\rm D}}$.

Integrating over the photon spectrum provided by the bosons produced at
redshift $z$, the present-day spectrum of decay photons is found to be

$\frac{dn_{\gamma}}{d\omega}=\int_{0}^{\infty}\\!\\!dz(1+z)n^{\prime}_{\rm
cc}(z)\int_{\omega_{z}}^{\infty}\\!\\!dE_{z}f_{\rm
D}(E_{z})\,\frac{2}{E_{z}}\,F_{a}(E_{z}),$ (62)

where $\omega_{z}=(1+z)\omega$ and the factor of 2 represents two decay
photons on the interval 0–$E_{z}$. The first factor $1+z$ now represents the
Jacobean between the detected and emitted photon energy—the spectrum of decay
photons is evaluated at the fixed redshift $z$ of emission because in our
ultra-relativistic approximation the relation between energy of parent boson
and decay photons is independent of redshift. Our result, Eq. (62), agrees
with the corresponding expression derived in Ref. [83] for the nonradiative
two-body decays of the SN relic neutrinos.

Finally we need the fraction of bosons that has decayed between the epoch of
emission at redshift $z$ and today

$f_{\rm D}(E_{z})=1-\exp\left[-\int_{0}^{z}d{z_{\rm
D}}\,\frac{t^{\prime}({z_{\rm D}})}{\tau_{a}}\,\frac{m_{a}}{E_{z_{\rm
D}}}\right],$ (63)

where $E_{z_{\rm D}}=E_{z}(1+{z_{\rm D}})/(1+z)$ is the energy of a boson at
the decay redshift ${z_{\rm D}}$ if it had energy $E_{z}$ at emission.
Therefore, the integral can be written as

$\int_{0}^{z}d{z_{\rm D}}\,\frac{t^{\prime}({z_{\rm
D}})}{\tau_{a}}\,\frac{m_{a}}{E_{z_{\rm
D}}}=\frac{m_{a}}{E_{z}}\,\frac{1}{H_{0}\tau_{a}}\,g_{\rm D}(z).$ (64)

Here the function of redshift

$g_{\rm D}(z)=(1+z)\int_{0}^{z}\frac{dz_{\rm D}}{(1+z_{\rm
D})^{2}\sqrt{\Omega_{\rm M}(1+z_{\rm D})^{3}+\Omega_{\Lambda}}}$ (65)

depends only on the assumed cosmological model.

For bosons that decay within a Hubble time, $f_{\rm D}=1$. In the opposite
limit of long-lived bosons the exponential can be expanded, so the spectrum of
decay photons is

$\frac{dn_{\gamma}}{d\omega}=\frac{m_{a}}{H_{0}\tau_{a}}\int_{0}^{\infty}\\!\\!dz(1+z)n^{\prime}_{\rm
cc}(z)g_{\rm
D}(z)\int_{\omega_{z}}^{\infty}\\!\\!dE_{z}\frac{2F_{a}(E_{z})}{E_{z}^{2}}.$
(66)

So we needed to replace $f_{\rm D}(E_{z})\to(m_{a}/H_{0}\tau_{a})\,g_{\rm
D}(z)/E_{z}$.

###  VI.2 Cosmic core-collapse rate

One central ingredient is the cosmic core-collapse rate $n^{\prime}_{\rm
cc}(z)$. The starting point is the comoving star-formation rate
$\dot{\rho}_{*}(z)$ for which several groups provide results in the form of
analytic approximation functions [84, 85, 86, 87]. Widely used is the one of
Yüksel et al. [84] which is piecewise linear in $\log z$ of the form

$\dot{\rho}_{*}(z)=\dot{\rho}_{0}\left[(1+z)^{a\eta}+\left(\frac{1+z}{B}\right)^{b\eta}+\left(\frac{1+z}{C}\right)^{c\eta}\right]^{1/\eta},$
(67)

where $\dot{\rho}_{0}=0.02\,M_{\odot}\,{\rm Mpc}^{-3}~{}{\rm yr}^{-1}$,
$a=3.4$, $b=-0.3$, $c=-3.5$, $B=5000$, $C=9$, and the somewhat arbitrary
smoothing parameter is $\eta=-10$.

This star-formation rate is converted to a cosmic core-collapse rate $R_{\rm
cc}(z)=k_{\rm cc}\dot{\rho}_{*}(z)$ with the coefficient $k_{\rm
cc}=(135M_{\odot})^{-1}$ as explained in Ref. [82] and very similar to
$(143M_{\odot})^{-1}$ of Ref. [77]. We finally multiply with $t^{\prime}(z)$
of Eq. (60) to obtain $n^{\prime}_{\rm cc}(z)$ shown in Fig. 9 (top line).

Figure 9: Cosmic core-collapse rate $n^{\prime}_{\rm cc}(z)$ from top to
bottom according to (1) Yüksel et al. [84], (2) Mathews et al. [85], (3)
Robertson et al. [86], and (4) Madau and Dickinson [87].

A similar representation was obtained by Mathews et al. [85] with somewhat
different fit parameters, also shown in Fig. 9 (second line from top). Their
stated uncertainty is nearly a factor of 2, depending on redshift.

A somewhat different result and representation was found by Robertson et al.
[86] who also used slightly different cosmological parameters. Transformed to
our reference cosmology their result is the third line from top in Fig. 9.
They also state a large uncertainty range, comparable to that of Mathews et
al. [85].

We finally consider explicitly the star-formation rate of Madau and Dickinson
[87] who provided

$\dot{\rho}_{*}(z)=\frac{0.015\,M_{\odot}}{{\rm Mpc}^{3}~{}{\rm
yr}}\,\frac{(1+z)^{2.7}}{1+[(1+z)/2.9]^{5.6}}.$ (68)

It leads to the lowest core-collapse rate in Fig. 9.

The total number of past core collapses is $n_{\rm cc}=1.05$, 0.83, 0.68, and
$0.58\times 10^{7}~{}{\rm Mpc}^{-3}$ in the four cases from top to bottom. In
other words, the spread between these different cases as well as the internal
uncertainties stated by some of them roughly span a factor of 2. Indeed, one
of the goals of measuring the DSNB in the gadolinium-enhanced Super-Kamiokande
and in the forthcoming JUNO scintillator detectors is to settle the overall
normalization of the cosmic star-formation rate. To represent the uncertainty
from $n_{\rm cc}$ we will often express our results in terms of the parameter

${n^{\rm cc}_{7}}=\frac{n_{\rm cc}}{10^{7}~{}{\rm Mpc}^{-3}}.$ (69)

It varies between 0.58 and 1.05 between the lowest and highest of the cited
star-formation rates, but of course a broader range can be considered.

###  VI.3 Generic limit

####  VI.3.1 Short-lived bosons

As a first simple case we consider bosons that decay so fast relative to the
Hubble time that we can set the decay fraction $f_{\rm D}=1$. Moreover, to
derive a first estimate of the expected $\gamma$-ray signal we assume that a
fraction $\zeta_{a}$ of a typical SN energy release of $E_{\rm SN}=3\times
10^{53}~{}{\rm erg}$ is emitted in ALPs with average energy $E_{\rm av}$. The
spectrum will be quasi-thermal, so as a rough description we model it as a
Maxwell-Boltzmann distribution of the form

$F_{a}(E)=\zeta_{a}\,\frac{E_{\rm SN}}{E_{\rm
av}}\,\frac{E^{2}}{2T^{3}}\,e^{-E/T},$ (70)

where $E_{\rm av}=3T$. We actually anticipate that the detailed energy
distributions of the emitted boson has a very small impact on the results.

The SN redshift distribution is concentrated at $z=1$, so for now we assume
that all SNe occur at exactly this redshift, i.e., $n_{\rm cc}^{\prime}=n_{\rm
cc}\delta(z-1)$, leading to

$\frac{dn_{\gamma}}{d\omega}=\zeta_{a}\,\frac{E_{\rm SN}}{E_{\rm av}}n_{\rm
cc}\,\frac{2(T+2\omega)}{T^{2}}\,e^{-2\omega/T_{\alpha}}$ (71)

as our prediction for the present-day cosmic $\gamma$ density.

To compare this diffuse $\gamma$ flux with the measurements summarized in Fig.
10 we observe that in our region of interest of 2–200 MeV the spectrum is
essentially flat when multiplied with $\omega^{2}$, i.e., the measured flux in
this range is approximately

$\omega^{2}\frac{d\Phi_{\gamma}}{d\omega}\simeq 2\times 10^{-3}~{}{\rm
MeV}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}~{}{\rm ster}^{-1}.$ (72)

Notice that our prediction is on the $\gamma$ density, whereas the flux per
solid angle shown in Fig. 10 requires a factor of the speed of light (which is
$c=1$ in natural units) and to be divided by $4\pi\,\rm ster$, so our
prediction is translated as
$d\Phi_{\gamma}/d\omega=(dn_{\gamma}/d\omega)/4\pi$ and then is a flux per
ster. So in this form our prediction is

$\omega^{2}\frac{d\Phi_{\gamma}}{d\omega}=\zeta_{a}\,\frac{E_{\rm SN}}{E_{\rm
av}}n_{\rm cc}\,\frac{2(T+2\omega)\omega^{2}}{4\pi T^{2}}\,e^{-2\omega/T}.$
(73)

This is again a quasi-thermal shape with a maximum at $\omega_{\rm
max}=T(1+\sqrt{3})/2\simeq 1.37\,T$. Inserting this in Eq. (73), using $E_{\rm
av}=3T$, and the fudge factor ${n^{\rm cc}_{7}}$ defined in Eq. (69) we find

$\displaystyle\kern-20.00003pt\omega^{2}\frac{d\Phi_{\gamma}}{d\omega}\Big{|}_{\rm
max}$ $\displaystyle=$ $\displaystyle\zeta_{a}\,E_{\rm SN}n_{\rm
cc}\,\frac{7+4\sqrt{3}}{12\pi}\,e^{-(1+\sqrt{3})}$ (74) $\displaystyle=$
$\displaystyle\zeta_{a}\,{n^{\rm cc}_{7}}\,46.2~{}{\rm MeV}~{}{\rm
cm}^{-2}~{}{\rm s}^{-1}~{}{\rm ster}^{-1}.$

The nice key point of this expression is that it does not depend on the
assumed $T$ and that the observed spectrum is flat in the region of interest.

Figure 10: The extragalactic background light (EBL) over a large range of
energies according to Ackermann et al. [88]. (Figure reproduced with
permission.)

Comparing this prediction with Eq. (72) we finally find our generic limit

$\zeta_{a}\lesssim 0.43\times 10^{-4}\big{/}{n^{\rm cc}_{7}}$ (75)

independently of the assumed average energy of the emitted bosons. This result
applies both to bosons emitted near the PNS surface in the trapping limit with
relatively small energies or for those emitted from the inner core with much
larger temperatures.

One key approximation was a fixed redshift of emission $z=1$. As a next step
we consider the extreme opposite case provided by the broad redshift
distribution derived from the Dickinson and Madau (DM) star-formation rate
(bottom curve in Fig. 9) that has an average core-collapse redshift of 1.51.
The average photon energy in the earlier case was
$\langle\omega\rangle=(3/4)\,T=0.75\,T$, whereas in the DM case it hardly
changes to $0.71\,T$. The location of the maximum of
$\omega^{2}d\Phi_{\gamma}/d\omega$ changes even less from $\omega_{\rm
max}=T\,(1+\sqrt{3})/2=1.366\,T$ to $1.340\,T$. The maximum itself is 0.80 of
the earlier value, so finally we find

$\zeta_{a}\lesssim 0.54\times 10^{-4}\big{/}{n^{\rm cc}_{7}}.$ (76)

Using this very different $n_{\rm cc}(z)$ distribution causes only a minimal
modification, much smaller than many other uncertainties such as our rough
representation of the observational data. We conclude that in practice the
exact core-collapse redshift distribution is irrelevant as well as the exact
spectrum of emitted bosons. The main uncertainty derives from the integrated
core-collapse rate $n_{\rm cc}$ where the published range spans a factor of 2.

####  VI.3.2 Long-lived bosons

For bosons that decay more slowly—and this will be the case for some of the
interesting mass range for muonic bosons—we expand the decay fraction $f_{\rm
D}$ and use the predicted photon density of Eq. (66). Otherwise we can go
through the same steps where the main difference is that the result now
depends on the emission spectrum which sets the scale for the Lorentz factor
$m_{a}/E_{a}$ in the decay rate. Using the Dickinson-Madau version for the
redshift distribution we find

$\zeta_{a}\,\frac{m_{a}}{E_{\rm av}H_{0}\tau_{a}}<0.58\times
10^{-4}\big{/}{n^{\rm cc}_{7}}$ (77)

as a generic limit.

###  VI.4 ALPs

To make contact with the previous literature we can immediately apply these
results to ALPs. Our cold reference model emits a total energy in the form of
ALPs of $E_{{\rm SN},a}=G_{10}^{2}\,2.5\times 10^{49}~{}{\rm erg}$ [see text
above Eq. (55)] and $G_{10}^{2}\,2.8\times 10^{50}~{}{\rm erg}$ for the hot
model [see text above Eq. (56)]. Relative to our cosmological average SN with
$E_{\rm tot}=3\times 10^{53}~{}{\rm erg}$ this is a fraction of
$\zeta_{a}=G_{10}^{2}\,0.83\times 10^{-4}$ and $G_{10}^{2}\,9.3\times 10^{-4}$
respectively. With the limit Eq. (76) for short-lived bosons that uses the
Dickinson-Madau redshift distribution we find

$G_{a\gamma\gamma}<0.81\,(0.24)\times 10^{-10}~{}{\rm
GeV}^{-1}\Big{(}\frac{1}{{n^{\rm cc}_{7}}}\Big{)}^{1/2},$ (78)

if all cosmic SNe emit on average as many ALPs as our cold (hot) reference
model. This bound agrees well with Ref. [34] where
$G_{a\gamma\gamma}<0.5\times 10^{-10}~{}{\rm GeV}^{-1}$ was stated for
$m_{a}=5~{}{\rm keV}$.

For long-lived bosons we use instead Eq. (77) with the same ALP fractions from
the cold (hot) models. Moreover, the decay time is

$H_{0}\tau_{a}=H_{0}\frac{64\pi}{G_{a\gamma\gamma}^{2}}\frac{1}{m_{a}^{3}}=\frac{3.02\times
10^{-5}}{G_{10}^{2}}\left(\frac{10~{}{\rm keV}}{m_{a}}\right)^{3}.$ (79)

With the average energies 89.9 (136.8) MeV we then find

$G_{a\gamma\gamma}<0.65\,(0.39)\times 10^{-11}~{}{\rm
GeV}^{-1}\Bigl{(}\frac{0.1~{}{\rm MeV}}{m_{a}}\Bigr{)}\Big{(}\frac{1}{{n^{\rm
cc}_{7}}}\Big{)}^{1/4}.$ (80)

For the cold model and ${n^{\rm cc}_{7}}=0.58$ these bounds are shown in gray
in Fig. 12. The short and long-lived cases cross at $m_{a}\simeq 10~{}{\rm
keV}$.

###  VI.5 Muonic bosons

####  VI.5.1 Simplified estimates

We now turn to our main case of interest, muonic (pseudo)scalars that are
emitted from SN cores by the photo production processes of Sec. V.4 and then
of course decay by the loop-induced two-photon coupling. The analysis is fully
analogous, except that the SN emission now derives from photoproduction on
muons. We will go through the steps for both our cold and hot Garching
reference models as if they were an average cosmic core collapse to obtain a
generous range of possibilities. While the main uncertainty depends on which
model is taken to be representative, we first repeat the previous study on the
impact of the $n_{\rm cc}$ redshift distribution as well as the detailed boson
emission spectrum.

So beginning with our simplified estimate once more we assume $n_{\rm
cc}^{\prime}=n_{\rm cc}\delta(z-1)$ and a quasi-thermal boson emission
separately for the cold and hot SN model. For the short-lived case we write
the analog of Eq. (73), but using the energy release $E_{\rm
SN}=N_{a,\phi}\times\langle E_{a,\phi}\rangle$, where the total number of
emitted bosons and their average energy are given in Tab. 3. Comparing with
the measured flux in Eq. (72) the limit on the Yukawa coupling is then

$g_{\rm a,\phi}<2.8\times 10^{27}\Big{(}\frac{1}{{n^{\rm
cc}_{7}}}\Big{)}^{1/2}\Big{(}\frac{1}{N_{a,\phi}}\Big{)}^{1/2}\Big{(}\frac{\rm
MeV}{\langle E_{a,\phi}\rangle}\Big{)}^{1/2},$ (81)

and therefore, using the values reported in Tab. 3, the bounds from the cold
(hot) model are

$\displaystyle g_{\phi}$ $\displaystyle<$ $\displaystyle 0.32\,(0.11)\times
10^{-10}\Big{(}\frac{1}{{n^{\rm cc}_{7}}}\Big{)}^{1/2},$ (82a) $\displaystyle
g_{a}$ $\displaystyle<$ $\displaystyle 0.72\,(0.21)\times
10^{-10}\Big{(}\frac{1}{{n^{\rm cc}_{7}}}\Big{)}^{1/2}.$ (82b)

As expected they are indeed roughly two orders of magnitude stronger than
those from the SN 1987A cooling argument.

Now we study the opposite limit of long lived bosons. In this case, following
the steps of Sec. VI.3.2, the condition on the Yukawa coupling are found to be

$\displaystyle\kern-10.00002ptg_{\phi}$ $\displaystyle<$ $\displaystyle
1.11\times 10^{7}\Big{(}\frac{0.1\,\rm
MeV}{m_{\phi}}\Big{)}\Big{(}\frac{1}{N_{\phi}}\Big{)}^{1/4}\Big{(}\frac{1}{{n^{\rm
cc}_{7}}}\Big{)}^{1/4}\\!,$ (83a) $\displaystyle\kern-10.00002ptg_{a}$
$\displaystyle<$ $\displaystyle 0.91\times 10^{7}\Big{(}\frac{0.1\,\rm
MeV}{m_{a}}\Big{)}\Big{(}\frac{1}{N_{a}}\Big{)}^{1/4}\Big{(}\frac{1}{{n^{\rm
cc}_{7}}}\Big{)}^{1/4}\\!.$ (83b)

Therefore using the $N_{\phi,a}$ values reported in Tab. 3 the limits for the
cold (hot) model are

$\displaystyle\kern-10.00002ptg_{\phi}$ $\displaystyle<$ $\displaystyle
0.34\,(0.21)\times 10^{-10}\,\Big{(}\frac{0.1\,\rm
MeV}{m_{\phi}}\Big{)}\Big{(}\frac{1}{{n^{\rm cc}_{7}}}\Big{)}^{1/4}\\!,$ (84a)
$\displaystyle\kern-10.00002ptg_{a}$ $\displaystyle<$ $\displaystyle
0.46\,(0.27)\times 10^{-10}\Big{(}\frac{0.1\,\rm
MeV}{m_{a}}\Big{)}\Big{(}\frac{1}{{n^{\rm cc}_{7}}}\Big{)}^{1/4}\\!.$ (84b)

####  VI.5.2 Full numerical distributions

In order to verify the quality of these simple approximations, we now proceed
to directly evaluate Eq. (62) in full generality, using the numerical boson
fluxes together with the complete cosmic core-collapse rate $n^{\prime}_{\rm
cc}(z)$. To show the maximum plausible effect we use the $n^{\prime}_{\rm cc}$
of Madau and Dickinson who provided a total core collapse rate of ${n^{\rm
cc}_{7}}=0.58$.

In Fig. 11 we show the resulting limits for the scalar (solid red curve) and
pseudoscalar case (solid black curve). In the same plot the dotted curves are
instead the results from the simple analytical recipe of Sec.VI.5.1, rescaled
with the corresponding total number of past core collapses ${n^{\rm
cc}_{7}}=0.58$. We see that the agreement in the two opposite limits, short-
lived and long-lived bosons, is excellent, with a smooth and fast transition
between the two around $m_{a,\phi}\simeq 0.1$ MeV.

Figure 11: Bound on the Yukawa coupling for the scalar (red curves) and
pseudoscalar (black curves) from the diffuse SN $\gamma$ rays. The solid
curves are obtained numerically computing Eq. (62) for the cold SN model and
the cosmic core-collapse rate according to Madau and Dickinson [87]. The
dotted curves instead are obtained with the analytical prescription of
Sec.VI.5.1, normalized to a total number of past core collapses ${n^{\rm
cc}_{7}}=0.58$.

It was already clear in the ALP case and we confirm once more that the
detailed energy distributions of the emitted bosons as well as the exact
redshift distribution of the cosmic collapses modify the limit on the 10–20%
level. The real uncertain parameters of interest are the total core collapse
rate $n_{\rm cc}$ and the coupling-strength dependent total energy emitted by
an average cosmic SN. In our summary plot we have used the coldest Garching
model to be conservative although we are afraid that we may be taking
conservatism to extremes because the cosmic distribution must involve a range
of progenitor masses. Deriving a plausible cosmic average of this type would
be an interesting exercise along the lines of what has been done for the DSNB
[79].

##  VII From Colliders to Cosmology

We now briefly comment on other bounds of interest for the considered range of
masses and couplings. In particular, we discuss bounds from cosmology and
colliders which are relevant for the strong coupling regime in the
neighborhood of the $g_{\mu}{-2}$ inspired values.

###  VII.1 Cosmology

Cosmology is a powerful tool to constrain new light particles. In particular,
they can change the cosmic radiation density that normally consists of the
cosmic microwave background (CMB) at a present-day $T=2,726$ K and the cosmic
neutrino background (CNB) with $T_{\nu}=(4/11)^{1/3}T_{\rm CMB}=1.946~{}{\rm
K}$, although neutrinos today are nonrelativistic and contribute to dark
matter, not to radiation.

The cosmic radiation density is traditionally expressed in terms of $N_{\rm
eff}$, the effective number of neutrino species, defined by the CMB radiation
density today [89],

$\rho_{\rm rad}=\Big{[}1+\frac{7}{8}\Big{(}\frac{4}{11}\Big{)}^{4/3}N_{\rm
eff}\Big{]}\rho_{\rm CMB}.$ (85)

The standard-model prediction is $N_{\rm eff}^{\rm SM}=3.045$ [90, 91, 92, 93]
where the difference from 3 derives from small deviations from equilibrium
when neutrinos freeze out. The Planck collaboration recently reported, within
the framework of the standard $\Lambda$CDM cosmology, the restrictive range
$N_{\rm eff}=2.99\pm 0.34$ at $95\%$ CL [94].

While new particles usually add to $N_{\rm eff}$, our case is different and
reduces it. Early on, the new bosons are in equilibrium with muons, providing
more radiation. However, if they decay radiatively after neutrino decoupling,
they will heat the photons so that later the CNB will yet colder than the CMB,
an effect that reduces $N_{\rm eff}$ [18]. One additional thermal boson
increases the entropy degrees of freedom before the adiabatic disappearance of
electron-positron pairs and before boson decay to $2_{\gamma}+1_{\rm
boson}+\frac{7}{8}4_{e^{+}e^{-}}=\frac{13}{2}$. After $e^{+}e^{-}$ and boson
disappearance, we are left with the 2 photons, so the ratio between them is
now 4/13, instead of 4/11. To write Eq. (85) in the traditional form we need
to replace $(4/11)^{4/3}N_{\rm eff}$ with $(4/13)^{4/3}3.045$ and find $\Delta
N_{\rm eff}=-0.608$. This large negative deviation is strongly excluded by
Planck.

This powerful argument applies to bosons decaying after neutrino freeze out,
i.e. for masses $m_{a,\phi}\lesssim 2$ MeV, as shown in the upper left panel
in Fig. 3 of Ref. [18]121212See also Fig. 6 of Ref. [95], although in that
case the limit extends to slightly larger masses because the vector has three
degrees of freedom instead of one.. For larger masses the density of the new
boson is strongly Boltzmann suppressed by the time of neutrino decoupling,
making the new degrees of freedom harmless.

Finally, one can also ask for which range of the Yukawa coupling the above
argument holds. In fact, we have assumed the new boson to be in thermal
equilibrium with the standard model (SM) bath in the early universe, but for
small enough Yukawa couplings this will not be the case. In order to determine
the lower limit on the couplings, we can compare the interaction rate with the
Hubble expansion rate for $T\simeq m_{\mu}$. In a radiation dominated universe
the Hubble rate is [30]

$H(T)=1.66\sqrt{g_{\star}(T)}\frac{T^{2}}{M_{\rm Pl}},$ (86)

where $g_{\star}(T)$ are the relativistic degrees of freedom [of
$\mathcal{O}(10)$ for $T\simeq m_{\mu}$] and $M_{\rm Pl}=1.22\times 10^{22}$
MeV is the Planck mass. We then determine the Yukawa coupling such that

$n_{\mu}2\sigma_{\phi,a}\big{|}_{T=m_{\mu}}=H(m_{\mu}),$ (87)

where we use the semi-Compton unpolarized cross sections in Eq. (5) for
$\omega=3m_{\mu}$, and the factor of two comes from photon polarizations. We
then find for the (pseudo)scalar case that for $g_{\phi,a}<2.0\,(2.9)\times
10^{-8}$ thermal equilibrium is not reached and therefore no limits from
$N_{\rm eff}$ can be placed. While this is not a strict limit, in this range
of parameters the SN arguments strongly dominate. Therefore CMB bounds are
particularly relevant only for the strong coupling regime and for masses below
$2$ MeV.

For very small couplings, other interesting cosmological bounds come from BBN
and CMB spectral distortions [96, 18]. However, for the masses and lifetime of
interest for our $\gamma$-ray constraints, these cosmological bounds are
strongly model dependent, e.g. on the assumed temperature in the dark sector
and the reheating temperature of the universe. Of course also the CMB bound
from $N_{\rm eff}$ may be easily circumvented, for example with the addition
of new light degrees of freedom increasing $N_{\rm eff}$.