Updated constraints on axion-like particles from temporal information in supernova SN1987A gamma-ray data

We revise gamma-ray limits on axion-like particles (ALPs) emitted from supernova SN1987A based on Solar Maximum Mission data. We improve and simplify the computation of the expected gamma-ray signal from ALP decays, while also extending it to non-instantaneous ALP emission. For the first time we make use of the temporal information in the data to update the associated ALP-photon coupling limits. For ALP decays, our updated likelihood only mildly affects the limit compared to previous works due to the absorption of gamma rays close to SN1987A. However, for ALP conversions in the Galactic magnetic field, temporal information improves the limit on the ALP-photon coupling by a factor of 1.4.


Introduction
Axion-like particles (ALPs) can arise as (pseudo-)Nambu-Goldstone bosons, associated with the breaking of a global U(1) symmetry (see e.g.refs [1,2] for reviews).They may also appear in string theory compactifications [3,4], residing in the so-called "axiverse" with masses potentially spanning many orders of magnitude [5][6][7].Recently calculated, explicit mass spectra in type IIB string theory confirm this general picture [8,9].In particular, ALPs from string theory need not be as light as their namesake, the QCD axion [10][11][12][13], but can be much heavier.In this work, we are specifically interested in masses up to the GeV scale.
In any case, constraining ALPs across different mass scales is evidently challenging.Matters are also complicated by the fact that ALPs -unlike QCD axions -need not solve the Strong CP problem, thus lacking a well-defined connection to QCD.However, thanks to the realignment mechanism [14][15][16][17][18], ALPs are still excellent dark matter candidates [e.g.19] and may also couple to photons.This enables a large ensemble of experimental searchers to look for them in the laboratory and using astrophysical and cosmological probes (see e.g.ref. [20] for a review of ALP searches).
Particularly important events in astrophysics are supernovae (SNe), such as the corecollapse supernova SN1987A in the Large Magellanic Cloud.As we revisit in this work, SN1987A is unique in allowing us to constrain ALPs emitted during the SN on very different mass scales: On the one hand, it is well known that the extreme path length in astrophysical settings can lead to strong constraints from particle decay [21].The possibility of ALPs decaying into photons has been used in the past to place limits on couplings of "heavy" ALPs in the keV-GeV range [22][23][24][25][26] and, similarly, on neutrinos [22,[27][28][29].
On the other hand, "light" ALPs with masses m a 1 neV may be converted e.g.inside the Galactic magnetic field into photons, which could then be detected.This has been used to place some of the most competitive limits on the ALP-photon coupling in this mass range [30][31][32].
Apart from considering individual SNe, it has also been pointed out that ALP emission from all past SNe gives rise to a diffuse supernova axion background [33], which can be searched for via the conversion or decay of ALP into gamma-ray photons [e.g.34,35].
In this work, we improve the ALP-photon limit by including the available temporal information contained in the Solar Maximum Mission (SMM) gamma-ray data.We first describe the construction of our updated likelihood in section 2 and appendix A. In particular, the theoretical computation of the expected ALP-induced gamma-ray flux from SN1987A is revisited in section 2.3 and appendix B, where we make further progress in the analytical formalism and extend the previous results to non-instantaneous ALP emission.We present limits derived from our updated likelihood in section 3, comparing them to previous works and discussing the differences.Finally, we conclude with an outlook and additional comments in section 4.
The digitised data sets and computational routines used in this work are available on Github at https://github.com/sebhoof/snax.

Constructing the updated likelihood
Construction of the likelihood function (presented in section 2.4) requires us to understand the available data and instrument response (section 2.1), to select a sensible background model (section 2.2), and to compute the expected gamma-ray signals (section 2.3).
Neutrino data.In addition to the telescope observations, three neutrino observatories saw a neutrino burst around the same time: Kamiokande II [41,42], Irvine-Michigan-Brookhaven (IMB) [43,44], and Baksan [45,46].The neutrino data has been analysed in various studies [e.g.47,48], which find that the measurements are consistent with the first neutrino measured in each detector arriving simultaneously (with an uncertainty of less than a second).Since the IMB detector had by far the most accurate clock, the weighted average of arrival times of the first neutrino is essentially identical to the IMB value, which is t ν = (27 341.37 ± 0.05) s after 00:00:00 UTC on 23 February 1987 [44].SMM/GRS data.Around the time of the neutrino burst, the gamma-ray spectrometer (GRS) [49] aboard the SMM satellite was operational.As shown in ref. [27,Fig. 2], and as described in ref. [28], the GRS took around 223 s of data following the arrival of the first neutrino.Afterwards the GRS went into calibration mode for about 10 min before taking data for another 15 min or so.After this second data-taking interval, the detector was switched off while transiting through the South Atlantic radiation anomaly.The authors of ref. [28] decide against using the data from the second data-taking interval due to concerns about the background model.While we think that it would still have been interesting to analyse it, we were unfortunately unable to obtain additional data despite a number of enquiries.Only GRS data associated with solar flares appears to have been designated for long-term storage. 1 This is also unfortunate in light of the slight discrepancies between the digitised data sets (up to 5% shifts), which we discuss together with our digitisation procedure in appendix A.
In figure 1 we show the data for two of the available energy bands, which we digitised from the literature (see appendix A for details).We do not include data from the 4.1-6.4MeV band since it has a negligible effect on our results due to its narrow range.To a lesser extend this is also true for the 10-25 MeV band, as illustrated by the benchmark models in figure 1.
The GRS effective area.The GRS was facing the Sun during the neutrino burst, meaning that gamma rays from SN1987A had to penetrate the walls of the spacecraft in order to reach the detectors.This reduced the effective detector area A eff,j for all energy bands [28].Still, observations of the 847 keV 56 Co line [50] demonstrate that the GRS was technically capable of detecting a gamma-ray burst at MeV energies.
Estimating the effective area is nonetheless one of the major sources of uncertainties in limits derived from the GRS data set.Monte Carlo (MC) simulations have been performed to estimate the effective area of the detector [51, Fig. 1(A)] (see also ref. [52,Fig. 1]).These MC simulations agree within 20-30% [51] with the Earth's gamma-ray albedo flux measurements [53].While this gives an estimate of the uncertainty of A eff,j under "normal" operating conditions, the difference in viewing angle during SN1987A may introduce additional uncertainties (see also the discussion in ref. [32, §4.2.4]).
The authors of ref. [27] quote effective areas of 115 cm 2 and 63 cm 2 for the 10-25 MeV and 25-100 MeV energy bands, respectively.While it is not clear in how far these values reflect the complications of the measurement, it appears that a reduction compared to the nominal effective area has been taken into account (cf.refs [51,52]).
In contrast the authors of ref. [28] quote a larger effective area of 90 cm 2 for the 25-100 MeV band, without providing details about the computation.The authors also argue that, due to the high gamma-ray energies, the effective areas should be close to their nominal values.Following this logic, we would have also expected larger A eff,j values for the other energy bands.
It is unfortunately not possible anymore to validate the calculations of A eff,j due to the lack of information provided.We use the values quoted in the earlier work, ref. [27].This is the more conservative choice, also allowing us to directly compare our results with most of the later literature.

Background model
As discussed in section 2.1, we assume that the arrival time of the first neutrino (t ν ) coincides with the travel time of light to SN1987A.The data can then be divided into an "off" and "on" measurement, where the "off" data is used to fit the background model nuisance parameters.
To analyse the available photon counting data, we use a Poisson likelihood.For the "off" data, we have (up to a constant) where b ij and n ij are the background model prediction and number of photon counts and in the ith time and jth energy bin, respectively, while i ν denotes the index of the time bin that contains t ν .Different background models were analysed in ref. [27] by considering photon count data on the day before and after SN1987A.For a timescale of 36 min before the satellite went into calibration mode, the authors conclude that a quadratic background model should be used to describe the data.However, over shorter timescales -such as the 6 min interval in figure 1 -we find that a linear model is sufficient to describe the data.We parameterise our linear ansatz for fitting the "off" data as where t i is the time at the centre of the ith time bin and ∆t = 2.048 s.Since the coefficients a (0) j and a (1) j only depend on the data in the jth energy bin, we can optimise the partial likelihoods for the jth energy bin independently instead of eq.(2.1).
We quote our best-fitting parameters for the background in table 1, along with the respective effective areas for each energy bin. Figure 1 shows the prediction of this background model for both the "off" and "on" regions of the data (solid red lines).As in all previous works, we too find excellent agreement of all data with the background-only hypothesis.

Signal prediction
Let us now compute the expected number of gamma rays from ALPs emitted during SN1987A.While both conversion (see section 2.3.2) and decay (see section 2.3.3)processes come from Table 1.Overview of parameters for the SMM/GRS properties and data.For each energy bin, we quote the effective area A eff,j and the best-fitting parameters of the background-only model, derived from fitting the "off" data.the same ALP spectrum, the relevant mass scales at allowed couplings are separated by some twelve orders of magnitude.The limits can thus be derived independently with the same likelihood, simply replacing the expressions for the expected number of photons s ij in the ith time and jth energy bin.
In the following we only consider Primakoff production, so let us discuss how including the coalescence process would affect our results.To this end we can rely on coalescence rates computed in an upcoming study [63],2 based on previous SN simulations.The additional contribution to the ALP flux from the coalescence processes has two effects: one is that the coalescence process starts dominating for ALP masses m a 50 MeV, and we will thus obtain stronger limits on g aγ for these masses.The other effect is that the bounds extend to slightly higher masses, namely up to m a ∼ 250 MeV compared to m a ∼ 200 MeV for Primakoff production only.
For the Primakoff-induced ALP flux, we interpolate the normalisation constant C 1 , average energy E * , and exponent α, tabulated in ref. [32,Table 1], using cubic splines.The axion emission spectrum for axion energy E a and emission time t em is, in the limit of m a → 0, given by the parametric form [32,Eq. (2.11)] We then also fit the instantaneous emission spectrum to the parametric form proposed in ref.
[24, Eq. ( 7)] For the reference value of g aγ = 10 −10 GeV −1 , we find Ĉ2 = 2.03 × 10 77 MeV −1 , Teff = 31.3MeV, and κs = 17.3 MeV.While the value for C 2 is slightly lower compared to ref. [24] which is likely due to our different interpolation method for the coefficients in eq. ( 2.3) -the total number of emitted ALPs from fully integrating the spectra only differs by 1-2%.The spectrum for massive ALPs is then approximately obtained by replacing the Primakoff cross section for massless ALPs in eq.(2.4), σ 0 , with the corresponding expression for massive particles, σ(E a ; m a , g aγ , κ s ), as discussed in ref. [24, Eq. ( 9)].

ALP conversion signal
Apart from decays, the g aγ coupling also allows for the mixing of ALPs and photons in the presence of external magnetic fields.Reference [64] was the first to correctly describe the evolution of the ALP-photon system in general magnetic field configurations.
For this work, we consider the Galactic magnetic field (GMF), for which the models of either Jansson & Farrar (J&F) [65] or Pshirkov et al. (P+) [66] are typically used.We choose the J&F model to compare our results with the literature.Note that the g aγ limits from the J&F model are weaker by a factor of 2-3 compared to P+ [32], which makes the J&F model a conservative choice.
There exist several software codes for computing the overall conversion probability from ALPs emitted from SN1987A into photons detected by the SMM satellite.Amongst them are the ALPro [67] or gammaALPs [68] packages, from which we choose the gammaALPs as it includes the J&F model by default.
Since the exact GMF configuration at the time of SN1987A is unknown, the overall ALP-photon conversion probability can only be computed on average by simulating different field configurations.In particular, the gammaALPs code varies the GMF and splits up the line-of-sight between the SMM satellite and SN1987A into a sufficiently large number of "cells" of size L [69].Inside these cells, the local (transverse) magnetic field B is assumed to be constant.Splitting up the path into many cells also naturally implements the matrix formalism described in ref. [64], according to which the ALP conversion probability -to leading order in g aγ -is given by where ω pl is the plasma frequency of the medium inside the cell.If the typical size and strength of magnetic fields are known, MC simulations of the magnetic field can be used to obtain an average conversion rate P aγ from ALPs into photons.
For the benchmark case of m a 10 −11 eV and g aγ = 10 −10 GeV −1 , the authors of ref. [32, §3] find a conversion rate of P aγ = 0.09, which is about a factor of 1.6 larger than the result of gammaALPs.We find a similar discrepancy for the resulting fluence, suggesting that we can successfully replicate the remaining calculations.The authors of ref. [32] state that they closely follow ref.[70] in their computations, but do not provide additional details beyond this.The differences could be due to the further improvements added after the publication of ref. [70], which eventually led to the release of the gammaALPs code.The latter has also been cross-validated by an independent calculation for ref. [71]. 3ince low-mass ALPs are highly relativistic (E γ E a and t t em ) and convert into one photon each, we have where P aγ is computed by gammaALPs, and the ALP emission spectrum is given by eq. ( 2.3).Finally, note that recently connections of eq.(2.6) to Fourier analysis have been explored [72,73], which is also used in the ALPro code.Furthermore, the simple "cell" model may be replaced by Gaussian random fields or full magnetohydrodynamic simulations [74,75].

ALP decay signal
The expected signal from decaying ALPs [23][24][25][26]76] or sterile neutrinos [28,29,77] has been calculated before, with various degrees of analytical and numerical methods such as MC simulations/integration or quadrature.We include a number of improvements compared to previous works, as described in detail in appendix B.Here we only quote the final result, according to which the expected number of photons s ij is given by 8) min = max {t a (t i ) , r env /β} , t (i) max = min t a (t i ) , t a (t geo ) , (2.10) (2.11) (2.12) Note that we assume an instantaneous ALP emission for the decay limits, i.e. ignore the ALP emission time t em considered in appendix B. The non-instantaneous ALP emission becomes more relevant for decays close to SN1987A, which are however already strongly excluded by data or happen within r env .For ALPs that decay further away from SN1987A, the expected signal in the first few minutes becomes relatively flat (see e.g.figure 1), meaning that a non-instantaneous ALP emission only affects the signal in the first few out of the 109 bins in the "on" region of the data.

Updated likelihood
In addition to the nuisance likelihood L off introduced in section 2, we are now in a position to write the complete likelihood as log L(m a , g aγ , a .We show our results with λ G (blue line, density plot), the rescaled limit from ref. [32] (dashed light blue line), and our results with λ P (black line).
Right: Limits from ALP decays at 3σ CL.We compare our results using λ G (blue line) and λ P (dotted light blue line) to those of ref. [24] (dashed light blue line) and the modified MC code from ref. [76] (black line, density plot).We also indicate the approximate scaling regimes of the limits (dashed grey lines; see main text for details).
where s ij is either the signal prediction from ALP conversions, computed in section 2.3.2, or decays, computed in section 2.3.3.Since we are only interested in limits in the (m a , g aγ ) plane, we can "profile out" the nuisance parameters by considering the log-likelihood ratio test statistic λ.For the Poissonian likelihood, this is where â(0) j , â(1) j are estimates to locally maximise L, i.e. given fixed m a and g aγ , while L is an estimate for the global maximum of L.
This in contrast with e.g.ref. [24], where a Gaussian approximation to the likelihood was used, assuming that all observed events in the 223 s measurement window and 25-100 MeV energy bin are background, which leads to where we use σ 2 b = 1393 to match ref. [24].

Results and discussion
Figure 2 compares our results to the literature to validate and explain the differences in our updated limits, presented later in this section.
The left panel of figure 2 compares our ALP conversion limits (black and blue lines) to ref. [32], where we rescaled their limit (dashed light blue line) to account for the difference in conversion rate P aγ , as discussed in section 2.3.2.We set λ ∆χ 2 = 1 since ref.[32] considers a limiting photon fluence of 0.6 cm −2 .While this corresponds to a 3σ upper limit for a 10 s window [27], the confidence level (CL) is only about 1σ for the full 223 s window (see also ref. [24]).
With the rescaling described before, we find excellent agreement with our computations and ref. [32] when using the simple Gaussian likelihood λ G (blue line).Once the timing information is included via λ P , the limit improves by a factor of about 1.4.This is not unexpected since ALP emission from SN1987A mostly happens over a time window of 20 s or so.Neglecting the time dependence essentially treats the signal as equally distributed across the whole time interval under consideration, while the actual standard deviation of the ALP emission time distribution is only about 4 s.In other words: neglecting the temporal information can "dilute" the signal when distributed across many time bins.This is also illustrated by the benchmark model (blue line) in figure 1.
The right panel of figure 2 then shows our results for ALP decay limits (solid and dotted blue lines) compared to the previous results from ref. [24] (dashed light blue lines).For this purpose, we set λ ∆χ 2 = 9.We also simplify the MC routines of ref. [76] along the lines of our derivations in appendix B and by introducing the popular MC integrator Python package vegas [78] as a more efficient integrator compared to brute-force MC simulations.The results are shown as the density plot and black lines in the right panel of figure 2. Despite our improvements, we can see that MC simulations still struggle to correctly capture the low-mass region.This is because, in the low-mass region, the acceptance fraction of the MC simulations are orders of magnitude smaller than for higher masses.The authors of ref. [76] state that the number of direct MC simulations should be O(10 7 ).However, they did not consider the low-mass region, where the number of MC simulations would have be increased according to the decrease in acceptance fraction.
In any case, our updated likelihood λ P for ALP decays does not result in stronger limits despite containing additional temporal information.This not due to e.g. the inclusion of the background nuisance parameters or other effects.In fact, the finer binning in time does not play much of a role since early ALP decay photons are reabsorbed in the envelope near SN1987A.For later decays, the temporal distribution of the arrival photons becomes relatively flat during the first few minutes, thus not containing any useful timing information.This is also illustrated by the benchmark model (dashed black line) in figure 1.
For a better understanding of how the ALP decay limits in the left panel of figure 2 arise, we also indicate their approximate scaling behaviour (dashed grey lines and text) for m a and g aγ , as previously discussed in ref. [24].In the regimes delimited by ) we set E a = E avg , where the average ALP energy is E avg ≈ 102 MeV for ALPs with sub-MeV masses.In the regime of too early decays, d a < r env (g aγ ∝ m −2 a ), we set E a = E 95 , where E 95 is the 95th percentile of the ALP energy distribution.For ALPs with sub-MeV masses, we find that E 95 = 208 MeV.
Finally, in figure 3 we show our updated limits using λ P (black lines) for the more standard 95% CL to allow a direct comparison with other limits in the literature.For ALP conversions (left panel), we also include the simpler λ G treatment to highlight that also the limits at this CL are still a factor of 1.4 stronger when using λ P .
Apart from complementing other limits, also note that ALP conversion after SN1987A excludes part of the parameter space where ALPs can explain the transparency of the Universe realistic and complete than previous approaches, thanks to the inclusion of a background model, one additional energy bin and temporal information in a Poisson likelihood.
The digitised Solar Maximum Mission data for our updated likelihood and the Python/C ++ code for computing the signal prediction are publicly available on Github at https://github.com/sebhoof/snax.
Despite the improved limits from the statistical analysis in this work, one should keep in mind that there are sizeable uncertainties coming from the supernova emission model, where the predicted signals could change by an order of magnitude.We also neglected the uncertainties from the distance to SN1987A (2.3%), a possible systematic shift of the photon counts depending on the choice of digitised data set (5%) and, most importantly, the effective area of the detector (estimated to be at least 20-30%).These effects could result in a relative uncertainty of ∼ 5-7% on the location of the g aγ limit.
With this said, future work lies in improving the ALP emission predictions from supernova modelling.Future nearby supernovae, such as the one predicted for the red supergiant Betelgeuse, will provide significantly more data and a drastically improved sensitivity to the ALP-photon coupling [e.g.24,90].Given that even the limited amount of SN1987A data leads to some of the strongest constraints on ALPs to date, this presents an exciting prospect.In such an event, the more detailed computation, analysis framework, and software code presented here will hopefully prove useful for probing ALP couplings across many orders of magnitude in ALP mass.
Before explaining how we obtain our consensus data of integer photon counts, we note that we found a rather large discrepancy in the 25-100 MeV energy bin, which cannot be explained by inaccuracies in the digitisation procedure.In the time range where we can compare them, the number of photons in O+ is about 40% larger than what we see in Ch+.
Table 2. Comparison of the digitised data sets from ref. [27,Fig. 4] ("Ch+") and ref. [28,Fig. 1] ("O+").We quote the photon counts for our digitised data sets, the relative deviation to the Ch+, and a cross check number of photons derived from information provided in Table 1 in ref. [28].It seems plausible to us that the authors of O+ did not have access to the actual photon count data but rather fluence data, to which they applied their higher value of the effective area of A eff,2 = 90 cm 2 (cf.section 2.1).To rectify this, we multiply the data in the 25-100 MeV energy range of O+ with a factor of 63/90 = 0.7 before proceeding.We can then use the following estimators for the O+ photon counts in each bin: (i) the digitised data point, (ii) the average of the upper and lower error bar, and (iii) the square of half the length of the error bar.Estimators (ii) and (iii) can be used since Fig. 1 of ref. [28] shows symmetrical error bars, suggesting that the authors use the Wald estimate for the uncertainty on n measured photons, i.e. n ± √ n for the 1σ interval.Indeed, at least two of these estimators give the same rounded integer value for all data points.Estimating the photon counts in each bin from Ch+ works in a similar way, except that the authors use asymmetric 2σ error bars, suggesting that they use a more rigorous approximation for their confidence intervals.We find that the commonly used approximation [93] agrees very well with the data.We thus use the following estimators: (i) the digitised data point, (ii) a fit to the upper, and (iii) to the lower error bar, using the approximation above.Again, at least two of these estimators give the same rounded integer value for all data points.

Data Energy band [MeV
How much are the digitised data sets from Ch+ and O+ in agreement?To answer this question, we compare the "on" (223.232s) and the overlapping parts of the "off" data sets (143.36 s) in table 2. We find that the O+ data (after correction for the effective area) is systematically lower by up to about 5%.This might indicate e.g. a difference in the plotting routines used in the publications.We checked that the difference is not due to a wrong calibration in our digitisation routines.As a result, trying to use data from both data sets -to make use of the longer Ch+ "off" data time window and at the same time the more finely binned O+ "on" data -is not conservative.The ∼ 5% elevated background levels in Ch+ would lead to stronger bound from the O+ "on" data; although we find that the g aγ limit from ALP conversion would only be a factor of 1.5 stronger compared to not including temporal information (factor of 1.4 when only using O+ data).We thus need to pick one of the two data sets and, since there is no definitive answer as to which digitised data set is more accurate, we decide to use the O+ data.

B Calculating the photon flux from astrophysical ALP decays
Here we provide a full derivation of the integral in eq.(2.8), making use of the advantages of previous computations [23-26, 28, 29, 76, 77] while using expressions valid for arbitrary decay times t a and a discussion of non-instantaneous ALP emission.For simplicity, we set c = = k B = 1, except when emphasising the difference between times and lengths by reinstating "c" as a factor.

B.1 Geometry and Lorentz boosts
Figure 4 shows the basic geometry of ALP decays after supernova SN1987A.Without loss of generality, we may choose all ALP and photon paths to cross the x -y-plane.In the ALP rest frame, the two decay photons are emitted back to back with energies of m a /2 each.The photon 4-momenta p ± γ,0 in the ALP rest frame transform to the lab frame, where the ALP is moving with speed β in x -direction, via the Lorentz boost Λ: where Since angles are defined via the 3-vector product, we find the emission angle x ± ≡ cos(θ ± ) from the x -component of p ± γ .The photon energy E γ is, in turn, given by the 0-component of p ± γ : where we used that E γ,0 = m a /2 and defined x ± 0 ≡ cos(θ 0 ).One consequence of the relativistic transformations was already pointed out in ref. [24], namely that the decay photons for highly relativistic ALPs are emitted in a narrow forward cone in the lab frame.This essentially means that e.g.geometries with ALP decays behind d, such as the one shown as a slightly transparent path in figure 4, may be neglected.While geometries with ALP decays behind d are possible for non-relativistic ALPs, they do not contribute much to the signal.This is due to their low speed β 1 compared to the short time window that we consider.
Finally note that, due to the relabelling symmetry of the two photons, we may pick either sign eq.(B.3) as long as we include an overall multiplicity factor of two in what follows.We choose to only discuss the "+" sign in eq.(B.3) to simplify the following derivations.While this choice does not correspond to the γ − photon shown in figure 4, it makes it easier to compare to previous results in the literature, e.g.ref. [77, §12.4].

B.2 Instantaneous ALP emission
We wish to obtain the signal prediction s ij in terms of photon counts by integrating the incoming photon flux over the ith time bin and jth energy bin for the effective detector area A eff,j .We also need to consider the spectral distribution of axion energies E a and photon emission angles x 0 (whose distribution is known in the axion rest frame).Similar to previous works, we also assume that photons from axions decaying within the envelope of the SN are fully absorbed (cf.ref. [24]).This leads us to Apart from these experimental parameters, we need to integrate over all (unobserved) variables, viz. the axion emission energies E a , decay angles x 0 in the rest frame, and decay times t a .Since t will be related to the decay time t a , it is necessary to find an expression for t a (t).By applying the law of cosines to figure 4, and using that cos(π − α) = − cos(α), it follows for the path lengths involved that where we defined x ≡ cos(θ) for convenience and all quantities are measured in the lab frame, i.e. the reference frame of the observing spacecraft.We then define the measurement time t in terms of other travel times such that t = 0 coincides with the time measured after the arrival of the first (massless) neutrino, as discussed in the main text.
Replacing the photon path c t γ in eq.(B.5) using eq.(B.6), we obtain a quadratic polynomial in t a .Further rewriting the polynomial with the help of eq.(B.3) and E γ,0 = m a /2, we find that its two solutions are provided that the determinant is non-negative, which can be interpreted as a condition on t: To choose the physical solution for t a in eq.(B.7), we remind the reader that m a > 0 is required for ALPs to decay into two photons.As a consequence, t = 0 is only possible if t a = 0. Any decays with t a > 0 would lead to t > 0 due to the ALPs' subluminal speed β < 1.Since t + a (0) = 0 while t − a (0) = 0, t a (t) ≡ t − a (t) is the physical solution. 5e also note that, in parts of the literature, the linear expansion of t − a has been used, which is [77, §12.4 (B.9) However, we will see that the approximation in eq.(B.9) is not necessary and, in fact, late decays are relevant for parts of the parameter space.Knowing an expression for t a (t) then allows a change of variables t → t a in eq.(B.4).Together with the other unobserved variables, the relevant part becomes Further expanding the integrand of eq.(B.10) using the chain rule yields: where we used that -in our case -the total ALP lifetime equals the lifetime from photon decays, i.e. τ tot,0 = τ aγ,0 = 1/Γ aγ,0 .6Consider now the variable transform x 0 → E γ,0 .By using eq.(B.3), the resulting factor in the integrand combines with the remaining dE γ,0 /dE γ in eq.(B.12) to an overall factor of The transformation of the x 0 integral boundaries can be understood by writing them as Θ(x 0 + 1) Θ(1 − x 0 ) = Θ(1 − x 2 0 ).Using eq.(B.3) and E γ,0 = m a /2, one finds that We note that t a (t i ) < t a (t i ) due to t i < t i , 7 while r env /β < t a (t geo ) as long as r env < d/2.We can also derive conditions on E a by comparing the other two remaining combinations of possible t a limits, which also improves the numerical convergence of the integral.In practice, the easier condition on E a comes from t i < t geo in the sense that eq.(B.16) is only non-zero if the following weak condition holds: Another possible condition on E a may follow from demanding that r env /β < t a (t i ).However, this leads to a complicated inequality of a sixth order polynomial in E a , which we did not attempt to simplify further.
Regarding the remaining number of numerical integrals to be computed, eq.(B.16) is as convenient as expressions found in some previous works but without using any approximations.In particular, we do not assume highly relativistic ALPs (β → 1) or the asymptotic result for t a given in eq.(B.9).When some combination of these assumptions is made, or when the t ↔ t a integration is not performed, we recover the formulae previously derived in the literature [25,26,28,29,76,77].

B.3 Non-instantaneous ALP emission
When finite ALP emission times t em are considered, the geometry in figure 4 is left unchanged.As a consequence, eq.(B.5) need not be modified.However, we have to account for the additional time delay in eq.(B.6), which becomes t ≡ t a + t γ + t em − d/c , (B. 19) which gives rise to the condition t a ≥ t em , or Θ(t a − t em ) since ALPs cannot decay before they are emitted.We can then simply replace t → t − t em in all equations of appendix B.2.In particular, the ALP decay time now becomes Overall, the signal computation becomes slightly more involved, as one more integral (over t em ) appears.It is convenient to perform this as the innermost integral, keeping in mind that also one new conditions arises in t (i) min = max {t em , t a (t i − t em ) , r env /β} and t (i)  max = min t a (t i − t em ) , t a (t geo ) .(B.21) In the case of SN1987A, r env /βc ≥ r env /c > t em such that this new condition is trivial. 7Observe that ta = t − a in eq.(B.7) is a product of two terms containing t, t + d and the term in square brackets.Using d > 0, Ea > Eγ, and eq.(B.14), it follows that both these terms are monotonic in t, meaning that ta is monotonic in t.

Figure 4 .
Figure 4. Geometry of ALP (dashed blue line) decays into photons (red lines) after SN1987A.The solid lines are labelled with the variables discussed in the main text, while the slightly transparent lines represent an "extreme," improbable geometry.Credit for the modified SSM satellite picture: G. Nelson/NASA (JSC image library; public domain) and modified SN1987A image: NASA/ESA, P. Challis, R. Kirshner, and B. Sugerman (Hubble image library; CC BY 4.0).
Equation (B.14) can be interpreted as a lower limit of the E a integral since E a ≥ E γ + m 2 a /4E γ (cf.ref.[77,§12.4.5]).This replaces the previous lower limit E a ≥ m a since E γ + m 2 a /4E γ has a global minimum at E γ = m a /2 with value m a .Since A eff,j is an effective constant for the jth energy bin, we can put all ingredients together to find that = max {t a (t i ) , r env /β} , and t (i) max = min t a (t i ) , t a (t geo ) .(B.17) min