Axion-sourced fireballs from supernovae

New feebly interacting particles would emerge from a supernova core with 100-MeV-range energies and produce $\gamma$-rays by subsequent decays. These would contribute to the diffuse cosmic $\gamma$-ray background or would have shown up in the Solar Maximum Mission (SMM) satellite from SN~1987A. However, we show for the example of axion-like particles (ALPs) that, even at distances beyond the progenitor star, the decay photons may not escape, and can instead form a fireball, a plasma shell with $T\lesssim 1$ MeV. Thus, existing arguments do not exclude ALPs with few 10 MeV masses and a two-photon coupling of a few $10^{-10}~{\rm GeV}^{-1}$. However, the energy would have showed up in sub-MeV photons, which were not seen from SN 1987A in the Pioneer Venus Orbiter (PVO), closing again this new window. A careful re-assessment is required for other particles that were constrained in similar ways.


I. INTRODUCTION
The core of a collapsing star is one of the hottest and densest regions in the Universe. Therefore, it can be a factory of high-energy (around 100 MeV) particles beyond the standard model, in particular, feebly interacting particles (FIPs), such as sterile neutrinos, dark photons, new scalars, QCD axions and axionlike particles (ALPs), and many others. It comes as no surprise, nowadays, that-despite the small number of neutrino events detected at several neutrino experiments-supernova 1987A (SN 1987A) represented a bonanza for bounds on FIPs. A renowned example is the SN 1987A energy-loss bound [1][2][3], that limits the luminosity of a novel particle ϕ to be smaller than the neutrino luminosity, L ϕ ≲ L ν , evaluated at 1 second after the bounce [4,5]. While a large body of literature is dedicated to light particles, the last decade has seen an ever-growing interest in the high-mass part of FIP parameter space, above 1 keV, that could have an impact on cosmology [6][7][8][9], astrophysical transients [10,11], and play the role of dark matter mediator [12,13].
Heavy FIPs can be probed down to luminosities much smaller than L ν . Depending on their lifetime and decay channels, they can travel for distances that are either smaller or larger than the progenitor radius. If the mean-free path against decay is small enough, FIPs can decay in the mantle of the progenitor, lighting up the SN [14]. Since the explosion energy of SNe is much smaller than the energy released in neutrinos, the radiative decay of FIPs to charged leptons and photons must be suppressed [15,16]. This argument, applied to low energy SNe, limits the lifetime of radiatively decaying particles even further [17]. If the mean-free path against the decay is larger (or if the FIP decays to neutrinos), then the decay product of FIPs produced by SN 1987A could have been seen directly in the detectors online during the explosion [18][19][20][21]. For example, the decay to 100-MeV neutrinos would have been seen in Cherenkov detectors, so the lack of such events gives constraints superseding the energy-loss bound by one order of magnitude in coupling [22].
The focus of this paper will be on heavy ALPs (hereafter axions) that can decay to photons through the coupling − 1 4 g aγγ aFF. The daughter photons would have contributed to the diffuse cosmic γ-ray background [16,23,24] and, for SN 1987A, would have showed up in the Gamma-Ray Spectrometer (GRS) on the Solar Maximum Mission satellite [25]. The latter measured the fluence in the interval 4.1-100 MeV during a 223.2 s interval coincident with the SN 1987A neutrino burst. Its observations have been used to constrain radiatively decaying particles by looking for γ-ray emission [16,21,26,27].
For such constraint to apply, the photons produced in the axion decay should have been in the energy interval detected by SMM [18]. If axions with a mass m a ≳ 1 MeV were produced in SN 1987A, then after escaping from the progenitor they would have decayed to γ-ray photons with energy around 100 MeV. If the density of such photons was sufficiently large, then they would rapidly produce a fireball, creating a plasma of electrons, positrons, and photons, as sketched in Fig. 1. Several processes rapidly drive the plasma to a much lower temperature. After this rapid thermalization, the plasma would follow the evolution of a "standard" fireball [28,29], though with a much smaller initial temperature and, potentially, negligible baryon load. 1 The gas would first expand adiabatically, converting the temperature into bulk momentum, and then expand freely. The photons produced in axion decays from SN 1987A would have had an energy E γ ≲ 1 MeV. Most of these photons would not have been energetic enough to be detected by SMM, though we are still able to put constraints using the data of the Pioneer Venus Orbiter Satellite (PVO) [31], which had an energy window 0.2-2 MeV. Our main results are collected in Fig. 2. Our new PVO bound covers the region of the parameter space carved by the fireball formation and previously thought to be excluded by SMM.
The paper is structured as follows. In Sec. II we describe the first stages of the fireball, from its formation to thermalization. In Sec. III we analyze the period of expansion. Sections IV and V are dedicated to the impact of the fireball formation on respectively the SN 1987A and diffuse γ-ray bounds, as well as to the new bounds we place on the axion parameter space. Finally, we devote Sec. VI to our conclusions.

II. FORMATION OF THE FIREBALL
Axions are produced in the protoneutron star (PNS) at the center of the supernova, mainly via Primakoff emission and photon-photon coalescence. Their spectrum is parametrized in Refs. [21,27]. Both papers account for Primakoff, and neglect coalescence, which is reasonable at light masses. Since we want to compare the region of fireball formation with the bounds drawn in Refs. [21,27], we self-consistently use their expressions; coalescence would increase the total number of photons injected and enhance the possibility of fireball formation. Further, we are mostly interested in masses below 60 MeV, as confirmed by our results, where coalescence should have little impact (see, e.g., the Supplemental Material of Ref. [17]).
The axion spectrum is extracted from Ref. [27] as The numerical values for the parameters C 2 , κ s , T eff are all taken from Ref. [27]; the expression for the Primakoff cross section for massive particles, σ 0 ðE a ; g aγ ; κ s ; m a Þ, is taken from Ref. [21]. Only those axions that decay outside of the progenitor, with a radius of approximately R ¼ 3 × 10 12 cm, contribute to the photons that are detected at Earth. Therefore, the total energy injected by the axions is where the decay length is We can similarly determine the total number of axions injected N and the total radial momentum injected P. Notice that SN 1987A had an atypical small radius. For a larger progenitor radius, as could be the case of a future supernova-for example, Betelgeuse has a radius about one order of magnitude larger than SN 1987A progenitor, Sanduleak-69 202 [32]-a smaller fraction of axions decay outside of the mantle.

A. Geometry of the fireball
The massive axions propagate with a speed close to the speed of light, approximately radially, and decay FIG. 1. Schematic representation (not to scale) of the fireball sourced by axions from a supernova. Axions are produced in the protoneutron star, travel outside of the progenitor star, and decay exterior to its surface. The produced γ rays, rather than conserving their spectrum, form a fireball.
substantially far from the center of the supernova. The decay is in reality a continuous process; however, most of the energy from the axion decay is injected at a distance of the order of the decay length. The latter is larger than the typical thickness of the axion shell propagating from the supernova at small couplings and masses. Based on these facts, we can schematically model the propagation as a shell of axions, with a thickness of the order of Δ 1 ≃ 3 s, the typical timescale over which they are emitted from the supernova.
The typical decay radius for axions that decay outside the progenitor is determined by averaging over both the energy distribution and the radius of decay distribution, e −ðr dec −RÞ=lðE a Þ dr dec =lðE a Þ, where r dec > R; this leads to the average radius For simplicity, here and henceforth we will indicate the average of a quantity x over the energy fluence of the axions decaying outside the progenitor as To estimate the thickness of the photon shell, we determine the quadratic dispersion in the radius of the ; this is the mean time at which the bulk of the axions outside the progenitor decay. If the axion is produced in the PNS at a time t i , and decays at time t dec , then the photon position at time t is The assumption that photons propagate radially is valid for axions with masses much belowĒ a ≃ 100 MeV, the typical axion energy. We will be interested in typical masses below 60 MeV, since at larger masses we find no fireball formation. Therefore, accounting for the angle with which the photons are emitted in the decay, typically of the order m a =Ē a , would account for 10% corrections to the average radius of photon emission. Without loss of generality, we shift the origin of time and define t i to be uniformly distributed between −Δ 1 =2 and Δ 1 =2, where the time window is given a representative value Δ 1 ¼ 3 s; with this choice the average value of t i vanishes. At a fixed energy, the decay time is averaged over the exponential distribution e −½t dec −R=vðE a Þ=τðE a Þ dt dec =τðE a Þ for t dec > R=vðE a Þ. The average position of a photon at time t therefore is FIG. 2. Region of fireball formation in the axion (ALP) mass and coupling space (solid black line). We also show bounds from lowenergy supernovae (red) [17], γ rays due to axion decays from SN 1987A (blue) [21,27], and diffuse supernova background of γ rays from axion decays (green) [16,17]. The constraints from low-energy supernovae are the only ones unaffected by the fireball formation. Bounds from γ-ray decay at the Solar Maximum Mission (SMM) and from the diffuse supernova background should be reevaluated in the fireball formation region. On the other hand, we find that the entire region, which we hatch in blue, is robustly excluded by the Pioneer Venus Orbiter observations. Evaluating this att, we find, as we should, our previous estimate for the fireball radius, We now compute the average of the squared radius and define the thickness as Expanding the averages over the exponential distribution, this is The first term corresponds to the average thickness of the axion shell, while the other terms correspond to the thickness induced by the delay between the slower axions yet to decay and the photons. Notice that for ultrarelativistic axions the last two terms become arbitrarily small, since axions and photons move with arbitrarily close velocities. Therefore, its order of magnitude for strongly boosted axions, as we can approximately consider for masses below 60 MeV that are of interest here, is rm 2 a =Ē 2 a . The thickness of the photon shell will therefore be determined by the largest among the thickness of the axion shell and the approximate spread rm 2 a =Ē 2 a . Notice that in the thickness of the shell the corrections due to the noncollinear emission of photons in axion decay might lead to corrections by factors of order 1, since these are suppressed by one factor of m 2 a =Ē 2 a and must be compared with terms of the same order of magnitude. On the other hand, we will find later that the thickness of the shell in the end never enters the final results presented here. Therefore, we stick to this simple estimate, which might be corrected by factors of order unity, since none of our results are affected by it.

B. Thermalization via pair production
If the photon density is large enough, then it can form a plasma by producing pairs of e AE . The pair production process (and subsequent Compton and Rutherford scattering) leads to kinetic equilibrium, with both γ and e AE acquiring thermal spectra. Since pair production is not able to change the total number of particles, chemical equilibrium cannot be fully reached, and both species develop a common temperature T i and (negative) chemical potential −μ. Furthermore, they behave as a tightly coupled fluid with a bulk Lorentz factor γ.
In the following, we will refer with primes to lab-frame quantities, and without primes to plasma-frame quantities. If the number density of photons n 0 0 ¼ 2N =4πr 2 Δ is too low, then they are not able to produce pairs of electrons and positrons. The condition for pair-production reactions to be efficient is where σ γγ→e þ e − is the pair production cross section, which is evaluated at the typical photon energies in the plasma frame, where the photons are isotropic. The cross section is therefore [33] where E γ is the center-of-mass energy of any of the two colliding photons. We average this cross section over a uniform distribution of the pitch angle between the two incoming photons, with all the photons having a constant energy m a =2 in the comoving frame. In reality, there will be photons from axions moving with different speeds with a relative angular distribution different from the isotropic one and with a relative boost with respect to one another; however, these small corrections typically change negligibly the shape or normalization of our curves corresponding to fireball formation due to the large power of the coupling appearing in these curves (see below). If Eq. (11) is satisfied, then a population of electrons and positrons is nearly immediately established in the plasma. Because the typical energies at this stage are much larger than the electron mass, the energy and number of particles are equally shared among the two populations of electrons and positrons and the population of photons. The corresponding fluid has equation of state p ¼ ρ=3, where p is the fluid pressure and ρ is the fluid energy density. Therefore, we may deduce its bulk velocity taking the ratio of the conservation of energy density (see, e.g., Ref. [34], Chapter 2) and momentum density Here, v is the bulk velocity and γ is the bulk Lorentz factor. After taking the ratio, we obtain We can also deduce the initial temperature of this plasma using the conservation of the total number of particles, which at this stage is valid and yields where n is the number density of the fluid. Since at this stage chemical equilibrium is not maintained, the chemical potential of photons, electrons, and positrons is equal to a common (negative) value Therefore, we can use the equation of state of a gas of relativistic Boltzmann particles: taking the ratio of Eqs. (14) and (16) we obtain the initial temperature of the plasma as

C. Thermalization via bremsstrahlung
With a sufficiently large concentration of e AE , bremsstrahlung reactions can become fast enough to drive the system towards chemical equilibrium, increasing the particle number. The formation of the fireball requires pair production to equilibrate, while the particle number dilution also entails bremsstrahlung. For both reactions to be fast, axions should be produced in large quantities, with a large fraction decaying outside the progenitor.
The condition for the initial population of electrons and positrons to be equilibrated by bremsstrahlung is where the factor 2=3 accounts for the fact that only the electron-positron part of the plasma interacts via bremsstrahlung and σ ee→eeγ ðT i Þ is the cross section for bremsstrahlung from ultrarelativistic electrons and positrons. The total bremsstrahlung cross section diverges due to the emission of soft photons with negligible energies that, however, do not contribute to the thermalization of the bulk of the plasma. We use as an estimate of the bremsstrahlung cross section the averaged energy emission rate per electron from Eq. (16) of Ref. [10], defined as where n e is the electron density, ρ e is the electron energy density, and dE γ =dt is the total energy emitted by bremsstrahlung per unit volume per unit time. For bremsstrahlung, different prescriptions are sometimes used in the literature, such as using the total cross section for emission of a photon with an energy larger than a typical electron energy; but this again leads to very similar results to the criterion used here. Figure 3 shows, for each axion mass, the minimum value of the coupling for which the energy density of photons is large enough to activate pair production and bremsstrahlung (the latter assuming, of course, that electrons and positrons are produced). For low axion masses, and correspondingly low average energies for the photons, pair production equilibrates at lower coupling than bremsstrahlung, as expected given that bremsstrahlung is suppressed by a higher power of the coupling. At large masses, as soon as electrons and positrons are produced, bremsstrahlung reactions become faster than pair production. The reason is that the average energy per particle, m a =2 in the rest frame, is so large that the pair annihilation cross section is suppressed as σ e þ e − →γγ ∼ α 2 =ðm a =2Þ 2 , whereas the bremsstrahlung cross section is σ ee→eeγ ∼ α 3 =m 2 e , up to logarithmic corrections. Therefore, for m a ∼ 2m e = ffiffiffi α p ∼ 10 MeV, bremsstrahlung is faster than pair production. The subsequent thermalization depends on whether pair production or bremsstrahlung is the faster process. The final state of the plasma is, however, independent of this. Nevertheless, we keep the discussion separate for these two cases for clarity. At large enough couplings and masses, the fireball does not form anymore; the reason is that the fraction of photons decaying outside the progenitor star becomes so low that pair production is not efficient. For this reason, we find no fireball formation for axions heavier than about 60 MeV.

D. Final state of the fireball
Since the total energy released outside of the progenitor, for a benchmark axion mass m a ¼ 40 MeV and coupling g aγγ ¼ 2 × 10 −10 GeV −1 , is E ¼ 1.22 × 10 50 erg, and the estimated volume of the shell is of order 4πr 2 Δ ¼ 3.60 × 10 38 cm 3 , thermalization tries to drive the plasma towards significantly low temperatures, of the order of T ≃ ð15E=4π 3 r 2 ΔÞ 1=4 ≃ 192 eV. 2 However, this state is never really reached. Thermalization via bremsstrahlung is interrupted much sooner, since at temperatures smaller than m e ¼ 0.51 MeV the electron population is depleted by the faster pair annihilation.
For m a < 10 MeV, in the intersection of the red and blue regions of Fig. 3, pair production is immediately the faster process. The role of pair production is therefore that of maintaining the chemical potential of electrons and photons equal to a common value μ γ ¼ μ e þ ¼ μ e − ¼ −jμj. The slower bremsstrahlung reaction (still faster than the timescale associated with the expansion of course) gradually drives this chemical potential towards 0, with the electrons and positrons emitting photons to dilute the excess average energy per particle. As the chemical potential drops, so does the temperature, since the number of particles is increasing and therefore the average energy per particle decreases. This emission of particles is interrupted when the electron number density is so low that bremsstrahlung cannot proceed any further. Therefore, the final state of the thermalization is a plasma of electrons and positrons with an equal chemical potential μ, with a temperature T and bulk motion γ determined by the initial energy and momentum released by the axions, and by the condition that bremsstrahlung is out of equilibrium γn e ðT; μÞv th σ ee→eeγ Δ ¼ 1: We use here the nonrelativistic expression for the bremsstrahlung cross section, since the number density of electrons and positrons starts to be kinematically suppressed only at temperatures lower than the electron mass. For the same reason, we also account for a factor v th corresponding to the typical thermal velocities of the electrons in the plasma frame, of the order of ffiffiffiffiffiffiffiffiffiffiffi ffi m e =T p . For the full rate of bremsstrahlung interaction, we use the expression [10] For m a > 10 MeV, in the intersection of the red and blue regions of Fig. 3, after the population of relativistic electrons and positrons is established, bremsstrahlung becomes immediately the faster process. Therefore, in this intermediate stage, the chemical potentials of electrons and photons are not bound to be equal. On the contrary, electrons and positrons start to radiate photons by bremsstrahlung rapidly in order to drive the chemical potential of the photons to zero. In doing this, they dilute the average energy and reduce the temperature of the plasma. Therefore, the initially equal chemical potential of electrons and photons now are unbalanced, with jμ γ j < jμ e þ j ¼ jμ e − j.
This state of affairs is interrupted once the temperature of the plasma drops below 5 MeV, at which point pair production becomes again the dominant reaction as we have seen above. Therefore, there is a second stage at which pair production rapidly drives the chemical potentials of electrons and photons, which had been unbalanced by bremsstrahlung, to a common value intermediate between the two, so now μ γ ¼ μ e þ ¼ μ e − ¼ −jμj. At this point, the subsequent evolution is identical to the case m a < 10 MeV: pair production is so rapid that it maintains the chemical potentials equal at all time, and bremsstrahlung slowly dilutes the number density by trying to bring μ to 0, and the chemical equilibration is interrupted once the bremsstrahlung freezes out, namely when the condition in Eq. (20) is reached.
Thus, the cases m a < 10 MeV and m a > 10 MeV differ only in their approach to thermalization, but the final state is determined by the same conditions: conservation of energy, conservation of radial momentum, and decoupling of bremsstrahlung. The first two conditions are enforced by simply determining the total energy and radial momentum density of the initial axion shell, and imposing it equal to the energy and radial momentum density of the relativistic fluid of photons and electrons. The three conditions determine the three quantities of the plasma immediately after thermalization, namely the Lorentz factor γ, the temperature in rest frame T, and the common chemical potential of photons and electrons μ=T.
The determination of the final state of the plasma is simplified by making two assumptions, both of which are verified by the final result. First, since the thermalization is interrupted because the electron density is suppressed, we can already guess that the plasma at the end of thermalization will be dominantly composed of photons. Therefore, in its energy density and pressure, we may neglect the electron-positron contribution and use again the equation of state of relativistic Boltzmann particles with p ¼ ρ=3. This means that Eqs. (13)-(15) are still satisfied.
The second assumption that we make is that, in the final state μ ≫ T; this implies that thermalization is interrupted much before chemical equilibrium is reached, so that μ does not lower significantly. With this assumption, photons behave as relativistic Boltzmann particles and electrons/ positrons behave as nonrelativistic Boltzmann particles. Taking the ratio of Eqs. (14) and (20), after replacing the cross section Eq. (21) and the corresponding energy density and pressure for photons and number density for electrons/ positrons, we find One may obtain an approximate solution to this equation in the limit T ≪ m e , which in turn requires the condition In this limit, the solution is approximately The conditions for applicability are not always verified in our region of interest, which is why for the results shown in the figures we numerically solve Eq. (22); however, Eq. (24) has some pedagogical value, since it shows that the final temperature is a small fraction of the electron mass of order 10% with weak dependence on the initial parameters of the fireball. One may compare this with the initial temperature of the plasma before bremsstrahlung sets it, in Eq. (17), which on the other hand is of the order of the initial average energy per particle in the comoving frame.
Finally, the chemical potential can be obtained from any of the equations Eqs. (13), (14), and (20); for example, using Eq. (14) we find Numerically, we find that this quantity is indeed much larger than 1 in all the region of fireball formation. Figure 4 shows the properties of the plasma formed after thermalization. The plasma moves with a bulk Lorentz factor that is larger at lighter masses due to the relativistic beaming in the photons emitted. The temperature reached is a fraction of the electron mass, with little dependence on the amount of energy injected; the chemical potential is larger close to the boundaries of the region, where the final state is farthest from the chemical equilibrium state. Consistently with our approximation scheme, the chemical potential never reaches values below about 8T. Therefore, at the end of the thermalization phase, the plasma is composed of a dominant population of photons and a subdominant population of e AE with a sub-MeV temperature.

III. FIREBALL EXPANSION
The subsequent evolution is driven by radial expansion. As long as pair annihilation remains fast, this expansion happens in the hydrodynamical regime. In the case of γ ≫ 1, this regime is well known from the analogous fireball formed in γ-ray bursts [28]. In our parameter space, we find γ to be always larger than 2.5, which means that the assumption of ultrarelativistic motion is reasonably satisfied; this approach is all the more justified since the radial expansion increases the bulk kinetic energy at the expense of the internal energy, leading to a rapid increase of γ. Furthermore, since the number of particles cannot change in this phase of expansion, the average energy per particle (in the laboratory frame) is not significantly affected by the details of the expansion. For γ ≫ 1, the thickness of the shell remains constant, since the fluid moves radially with the speed of light. The evolution of the thermodynamic quantities is ruled by the conservation laws of the entropy per particle, σ γ ðrÞ=n γ ðrÞ, where σ γ is the entropy density; the average energy per particle, γðrÞρ γ ðrÞ=n γ ðrÞ, with ρ γ being the rest-frame energy density; and the total number of photons, γn γ 4πr 2 Δ. Unprimed quantities refer to the frame comoving with the fluid. These conservation laws happen to be identical to the fireball equations when chemical equilibrium is maintained, as derived from the hydrodynamical equations, see, e.g., Ref. [28]. We are assuming the plasma to be dominated by photons so that the e AE component can be neglected. The conservation of entropy per particle implies that the ratio μ=T remains constant, since for photons σ γ =n γ is only a function of μ=T. In turn, the conservation of the average energy per particle implies the constancy of the product γT. Finally, the conservation of the photon number implies the scaling laws γ ∝ r and therefore T ∝ r −1 .
The expansion ends when pair annihilation runs out of equilibrium. At this point, photons moving at an angle θ with the radial direction will have a thermal distribution with an effective temperature where the subscript f means that the quantity must be evaluated at the moment of decoupling. Since the photons reaching Earth will mostly be forward, with an angle θ ≪ γ −1 , we may safely take the limit θ ≃ 0, and therefore photons reaching Earth will have an average temperature T 0 ≃ 2T f γ f for γ f ≫ 1. Using the constancy of Tγ during the hydrodynamic expansion phase, this can be evaluated at the beginning of the expansion phase, that is, immediately after thermalization. As it should be, this quantity is of the same order of the average energy per particle in the lab frame.
A. Remaining fraction of high-energy photons Figure 5 shows the average energy of the thermal distribution seen at Earth, namely three times the effective temperature, in the parameter space of interest. At low axion mass, the average energy can reach up to the order of MeV, mostly because of the higher bulk Lorentz factor due to the relativistic motion of the light axions. However, already at about 10 MeV, the average energy becomes sub-MeV. To qualitatively understand the impact on the bounds from the nonobservation of γ rays in Refs. [21,27], we look at the fraction of photons that, after thermalization, remain at an energy higher than 25 MeV, which is the threshold energy above which upper bounds on the fluence have been set. Therefore, we show in the middle panel of Fig. 5 the fraction of energy remaining in photons above E min ¼ 25 MeV, the minimum energy at which data were available from SN 1987A. Since the photon spectrum has a Boltzmann shape, namely proportional to E 2 e −E=2γT , this fraction is This fraction rapidly drops by more than four orders of magnitude at a mass larger than about 5 MeV. All the previous discussion relates to the bulk of the photons produced in the decay of the axion, at distances comparable with the mean decay length or the progenitor radius. However, at sufficiently large radii, the photons coming from the decay of the small fraction of surviving FIG. 5. Spectral properties of the γ-ray emission from the fireball. We show the average energy of the thermalized photons (left), and the fraction of the thermalized photons integrated above 25 MeV (middle). We also show the fraction of axion decaying far enough from the center that the produced high-energy photons cannot thermalize and reach the Earth with their original energies; the red line is the boundary of the region in which the pair production reaction is always kinematically allowed from the criterion E CM;min > m e . axions may not be able to thermalize with the dominant bath of lower energy photons, if the density of the photon shell has sufficiently rarefied. This small fraction of the total photon flux reaches Earth with the originally expected large energies.
To quantify the impact of this contribution, we estimate up to what radius the high-energy photons from late axion decays can thermalize with the bulk photon fluid. The dominant channel for thermalization is pair annihilation. While the temperature of the photon fluid is significantly lower than m e , the photons from axion decay have a typical lab-frame energyĒ a =2, withĒ a ∼ 100 MeV. In the rest frame of the decaying axion, the high-energy photons are isotropically distributed. Boosting to the laboratory frame, and then into the frame comoving with the plasma, we find that the center-of-mass energy for the collision of the two photons is where we have assumed an average energy 3T for the restframe photons in the plasma; here V ¼ is the typical axion velocity. As long as E CM > m e , photons from axion decay are kinematically allowed to collide with the photons from the bulk of the plasma. Notice that the minimum value that E 2 CM can take as the plasma expands and γ becoming larger is simply obtained by writing V ¼ 1 A sufficient condition for the reaction to be kinematically allowed is E CM;min > m e . The reaction will remain in equilibrium provided that where we estimate the pair production cross section as σ γγ→e þ e − ≃ πα 2 =s, with s the Mandelstam parameter of the collision. This estimate holds up to factors of order unity. Once the radius of the fireball has expanded by a factor x to a radius r pair ¼ xr, the condition for the freeze-out of the pair production reaction becomes where we identify by the suffix 0 the quantities at the beginning of the expansion. Expanding we can obtain the factor of expansion as This defines the maximum radius r pair ¼ xr out of which injected photons are not able to thermalize with the bulk of the plasma. At this radius, the fraction of energy injected is exponentially suppressed as We show this suppression factor in the right panel of Fig. 5: this should be interpreted as the fraction of energy produced by axions at a large enough radius that the decay photons are not able to thermalize with the remaining bath of low-energy photons. The red boundary identifies the region within which E CM;min > m e , and therefore our treatment above is applicable. We define the fraction of energy that is still injected above the photon energy of 4 MeV as f ¼ max ½f pair ; f spectrum .

IV. UPDATED BOUNDS FROM SN 1987A
Since the current bounds from SN 1987A rely on the energy window above 25 MeV [21,27], fireball formation can impact them. Unless the axion is very light and the photon energy remains large due to the large Lorentz factor, fireball formation reduces the photon flux above 25 MeV by more than ten orders of magnitude. (We comment on SMM data in the interval 4-25 MeV in the Appendix.) This already indicates that a region of the parameter space, previously excluded by nonobservation by SMM of the decay photons, may be ruled in-though, as we will see, this is not the case.
Thus, our first result is the region identified by the black solid line in Fig. 2, in which the bounds from both γ-ray decay from SN 1987A and diffuse supernova background need to be reevaluated, since the spectrum adopted in previous literature for γ rays should actually be replaced by a Boltzmann spectrum with a temperature 2γT. Notice that this region is highly complementary to bounds from energy deposition in low-energy supernovae. If axions produced in the PNS travel and decay in the mantle, the ejecta kinetic energy becomes too large. Assuming that less than 0.1 B [1 BðbetheÞ ¼ 10 51 erg] energy is deposited, as observed in low-energy supernovae, the solid red bounds are obtained [17]. A particularly conservative limit of 1 B leads to the dotted red bound. Since these bounds are calorimetric, they are not strongly affected-if at all-by the formation of a fireball. We notice that even within the red region, axion production has negligible impact on the inner dynamics of the SN core, where the axions themselves are produced. The cooling rate due to axion emission in this region of the parameter space is much lower than the one due to neutrino emission. Therefore, the axion flux we are considering is self-consistent.
While the limits from SMM do not apply straightforwardly anymore, we find that the entire region of fireball formation is directly excluded by the nonobservation of γ rays at the Pioneer Venus Orbiter (PVO). Launched in 1978 (last contact in 1992), PVO was part of NASA'S Pioneer Venus project and featured a γ-ray burst detector with two NaI photomultiplier detector units [35,36]. The sensitivity window of this detector was in a much lower energy range, between 0.2 and 2 MeV, which is just the region in which the average energy of the thermalized spectrum lies in the region of fireball formation. Moreover, PVO had a 4π acceptance and, being in orbit around Venus, it was obviously outside of the Earth's radiation belts, so the background was particularly low [37].
As reported by Ref. [37], at the time of SN 1987A no excess over the background was observed in any of the energy bins observed by the instrument. While the response functions of the experiment are not published in this paper, we may extract a coarse bound on the fluence from Ref. [38], which reports the detection of γ-ray bursts in the same energy window with a fluence of the order of about 10 −4 erg cm −2 . At the distance of SN 1987A of about 50 kpc, this means that a minimum total energy E of about 3 × 10 43 erg should have been detected. In the whole region of fireball formation, we find that the total energy E is always larger than 5 × 10 48 erg, namely more than five orders of magnitude higher than the approximate exclusion. Therefore, even without refining the calculation, we can already claim that the entire region of fireball formation is excluded on the basis of the PVO nonobservation of a γ-ray burst in the low energy range 0.2-2 MeV. In fact, the bounds from PVO may extend even outside the region of fireball formation, since the γ-ray spectrum from axion decay is more or less flat so that a non-negligible fraction of energy would be released in the low-energy window even without fireball formation; however, since in these regions bounds would be only complementary to the ones from SMM, we do not investigate this question further.
It is logical to ask if other detectors could probe the low-energy flux of the fireball from SN 1987A. 3 While the GRS onboard of the SMM could detect photons even down to tens of keV [40], no available public data exist of SN 1987A in this energy band. Other hard x-ray detectors that observed SN 1987A, such as a Jet Propulsion Laboratory balloon [41] and the Ginga satellite [42,43], did so only days after the collapse.

V. UPDATED BOUNDS FROM DIFFUSE GAMMA-RAY BACKGROUND
Photons originating from the decay of heavy axions produced by all past SNe would contribute to the diffuse cosmic γ-ray background. The idea of constraining neutrino radiative decays in this way dates back to the 1970s [44], and the argument has been revisited recently with applications to axions (see, e.g., [16,24]).
We here follow closely Ref. [16] to compute the expected photon flux ignoring the possibility of a fireball. Assuming that all SNe occur at z ¼ 1, the photon energy flux (energy per unit area, time, and solid angle) due to axion decays is where ζ a is a fudge factor proportional to g 2 aγγ , E SN ¼ 3 × 10 53 erg is the typical energy released by a SN, T eff is the effective temperature akin to the one appearing in Eq. (1), E av ¼ 3T eff is the axion average energy, and n cc ≃ 10 7 Mpc −3 is the core-collapse number per comoving volume per redshift interval. This spectrum has a maximum at ω max ¼ T eff ð1 þ ffiffi ffi 3 p Þ=2 ≃ 1.37T eff . The diffuse γ flux produced by axions should be compared with the measurements of the extragalactic background radiation [23]. For axions with masses and couplings not allowing the formation of a fireball, we can use the range of the extragalactic background spectrum 2-200 MeV [23], Therefore, the bound is found comparing Eq. (33) evaluated at ω max , with Eq. (34). One obtains [16] ζ a ≲ 0.43 × 10 −4 =n cc ; ð36Þ independently of the assumed average energy of the emitted bosons. This corresponds to a constraint on axions with a photon coupling g aγγ ≲ 10 −10 GeV −1 [16,17]. The result does not depend strongly on the assumed redshift dependence of the cosmic core-collapse rate n cc . The formation of the fireball impacts the diffuse γ-ray bounds: as the average energy of photons is driven below 1 MeV, the flux from axion decays needs to be compared to a much larger astrophysical background (see, e.g., Fig. 10 of Ref. [23]). A very conservative estimate can be obtained using the largest value of the extragalactic background light energy flux, This is 25 times larger than the value that is observed in the 2-200 MeV energy range. Therefore, the diffuse γ-ray bound is relaxed at most by a factor of 5 (g aγγ ≲ 5 × 10 −10 GeV −1 ) in the region corresponding to fireball formation. On the other hand, the bound rapidly degrades at m a ≃ 50 MeV [17], so we cannot commit to a precise evaluation of this region, and a dedicated analysis is needed to revisit this bound at large masses.

VI. DISCUSSION AND OUTLOOK
Heavy axions coupling to photons can be produced in the hot core of supernovae and decay back to photons with around 100 MeV. While it has always been assumed that such photons would travel freely, in part of the axion parameter space photons form a fireball, and their energy degrades to values below 1 MeV. Part of the parameter space previously thought to be excluded by observations of SN 1987A with SMM [21,27] is actually excluded by PVO data [31,37], as the energy of the photons at Earth would have been much lower than expected. PVO observations, previously applied only to neutrino radiative decays [37], are relevant since cosmological constraints apply in the region of interest only for a large reheating temperature [7,9]. Constraints arising from the extragalactic background light (see, e.g., Ref. [16] and references therein) also get relaxed. Other complementary probes in the regions might be, for example, very low luminosity SNe and their light-curve shape and spectral line velocities, that can potentially probe unexplored parts of the parameter space [17]. Another interesting probe could rely on past [45] and future neutron star merger observations, which we plan to explore in forthcoming work. We stress that our discussion applies to different models, such as axions coupling to charged leptons [16,26] and heavy neutral leptons featuring a dipole portal [46]. Therefore, constraints arising from SN γ-ray observations should be updated to include PVO data. Fireball formation is even easier for axions with an electron coupling, since they can decay to electron-positron pairs and bremsstrahlung can drive the thermalization without the need of producing pairs through two-photon annihilation.

ACKNOWLEDGMENTS
We thank Hans-Thomas Janka, Georg Raffelt, and Irene Tamborra for comments on the first draft of this paper.

APPENDIX: BEST REACH OF SOLAR MAXIMUM MISSION
In the main text we have discussed what is the fraction of photons remaining above 25 MeV, since Refs. [21,27] both use γ-ray data above this energy. In reality, SMM did have additional observations down to about 4 MeV. Therefore, it makes sense to also look at the fraction of photons remaining above 4 MeV; this would provide a guide as to how powerful the bounds from SMM could be by accounting also for these additional data which were not used in past approaches. Figure 6 shows this fraction of energy. We find this to be significantly higher than the fraction of energy injected FIG. 6. Fraction of the thermalized photons integrated above 4 MeV. above 25 MeV, which is expected since the tail of the spectrum decreases exponentially. Therefore, SMM may in principle still constrain a significant part of the region of fireball formation, though with different data and a different approach than what has been done in the past. However, we do not follow up on this question further, since the bounds from PVO robustly exclude the entirety of the region.