Decay photons from the ALP burst of type-II supernovae

We determine limits from SN 1987A on massive axion-like particles (ALPs) with masses in the 10 keV - 100 MeV range and purely coupled to two photons. ALPs produced in the core collapse escape from the star and decay into photons that can be observed as a delayed and diffuse burst. We discuss the time and angular distribution of such a signal. Looking into the future we also estimate the possible improvements if the red supergiant Betelgeuse explodes in a supernova event.


I. INTRODUCTION
Many beyond the Standard Model scenarios include new massive (pseudo-)scalarsdubbed axion-like particles (ALPs) -among their particle spectrum (see e.g. Refs. [1][2][3] and Refs. [4][5][6] for recent reviews). The name originates from their similarity to the axion of the Peccei-Quinn solution to the strong CP problem [7][8][9]. Contrary to the usual QCD axion, that also couples to gluons, ALPs may solely interact with two photons, where a denotes the ALP and g aγγ is its coupling constant with dimension of inverse energy, which is often linked to an underlying scale of new physics f a via g aγγ ∼ α 2π 1 fa . In this paper we will focus on a pure coupling to photons 1 as given by Eq. (1). In contrast to the QCD axion, in the more general case of ALPs there is no fixed relation between mass and coupling: these are taken as completely independent parameters.
The coupling to photons allows the production of ALPs in stellar cores via the Primakoff mechanism. Even if this coupling is extremely small, a sizeable amount of ALPs can be produced in the stellar bulk and have observable consequences. Here our aim is to study very weakly coupled ALPs produced in type-II supernovae (SN) and the effects of their decay into photons outside of the star. The objective is to determine which regions in the g aγγ −m a space are allowed (or not) by the exquisite sensitivity that gamma-ray satellites can have to measure this photon burst. As we will detail in Section II, this process is most effective in the 10 keV -100 MeV mass range, where masses are smaller than the temperature of the SN core but large enough as to make the decay rate outside the star sizeable.
Our main result is the excluded region labelled "SN decay" in the ALP parameter space of Fig. 1. In this region our arguments provide better limits than existing laboratory and 1 The main reason for this is simplicity, but also that this is an often considered test model. Nevertheless, other interactions -in particular new decay modes -could be included by an appropriate recasting using an adapted decay length and branching ratio into two photons along similar lines as discussed in Ref. [10].  [6,[11][12][13] with added limits from Refs. [14][15][16][17][18][19][20]). Our bound is shown in dark blue ("SN decay").
astrophysical constraints. ALPs in this parameter range can have a strong impact in cosmology if their relic density is as large as the thermal one, but this depends on the maximum temperature of the big bang [14] which is, at the moment, unknown (cf. Section V).
We focus on SN 1987A, which has already been exploited to derive a variety of limits on ALPs. Perhaps the simplest one arises from the energy loss implied by significant ALP emission, which would reduce the measured neutrino burst below the ∼ 10 s observed by neutrino detectors [21,22] (light green region labelled SN 1987A in Fig. 1). For very light ALPs with masses below m a < few × 10 −10 eV a better limit can be obtained by taking into account that ALPs emitted from the supernova can convert into photons in the magnetic field of the galaxy [23,24], but no gamma-ray signal was ever detected after SN 1987A [18,[25][26][27][28][29] (dark green region labelled SN 1987A) 2 . For heavier ALPs this does not work because the 2 For a future supernova the sensitivity could be improved employing Fermi-LAT [30].
reconversion into photons is strongly suppressed.
For sufficiently heavy ALPs with masses in the 10 keV -100 MeV region however, another process becomes possible: the decay into two photons. This possibility was analysed in Ref. [10] assuming that ALPs are produced in the SN core via a direct ALP coupling to nucleons. In this paper we consider a less model-dependent case in which the photon coupling, Eq. (1), is responsible for both ALP emission and decay. Moreover, we improve on several aspects of the calculation 3 . It is very important to notice that not all photons can reach Earth. We improve the estimate of the number that do so by employing a numerical simulation, in particular carefully treating the time delay and the angular distributions.
In addition to SN 1987A we also consider possible future SNe events. For concreteness, and also because it may produce spectacular effects, we entertain the core collapse of the red supergiant Betelgeuse. This is particularly interesting since its distance to Earth is only ∼ 200 pc (∼ 650 ly), therefore much closer than SN 1987A (at 51.4 kpc ∼ 170000 ly).
The paper is organised as follows. In Section II we discuss the essentials of the production mechanism, the subsequent decay and the geometrical features relevant for our analysis. In Section III we briefly describe our numerical simulation and the obtained results for the detectable fraction of signal photons as well as their distribution in time and angle. In Section IV we use SN 1987A to obtain concrete limits and discuss the potential sensitivity if Betelgeuse goes supernova. We summarize and conclude in Section V.

II. SETTING UP THE ANALYSIS
ALPs are produced in the core of the SN via the Primakoff process, γ + p → p + a, whereby a thermal photon is converted into a pseudo-scalar in the presence of the external electromagnetic field provided by the charged particles in the medium (plasma). The typical energies of the produced ALPs are of the order of the core temperature and are in the 3 In Section III A we briefly discuss the discrepancies between our results and those from Ref. [10]) ∼ 100 MeV range. Recently the associated ALP energy spectrum has been calculated with detailed account of the production process in core-collapse SNe in Ref. [18] and we will use the ensuing ALP production rate to estimate the ALP-originated photon fluence, i.e., the number of gamma-ray photons per unit area, on Earth.
Despite the core being extremely dense (∼ 10 14 g/cm 3 ), due to the smallness of the coupling to two photons, the ALPs escape the SN essentially unimpeded (cf. Eq. (15)) and are emitted isotropically. As already mentioned, this is an effective energy-loss mechanism that in itself constrains the parameter space of ALPs [21,22]. However, the couplings we are interested in are so weak that the energy taken by the emitted ALPs does not affect the core collapse. Rather, they contribute a tiny perturbation to the standard picture of the collapse.
Nevertheless, our ALPs can decay outside the SN producing 10 MeV photons, which is much larger than the typical photon energy in the outskirts of the SN. These gamma rays have no standard background, so we are sensitive to a much smaller flux.
The ALP decay rate is [31] which results in a decay length for the ALPs given by and, for the masses and couplings we are interested in, the decay length is large, but still allows for a sizeable number of ALPs to decay between SN 1987A and Earth (d SN ∼ 51 kpc).
Therefore, ALP decay is the relevant gamma-ray production mechanism to be considered and we may expect a flux of ALP-originated photons on Earth. Once we have the ALP production rate we may convolute it with the decay probability to obtain the fluence of ALPoriginated photons at the detector. By obtaining (upper) limits on the gamma-ray fluence shortly after the observation of SN 1987A, we may constrain g aγγ and m a by demanding that the number of ALP-originated photons arriving at the detector does not exceed what was observed [32,33].
Another important point to be considered is that, having significant masses and being typically produced with energies ∼ few × 10 MeV, ALPs emitted from the SN core have appreciable -but not enormous -Lorentz boost factors (γ = E a /m a ). This has important consequences. Firstly, the angle between the two photons and hence also between the original propagation direction α of the parent ALP is non-vanishing (see Eq. (23)), thus implying that the ALP-originated photons that reach Earth from the SN are not necessarily emitted along the SN-Earth direction, but rather at an angle as schematically shown in Fig. 2. Conversely, this implies that on Earth we would see the photons as if they were coming from a direction somewhat off the location of the SN, i.e., the signal is effectively smeared out over a halo.
We must also consider that the length travelled by the ALP before decay plus the length travelled by the decay photon reaching Earth is longer than the distance between Earth and the SN, d SN [10] (cf. Fig. 2). This leads to a time delay that, considering the distances involved, can be of the order of years even if the decay angles are not very large. The fact that massive ALPs have a velocity < c increases the delay even further. Instead of a signal that lasts as long as the duration of the SN explosion and associated neutrino burst (∼ 10 s), as was the case for nearly massless ALPs in Refs. [18,[25][26][27][28][29] for SN 1987A, the signal from massive ALPs may be spread out over much longer time scales.

A. Flux of massive ALPs
The subtleties of the core collapse of a progenitor star of mass around ∼ 18M and the associated ALP production in its interior have been thoroughly analysed in Ref. [18]. We adopt their results for the production rate of massless ALPs, which, for g aγγ = 10 −10 GeV −1 , can be fitted by For our purposes, it is enough to consider the ALP production as instantaneous. We fit the integrated emission with the formula where T is an effective temperature, C contains information on the volume and density of scatterers and σ 0 (E a ) is the Primakoff production cross section for massless ALPs off non-relativistic targets given by Here k s is an effective Debye screening scale that determines the range of the Coulomb potential created by the scattering centers in the plasma 4 . We find that Eq. Both results are compared in Fig. 3, which shows excellent agreement. This gives us confidence that our fit values are close to the relevant physical time and space averages.
The production rate for ALPs depends only on the properties of the plasma, i.e. the temperature and the Debye scale. This is true both for the massless, as well as for the massive case. In a real supernova both depend on time and location inside the core of the progenitor star. However, as we have just discussed, the full production in the massless case can be very well approximated by assuming a suitable effective temperature and Debye scale.
We will therefore use this effective temperature and Debey scale to perform an extrapolation to the massive case. While this is certainly only an estimate, we believe it to be fairly accurate in light of the excellent fit for the massless case.
We can now calculate the ALP production by including the mass in the Primakoff cross section σ 0 , which is found to be (cf. Ref. [14]) where β = E 2 a − m 2 a /E a is the ALP velocity in natural units. The effect of non-vanishing ALP masses can now be encoded in an energy-averaged suppression factor, which we define as and show in Fig. 4 below. In the formula and figure we neglect the energy dependence for simplicity. However, in our numerical simulation this energy dependence is included. Note where we expect a suppression of the massive ALP production relative to the massless case.
We can compute the ALP flux by rescaling the massless flux with this factor.
where d SN is the SN-Earth distance and a factor of two is included to account for the two photons emitted per ALP decay. Here dN a /dE a , the ALP spectrum, is the result of integrating the production rate over the time of the core collapse (∼ 10 s) [18].
For nearly massless ALPs with instant decay we then have a naive fluence However, as already indicated above, for massive ALPs with a finite decay length we have to take a variety of additional effects into account. Therefore we correct the fluence to where P total (m a , g aγγ ) = S(m a )P survival P decay P time P acceptance .
Equations (13) and (14) are only formal: the total probability and fluence cannot be factored in general, as photons with a given energy can originate from ALPs of different velocities depending on the decay angle. However, it shows the various effects that must be included in a realistic estimate of the sensitivity: • S(m a ) is the mass-dependent factor defined in Eq. (10) which corrects Eq. (7) for massive ALPs. It is responsible for the suppression of the production for heavy ALPs (m a 20 MeV).
• P survival gives the fraction of ALPs decaying outside the region effectively occupied by the progenitor star, since the photons originating inside this region may be absorbed or at least scattered.
• P decay takes into account that the ALP-originated photons from a decay at a distance greater than d SN typically do not reach the detector. Photons emitted backwards with respect to the direction of the parent ALP could in principle reach the detector, even if the ALP decays after travelling beyond Earth, but their number is typically quite small so we neglect them in most of the discussion.
• P time is the fraction of ALP-originated photons arriving within the measurement time of the detector. These are the only photons that can be effectively counted.
• P acceptance accounts for the fact that some detectors may have a finite angular acceptance. Photons from ALP decays will arrive within a range of non-vanishing angles with respect to the SN (cf. Fig. 2). Therefore, a detector with finite angular acceptance will only see a fraction P acceptance of all photons. Besides this, detectors usually have specific energy ranges where their sensitivity is optimal, e.g. for SN 1987A we work with photons in the 25 − 100 MeV energy range at the detector [32,33]; this cut is also described by P acceptance .
In Section III we will numerically determine the effects of P total , but it is worthwhile to briefly address the essential factors. To simplify the discussion we consider next a situation of an ALP with fixed energy and discuss the probability for the resulting photons to reach the detector. The probabilities for such a case will be denoted by P . To obtain the energyaveraged P total in Eq. (12) the contributing factors have to be appropriately folded with the energy distribution (spectrum), which we chose to do numerically with a Monte Carlo simulation.
Since we already discussed S(m a ) earlier, let us comment on the second and third effects.
The second factor is simply the survival probability for the ALP to leave the volume of the progenitor star before decay, where R , the effective radius, is much larger than the actual radius of the progenitor's core itself (∼ 50 km for SN 1987A [18]). Following Ref. [34] we shall take for SN 1987A. As can be seen from Eq. (3), at large masses (above a few tens MeV) and couplings (above ∼ 10 −10 GeV −1 ), ALP is typically smaller than R so that, in this region, the sensitivity is strongly suppressed. For such large masses one may expect that the bound weakens due to the absorption of ALPs and this is indeed observed in our simulations (cf. Figs. 8 and 9). On the upper-right corner one sees that the region where ALP R * is not excluded. The bound behaves as g aγγ ∼ 1/m 2 a , which is compatible with Eq. (15) (cf. Eq. (3)). We have conservatively assumed that any ALP decaying inside the effective radius of the progenitor star will be absorbed, thus not leaving an observable trace. Incidentally, the fact that decays within R are practically blocked inside the progenitor places an upper limit on the mass range available to our analysis.
The effect of the third and fourth factors are somewhat entangled. Let us consider the typical time delay of an ALP-originated photon. Compared to a massless particle travelling directly the distance between the SN and Earth, we observe a time delay We are interested in the probability that the ALP decays before d SN and is detected within a given time frame δt (from data provided in Ref. [32] we can extract limits for time intervals δt 220 s). For concreteness, we consider small masses (m a MeV) and, in this region, we may estimate the fraction of events occurring within a measurement time δt ∼ 100 s as follows. To incur a time delay ∆t δt the ALP must decay before a distance As long as L max ALP , the probability of measuring ALPs with time delays ∆t δt is given by P decay × P time ≈ L max / ALP , i.e., which is the dominant effect limiting the sensitivity.
Having in mind that the ALP production cross section includes an extra factor of g 2 aγγ , we see that, for a given detection time, the ALP-originated fluence F γ = F naive γ × P total behaves as ∼ g 4 aγγ m 2 a , thus implying that the bound goes like which is exactly the behavior we observe in our numerical results (cf. Section IV).
Finally, since the detector for SN 1987A had a half-sky field of view [32] Fig. 3 we may estimate P acceptance as the normalized area under the curve within the aforementioned ALP energy range. By doing so we find that P acceptance ∼ 0.7. As we shall see in Section IV B, for Betelgeuse P acceptance will be reduced more significantly due to the energy range accessible to the Fermi-LAT detector that is slightly too high compared to the temperature of the supernova.
The angular and time distributions are closely related. Let us first note that, due to the assumed isotropy of the SN explosion, the angular and time-delay distributions will be the same at any point on a sphere with radius d SN around the SN, cf. Fig. 2. Therefore, to obtain the aforementioned distributions, it is enough to look at the distribution in angle and time with which the photons hit the surface of this sphere. In other words: if a photon originating from an ALP emitted in a certain direction hits this sphere, one can always find a rotation that puts Earth into the path of this particular photon. Hence, due to the isotropy assumption, emission of an ALP in this rotated direction has the same probability.
Let us consider an ALP emitted from the SN and decaying after covering a length L 1 .
One of the photons is emitted at an angle α relative to the direction of the parent ALP and then, after travelling a distance L 2 , hits the sphere of radius d SN . From Fig. 2 we see that these quantities are related via L 2 2 + 2 (L 1 cos α) L 2 + L 2 1 − d 2 SN = 0, which can be solved for the distance travelled by the photon from the point at which the ALP decays, Moreover, from the law of sines one finds that the incidence (detection) angle of the ALPoriginated photon with respect to the line of sight is given by 5 It is important to differentiate between two regions in space: 0 ≤ L 1 ≤ d SN and L 1 > d SN .
For the former, corresponding to an ALP decaying between the SN and Earth, it is clear that, for α ≤ π/2, only the plus sign is physically meaningful. For an obtuse decay angle, due to the condition L 1 /d SN ≤ 1, the plus sign is again the only choice. Both situations imply an incident photon in the "frontal" hemisphere of the detector, which is assumed to be aimed directly at the SN. 5 For a large number of photons with this detection angle the angular halo is ∆φ 2θ.
For the outer region, L 1 > d SN , a photon emitted with α ≤ π/2 will not be able to reach Earth, so only backward decays are relevant. Here, both signs may result in acceptable solutions provided that sin α ≤ d SN /L 1 . This condition is necessary to guarantee that the photon crosses the sphere with radius d SN at least once. For angles satisfying this condition, the plus (minus) sign indicates the first (second) intersection of the photon with the sphere at r = d SN .
As already mentioned, the probability that a photon is emitted backwards is very small, since this can only happen if the parent ALP is not very boosted (either very heavy or low energetic) and, at the same time, travels beyond r = d SN . This is a highly unlikely scenario and most of the backward decays in the outer region do not reach Earth at allthese photons are therefore essentially lost.
To get an idea of the size of the effects discussed above, let us evaluate the time delay for an ALP with m a = 10 MeV and g aγγ = 10 −10 GeV −1 . Taking E a = 100 MeV as a representative value for its energy, the ALP emits a photon under an angle α ∼ γ −1 ∼ 6 deg. Using d SN = 51.4 kpc for SN 1987A and assuming that the ALP decays after L 1 ∼ ALP ∼ 0.13 pc, Eqs. (17) and (20) show that the overall time delay would be 6 (cf. Fig. 5) This example shows that the time delays may be significant, potentially spreading the signal over a period that is much longer than the duration of the SN explosion (∼ 10 s).
Repeating this exercise for points in the allowed region in parameter space shown in Fig. 7 we would get even larger effects. For m a = 1 MeV and g aγγ = 10 −12 GeV −1 , we find that the time delay is ∆t ∼ 3 × 10 8 s, whereas the angular halo is ∆φ ∼ 1 deg (cf. Fig. 6).
So far we have assumed a fixed emission angle of the photon with respect to the original ALP direction. Let us now justify this assumption. Since the ALP is a pseudo-scalar, in its rest frame, photon decay (emission) is equally likely in any direction, i.e., it is isotropic.
Taking the Lorentz boost that brings the ALP from its rest frame into ours, where it travels with finite speed β, the originally isotropic angular distribution is distorted and is translated into an anisotropic one. To see this explicitly we consider the angular distribution for the separation angle between the two photons ψ = α 1 + α 2 . We find which is peaked at sin(ψ/2) = γ −1 . In the β → 0 limit, the laboratory frame becomes also the ALP rest frame, with the photons emitted back to back, resulting in an increasingly peaked distribution around ψ = π. angular distributions in the context of SN 1987A. We denote as "valid events" the events that pass all the cuts and reach Earth, that is, these are the detected photons.
In order to estimate the size of these effects we first generate the ALP energy distribution (spectrum) for each {m a , g aγγ }-pair based on the massive Primakoff cross section, Eq. (9).
In our numerical simulation roughly 10 7 ALPs are produced per {m a , g aγγ }-pair. For this value the results were stable. The {m a , g aγγ } parameter space itself is scanned in steps of 0.1 (in log scale). This coarse graining produces slightly visible kinks in our limit which are however in line with the level of precision we are aiming for. Using the geometry displayed in ALPs. For masses below a few tens of MeVs it has no significant effect.
After covering the distance L 1 , the ALP decays in two photons. The decay is isotropic in the rest frame of the ALP, but not in the common Earth-SN frame. This can be taken into account by applying an appropriate boost to the ALP and its by-products, whereby we re-obtain the expected focusing of the decays in the forward direction, cf. Eq. (23).
In the ALP's rest frame the ALP-originated photons have equal energies, E a /2 each, but, due to the boost, in our frame their energies are distributed with some spread around this value. Since the detectors have in general a limited energy acceptance, we impose here our second important physical cut by limiting the valid events in the simulation to (boosted) photons with energies in the range [E − , E + ], with E ± determined by the specific detector under consideration. In the case of the original measurements from SN 1987A, the optimal energy range was for gamma rays in the interval 25 − 100 MeV [18,32].
Both simulations, for the time delay and angular distributions, take the aforementioned aspects into account. We highlight again that, in our numerical simulations, the ALP production is taken as being instantaneous, i.e., all ALPs are produced at the same time in the core of the progenitor. We will return to this point in Section IV. Below we present a few representative examples, as well as discuss their most important physical features.

B. Time distribution
In Section II C we discussed the path covered by the ALP and the ensuing photons and we showed that the combined trajectory leads to time delays often longer than the ∼ 10 s duration of the neutrino burst associated with SN 1987A [32], cf. Eq. (22). As mentioned before, backwards photons reaching Earth are very rare, but these will be counted, despite of their relatively small contribution to the overall number of events.
The time-delay simulation follows the steps indicated in Section III A: a number of ALPs is generated with the energy distribution from Fig. 3 and travel a distance L 1 , which is statistically determined by ALP . In the sequence they decay into two photons that cover distances L 2 until detection (cf. Eq. (20) This halo could, in principle, cover a relatively large area of the sky 9 . Considering the unconstrained angular distribution for a given {m a , g aγγ }-pair we obtain distributions similar to the one in Fig. 6. This example shows that it is possible to reach sizeable maximal angular openings in the non-excluded region in parameter space -in this case it is up to ∆φ ∼ 1 • , which is about two times larger than the angular diameter of the Moon.
The fraction of events within the angular acceptance is relevant to determine the sensitivity. Therefore, we are looking for the angular windows ∆φ such that the interval [θ 0 − ∆φ/2, θ 0 + ∆φ/2] centered around θ 0 contains the desired fraction f ang of valid events.
Here we have defined θ = θ 0 = 0 as the direction of line-of-sight to the SN. This is shown for SN 1987A in Fig. 7, where we present the contours of constant angle (for convenience expressed in logarithmic scale as log 10 (∆φ/deg)) in the g aγγ − m a plane for f ang = 90%.
Let us now discuss some features of Fig. 7. combination sin θ ∼ ALP sin α, cf. Eq. (21). It is noteworthy that, in the allowed regions -those with ALP-originated fluence below the experimental upper limit -it is possible to find relatively broad angular windows, reaching ∼ 10 • for large masses and small couplings.
Also interesting is the presence of vertical lines. To understand these, let us consider a fixed detection angle and a roughly fixed ALP energy. For heavy, not very boosted ALPs, the decay angles may vary within a relatively large range, so that many different paths will end up having the same detection angle -even though they have a variety of travelled distances (L 1 ∼ ALP ) and decay angles.
In a sense, the decay length compensates for the freedom in the decay angles. As one goes to smaller masses (i.e., the ALPs are more boosted), the decay angles are quickly more constrained, thus leaving less room for the decay length to compensate. In this low-mass region the detection angle is then dominated by the maximal decay angle ∼ γ −1 , which is independent of the coupling constant (cf. Eq. (4)), hence the saturation and vertical drop-off observed in Fig. 7.
Moreover, in Fig. 7 we have superimposed the boundary of the excluded region where the ALP-originated fluence exceeds the experimental upper limit (cf. Fig. 8). We note that the maximal angular openings within the excluded region are ∼ 0.1 • in the 0.1 − 1 MeV range.

IV. LIMITS FROM SUPERNOVAE
A. SN 1987A The supernova from 1987, whose progenitor Sk -69 202 was a blue supergiant (∼ 18M ) [33], As mentioned in Section III A, we assume that all ALPs are produced in the core of the progenitor at the same time. In practice, according to Ref. [18], this process happens in a time frame of ∼ 10 s. In contrast to the nearly massless case studied in Refs. [18,[25][26][27][28][29][30], for the photons from the decay of massive ALPs the time-delay distributions may be quite broad (cf. Sections II C and III B).
It is therefore advantageous to use a longer time window after the first neutrino recorded.
To do so we look at the full time window of Ref. [32], δt 223 s, and consider the 3 σ statistical fluctuation on the fluence in this period. Since no excess number of events are recorded compared to the control region we use N = 1393 events and the estimate σ = √ N . The upper bound on the fluence for the extended observation time is F exp γ (223 s) ≤ 3 × σ/A eff = 1.78 γ · cm −2 , with the effective area A eff = 63 cm 2 [32]. This effective area takes into account the full field of view of the detector. In the region in parameter space under considerationr the angular distribution is still very narrow -as in the case of photons from ALP conversions in magnetic fields -and we simply follow the same procedure as in Ref. [18,32].
With this upper limit on the fluence we are able to derive the bound presented in Fig. 8 that is also included in the overview Fig. 1. As mentioned in Section II B, the shape of the bound in the low-mass region behaves as g aγγ ∼ m −1/2 a . The non-excluded region increases as the mass decreases due to the decay length: as the masses get smaller, the ALPs are able to survive statistically longer, until the point where they decay predominantly behind Earth, so that, being extremely boosted, very few photons reach us on average, thus suppressing the bound accordingly.
The above behavior is sustained until masses of order m a ∼ 20 MeV and then a turnup takes place, as anticipated in the discussion of the suppression factor, cf. Fig. 4. This is attributed to the size of the physical parameters entering the (massive) Primakoff cross section, Eq. (9), namely, the effective temperature and the Debye screening scale, which are both O(10 MeV). This is the point where the production rate "feels" that the ALPs in the final state are actually heavy, thus consuming a portion of the energy available to convert it into rest mass for the ALPs. This reduces the number of ALPs produced and the fluence is correspondingly suppressed, thus causing the bound to recede.
Another interesting feature is the impact of the effective radius, which is visible on the upper-right region of Fig. 8. This region is characterized by very small decay lengths ALP R * . In this region the ALP-originated photons are absorbed and cannot be detected on Earth, thus explaining the allowed region on the upper-right area (cf. Section II B).
The energy range of the photons used to obtain the bound in Fig. 8 This is more or less the optimal range. Since each ALP decays in two photons, the energy of each photon is distributed around E a /2 and, given that the ALP spectrum achieves a maximum around E a ∼ 80 MeV, most of the ALP-originated gamma rays will be produced with E γ ∼ 40 MeV. In this sense, the optimal energy range for gamma-ray detection should include this value -and possibly even lower ones -in order to cover the range where ALP production is largest.
Let us briefly comment on the uncertainties involved in our limit. There are two uncertainties related to the production process of massive ALPs. One is our simplified modelling of the plasma with an effective temperature and Debye scale. As we have argued in Section II A, we think that this is a relatively small effect. Moreover, there are uncertainties in the modelling of the supernova itself. These have been discussed in Ref. [18] and is assumed to be also reasonably small. In our figures we include only the statistical 3σ uncertainty, but the systematic uncertainties should be kept in mind.

B. Betelgeuse
Betelgeuse is a red supergiant with similar mass as the progenitor from SN 1987A, ∼ 18 M , located in the Orion constellation around 200 pc (650 ly) from Earth. It is expected to explode in a SN event in the next few hundred thousand years [35]. Due to its proximity, it is one of the brightest objects in the night sky and, should it transition into a SN, its associated ALP-originated gamma-ray flux would be much more intense than the one from SN 1987A. Besides this, the gamma-ray instruments have improved in the last decades, so we expect that the overall sensitivity will be significantly better, thus allowing us to set stronger bounds on the ALP parameter space.
Currently, one of the best detection possibilities would be the Fermi-LAT, whose pointsource sensitivity after an observation time of one year (∼ 3 × 10 7 s) is 3 × 10 −9 γ · cm −2 · s −1 for incident photons with energies E γ > 100 MeV [36]. The bounds shown in Fig. 9 behave as expected in the ∼ 1 − 25 MeV-mass region.
However, for masses below ∼ 1 MeV, the g aγγ ∼ m −1/2 a behavior changes to g aγγ ∼ m −1 a . This is due to two intertwined factors: the relatively shorter distance to Betelgeuse and the observation times. (3) and the extra g 2 aγγ from ALP production this gives F γ ∼ g 4 aγγ m 4 a , which leads to From Fig. 9 we see that the change in behavior takes place at different points for different observation times. It is therefore important to estimate where these changes occur. Firstly, we notice that we are limited to a finite observation time δt. Hence, the ALP-originated photons that are effectively counted at the detector need to arrive within this time period, i.e., ∆t ≤ δt and the maximal time delay allowed is δt. Since m a E a we have very boosted s . The fact that m increases with δt is reasonable, since heavier ALPs move slower, so longer observation times are sensitive to larger masses. One must note however that here E a must be such that the photon energies satisfy the constraints of the detector (for E a m a we have E a ∼ 2E γ ).
In fact, this expression roughly matches the transition points for Betelgeuse with the Fermi-LAT detector shown in Fig. 9. For SN 1987A, however,m < 10 keV and the low-mass behavior is not visible 10 in Figs. 8 and 9.
One may wonder why for very small masses an increased observation time does not lead to an improvement in the limit (cf. Fig. 9), but instead actually to a slight weakening. The reason for this is rather simple. In this region the time delay is actually quite small and all photons that will ever arrive already do so in a time smaller than the smallest chosen observation time. Due to the presence of a non-vanishing background for longer observation times we have more background events without any gain in signal events, hence the limit is weaker. . This also causes the displacement of the constant-angle contours in Fig. 10.
For a given {m a , g aγγ }-pair, the angular window for Betelgeuse is expected to be ∼ 200 larger than for SN 1987A, cf. Eq. (21).
Furthermore, the angular acceptance of the Fermi-LAT detector does not strongly constrain P acceptance . However, as mentioned in Section II B, this factor also takes the energy range of the detector into account: for Fermi-LAT we have E γ > 100 MeV. Keeping in mind that E γ ∼ E a /2, from Fig. 3 we see that E a 200 MeV is far from the peak of the ALP production, thus causing the sensitivity to drop by a factor of P acceptance ≈ 0.06. Future experiments like e-ASTROGAM [37], ComPair [40], or PANGU [41] will hopefully be able to improve on this aspect. Yet, even with this reduction, the ALP-originated gamma-ray flux from Betelgeuse is significantly larger due to the closer distance.
As a final remark, we would like to note that a possible ALP burst from Betelgeuse going supernova should not be dangerous to us on Earth. Recent analyses indicate that Betelgeuse would release ∆E ∼ 10 53 erg of energy -similar to SN 1987A -, but the X-and gamma-ray emissions would not be large enough to penetrate Earth's atmosphere [42]. On the other hand, the ALP-related energy release is ∆E ALP ∼ 10 50 erg gaγγ 10 −10 GeV −1 2 , which is much smaller that ∆E. Hence, we do not expect the ensuing ALP-originated gamma-ray burst to pose any harm to us on Earth.

V. CONCLUDING REMARKS
In this paper we have derived new limits on massive ALPs purely coupled to two photons from the supernova explosion SN 1987A. The process we use for our limits is the Primakoff production inside the core of the collapsing progenitor and subsequent decay into two photons. In the 10 keV -100 MeV mass range these limits improve upon existing laboratory and astrophysical limits (see Fig. 1).
Our limits overlap with the cosmological limits discussed in Refs. [14,19], which are based on the effects of the decay of early universe relic ALPs on several CMB and BBN observables. While cosmological limits from thermal production of ALPs in the early Universe are potentially stronger, they are model dependent.
The grey region of Fig. 1 is excluded assuming that ALPs were in thermal equilibrium with the rest of primordial plasma and that the expansion of the universe afterwards was dominated by the relativistic degrees of freedom of the SM. These assumptions set the relic abundance of ALPs and require a sufficiently large maximum temperature of the Universe, where g * and g q are the energy and electric charge effective number of relativistic species, respectively (see Ref. [14]). Of course, achieving such a large temperature depends on the cosmological model considered. Many models feature lower reheating temperatures, which would not be large enough to produce the thermal abundance assumed in the constraints.
FIG . 11. The contour plot shows the average flux (in units of γ · cm −2 · s −1 ) of ALP-originated photons from SN 1987A between 30 and 40 years after its explosion with detection energy E γ ≥ 5 MeV. We see that for a detector with the same point-source sensitivity as Fermi-LAT (3×10 −9 γ · cm −2 ·s −1 ), but with a significantly lower energy threshold, a sizeable area of parameter space could be accessed. For comparison we show in blue the limit obtained from the GRS, cf. Fig. 8.
A very conservative lower limit is set by standard BBN, which requires T RH 20 MeV, implying that only for g aγγ > 7×10 −8 GeV −1 the bounds are independent of the cosmological model. The constraints discussed in this paper, even if superficially weaker, do not suffer from this model dependence and imply a more robust exclusion.
As already discussed in Ref. [10], a large fraction of the ALP-originated photons arrives significantly delayed compared to the neutrinos from the ∼ 10 s-long burst. Indeed, depending on the parameter values, we have found that the signal can be spread out over years.
Although this dilutes the signal, it also provides opportunities as photons may be observed today or even in the near future with better instruments than were available in 1987.
In fact, we have also looked for a possible ALP-originated signal from SN 1987A between today (i.e., 30 years after) and ten years from now. Unfortunately, for Fermi-LAT with its energy cut of E γ > 100 MeV we have found that no signal is detectable. This is mainly due to the high energy required for photons to be accepted in that detector. This requirement reduces the signal twofold, as the SN-ALP spectrum is strongly suppressed at energies 200 MeV and the highly energetic ALPs are typically less delayed and arrive earlier (i.e. before today).
If, however, we lower the energy cut to, say, E γ > 5 MeV, it is possible that some ALPoriginated photons reach us in the next ten years. In Fig. 11 we show the average flux of ALPs between 30 and 40 years after the SN 1987A event. With such an improved gamma-ray detector a considerable amount of interesting parameter space could be probed.
Notably, along with the time delay, the signal will also be spread in angles away from the line of sight. While this effect is not very large for SN 1987A, it can become important for future, closer, supernovae which may be observed with gamma-ray instruments with better directional readiness and angular resolution.
As an example we discuss the possibility of Betelgeuse going supernova in the future, leading to a significantly improved sensitivity. We view this as good motivation to investigate potential improvements also in other limits based upon SN 1987A that could be improved by a future SN observation and also consider the experimental "readiness" to maximally exploit such an opportunity.