First Constraints on the Photon Coupling of Axion-like Particles from Multimessenger Studies of the Neutron Star Merger GW170817

We use multimessenger observations of the neutron star merger event GW170817 to derive new constraints on axion-like particles (ALPs) coupling to photons. ALPs are produced via Primakoff and photon coalescence processes in the merger, escape the remnant and decay back into two photons, giving rise to a photon signal approximately along the line-of-sight to the merger. We analyze the spectral and temporal information of the ALP-induced photon signal, and use the Fermi-LAT observations of GW170817 to derive our new ALP constraints. We also show the improved prospects with future MeV gamma-ray missions, taking the spectral and temporal coverage of Fermi-LAT as an example.

Introduction.-Theextreme astrophysical environments in the vicinity of black holes (BHs), neutron stars (NSs), magnetars, and binary BH and NS mergers have recently emerged as a new tool for probing light darksector physics, complementary to and beyond the traditional arena of stellar and supernova environments [1].Much of this recent progress is driven by data from across the electromagnetic spectrum, as well as neutrinos and gravitational waves (GWs), together with the exciting prospects of multimessenger studies [2,3].In particular, the discovery of the NS merger event GW170817 [4] has opened a new window to beyond-the-Standard Model (BSM) particle searches, such as axions and axion-like particles (ALPs) [5][6][7][8][9][10][11], CP-even scalars [12], and dark photons [13].The purpose of this letter is to utilize the multimessenger studies of GW170817 [4,14,15] to derive new constraints on ALPs.
Our main idea is depicted in Fig. 1.A generic feature of the QCD axion [16][17][18][19][20][21], or any pseudoscalar ALP [22,23] is its coupling to photons g aγγ , which is central to most of the ALP searches [24,25].The relevant Lagrangian is where a is the ALP field, m a is its mass, F µν is the electromagnetic field strength tensor and F µν its dual.ALPs can then be produced via the Primakoff and photon coalescence processes in the NS merger.Depending on their mass and coupling, they could escape the merger environment and subsequently decay into two photons.This ALP-induced photon spectrum, which mostly falls in the (soft) gamma-ray or (hard) X-ray regime, is compared to the flux upper limits obtained by X-ray and gammaray observations of GW170817 for various energy ranges and time intervals relative to the GW signal [15,[26][27][28][29] FIG.
1.An artist's rendition of our main idea.The ALP (dashed line), after being produced in the NS merger, escapes and decays outside the merger environment into photons, which can be detected by the Fermi satellite (or future MeV gamma-ray telescopes).
to derive the first merger constraints in the (m a , g aγγ ) parameter space.
ALP production.-Thehot, dense matter in a NS merger remnant would efficiently produce ALPs coupling to photons.The two main production mechanisms are the Primakoff process γ + p → a + p (p being a proton), where electromagnetic scattering essentially converts a photon into an ALP, and photon coalescence γ + γ → a.These processes occur in medium, where the properties of the photon are modified due to the high densities [1].The photon picks up a longitudinal mode (though we will neglect this degree of freedom in this work), and the two conventional transverse modes are altered -in essence, the photon picks up a mass m γ which is equal to the plasma frequency ω pl .Our Primakoff calculation includes the recoil of the charged particle, as well as the electromagnetic screening of the exchanged photon, while the photon coalescence calculation is computed in the usual manner.In both processes, we include the full distribution functions -Fermi-Dirac or Bose-Einstein.More details pertaining to electromagnetism in a dense medium and the local production rate calculations are presented in the Supplemental Material [30].
The production spectrum is then integrated over the entire merger remnant.We use one of the merger profiles calculated in Ref. [31] (available on Zenodo [32]), which use the ALF2 equation of state [33,34].We include the gravitational redshift and trapping factor that prevents ALPs with sufficiently low kinetic energy from escaping the deep potential well of the merger remnant.In the region of parameter space that we study, the ALP mean free path in the merger remnant is long enough that ALP trapping due to inverse Primakoff or two-photon decay is negligible.
After the two NSs collided in GW170817, the resulting merger remnant [35] survived for about one second before collapsing to a BH [36].While ALPs would be produced in the (cold) constituent NSs before they merge, the production rate strongly increases with temperature [37].When the two NSs merge, their temperature rises from keV [38] to tens of MeV [39].Thus the ALP production rate from the binary NS system as a function of time will dramatically upsurge once the stars merge and heat up.We assume that the merger remnant produces ALPs at a constant rate for one second, starting from the time of collision (heating occurs within a few milliseconds of this time) and ending when the remnant collapses to a BH.Since the merger profiles provided in Ref. [31] are time-independent configurations, we use the same profile for the entire lifetime of the remnant.We have checked that, in the parameter space we eventually constrain, the ALP emission does not substantially cool the remnant.
The time-integrated ALP production spectrum (as observed at infinity) is plotted in Fig. 2 for several different values of the ALP mass.The Primakoff and photon coalescence contributions are displayed separately.The ALP-photon coupling is set to g aγγ = 10 −10 GeV −1 , and gravitational redshift and trapping are included.Without gravitational trapping, the spectra would start at ω a = m a , representing a zero kinetic energy ALP.However, such an ALP could not escape the potential well of the merger remnant, and thus would not reach an observer at infinity.Gravitational trapping is much more severe for high-mass ALPs, for example, the m a = 1 GeV case in this figure.The dominant production mechanism for low-mass ALPs is Primakoff, and for high-mass ALPs is photon coalescence.The switchover occurs somewhere around m a ≈ 100 MeV.For low-mass ALPs, the production peaks around energies of 100 MeV, which is due to the temperature T of the environment in which they are produced.For ALPs with masses m a ≳ 3T , the dominant contribution to their energy comes from their rest mass, so their spectrum peaks at energies just above their rest mass (which is why gravitational trapping so diminishes their spectrum).Photon Spectrum.-Thephoton flux F γ observed at Earth is given by The various parameters in Eq. ( 2) are the photon energy ω γ ; the Earth-merger distance D; the photon arrival time t compared to GW arrival time (with t = 0 being the arrival time of the GW); the ALP-to-photon decay angle z = cos α in the ALP frame (see Fig. S6 of the Supplemental Material); the distance L traveled by the ALP before decaying; the distance L γ traveled by the photon before reaching the detector; the ALP spectrum N a as a function of the ALP energy ω a (cf.Fig. 2); the Jacobian Jac(ω a , ω γ ) of the transformation from ALP energy to photon energy; the ALP velocity β a ; the (boosted) ALP decay length ℓ a ; and the minimal distance R ⋆ that ALPs must reach before decaying such that the emitted photons are unlikely to be absorbed by the ejecta in the outskirts of the remnant.The different contributions to Eq. ( 2) are explained in detail in the Supplemental Material.We take benchmark values of D = 40 Mpc and R ⋆ = 1000 km for a GW170817-like NS merger (see Ref. [13]).An important feature of Eq. ( 2) pertains to the temporal dependence of the photon signal, which is critical for multimessenger studies.Although the production of ALPs from the merger environment lasts for 1 sec and is taken to be time-independent, the photon flux from ALP decay has a nontrivial time dependence, as shown in Fig. 3, coming from the decay geometry (cf.Fig. S6).Two benchmarks (both currently unconstrained by other methods) for the ALP mass and coupling are chosen: Benchmark 1 (m a = 398 MeV, g aγγ = 5.01 × 10 −11 GeV −1 ) and Benchmark 2 (m a = 200 MeV, g aγγ = 1.41 × 10 −12 GeV −1 ), in the left and right panels of Fig. 3, respectively.The top (bottom) panels of Fig. 3 display the flux as a function of the arrival time t (photon energy ω γ ), and the contours of different colors correspond to spectral (temporal) snapshots of photon energy (arrival time).The time window corresponding to Fermi-LAT data [28] is shown by the vertical gray shaded region in the top right panel.
Several observations are in order: (i) For ALPs that are longer-lived, corresponding to smaller values of g aγγ and/or m a (Benchmark 2), the photon flux reaches a temporal plateau around t ∼ 1 sec that persists up to t ∼ 10 4 sec for the photon energies displayed (top right panel in Fig. 3).For this time span t ∼ (1 − 10 4 ) sec, the spectral peak lies at ω γ ∼ 300 MeV corresponding to ω 2 γ d 2 F γ /dω γ dt ∼ 10 −5 MeV/cm 2 /sec (bottom right panel).At these early times, a gamma-ray telescope that has its peak sensitivity around ω γ ∼ O(300) MeV would be most suitable for multimessenger studies.The photon signal typically diminishes around t ∼ O(10 5 ) sec, signifying the latest times that a correlated photon-GW study of the merger is likely to be competitive.At such late times, the flux becomes softer and diminishes, with the spectral peak at ω γ ∼ 50 MeV corresponding to ω 2 γ d 2 F γ /dω γ dt ∼ 10 −10 MeV/cm 2 /sec (see the t = 10 5 sec curve in the bottom right panel); therefore, an instrument with peak sensitivity around ω γ ∼ 50 MeV would be more effective there.(ii) For ALPs with shorter lifetimes corresponding to larger values of g aγγ and/or m a (Benchmark 1), the photon flux reaches a temporal plateau and starts diminishing very early, around t ∼ O(1) sec (see top left panel).For the time span t ∼ (0.1−1) sec, the spectral shape is broad between ω γ ∼ (100−600) MeV, corresponding to MeV/cm 2 /sec (see bottom left panel).At later times t ∼ 10 sec, the spectrum softens considerably, peaking at ω γ ∼ 120 MeV, corresponding to ω 2 γ d 2 F γ /dω γ dt ∼ 10 −5 MeV/cm 2 /sec.From these benchmarks, two lessons for multimessenger studies of ALPs can be learned.Firstly, timing and early coordination with the GW signal (within the first second) are most critical for short-lived ALPs, although the studies generally lose competitiveness for ALPs of all lifetimes beyond 10 5 sec.Interestingly, Fermi-LAT data from GW170817 [28] is just at this threshold.Secondly, instruments that are able to collect photon data early should have peak sensitivity in the hundreds of MeV; conversely, at late times, instruments with softer sensitivity in the tens of MeV are preferred. Constraints.
Fig. 4 displays the main results of our study.The purple shaded region is our 95% CL exclusion derived using Fermi-LAT data [28].Here the merger remnant was assumed to survive for 1 sec, but varying this, and hence the ALP emission duration, between 0.72 and 1.29 sec (the range estimated by Ref. [36]) leads to very little change in the exclusion region.Furthermore, this constraint region depends little on the merger profile employed (see Fig. S3).Benchmark 1 (2) is shown with a red square (purple triangle) in Fig. 4. Some future projections are also displayed.The first projection, labeled as "(stacked)" (dot-dashed purple curve), corresponds to constraints coming from a stacked analysis of several mergers, all taken to be at D = 40 Mpc.The number of mergers in the stacked analysis depends on the predicted merger rate, which is (10−1700)/Gpc 3 /yr [45].Within 40 Mpc, this rate predicts a merger every (2.2−370) years; taking a mission lifetime of ∼ 20 years, we take the optimistic case of 9 mergers in the stacked analysis.The second projection, labeled as "(×10 closer)" (dashed purple curve), is for a single future merger that occurs at a distance ten times closer to Earth (D = 4 Mpc), which is expected to happen every (2200−370,000) years based on the same merger rate.A more likely projection is shown by the short-dashed black curve, obtained for the current Fermi-LAT flux upper limit [28], but with an early observation window of 0.1−100 sec.We also show a more futuristic projection (long-dashed black curve) for a hypothetical instrument with 100 times better flux sensitivity than current Fermi-LAT limit and in a broader energy window of 1 MeV-1 GeV.This can happen in principle, given the whole array of proposed MeV gamma-ray missions, such as AMEGO-X [46], e-ASTROGAM [47], PANGU [48], AdEPT [49], APT [50], GAMMA-400 [51], and GECCO [52].
For comparison, the existing astrophysical constraints are also shown in Fig. 4, including those from horizontal branch (HB) star evolution (shaded red) [53], the diffuse γ-ray background from supernovae (shaded orange) [37], the SN1987A energy loss constraint from neutrino signal (shaded magenta) [37], calorimetric constraints from lowenergy supernovae (shaded green) [37], and the SN1987A (shaded blue) [54,55] and SN2023ixf (shaded gray) [56] analogs of the γ-ray constraints.Also shown are the fireball formation (dashed brown) and exclusion (shaded brown) regions [57,58].While the fireball region is taken directly from Ref. [58], and thus is based on different merger profiles than what we use here, the lower edge of the fireball region is almost independent of the merger profile (see Figs. 3 and S2 in Ref. [58]), and therefore, is unlikely to affect our exclusion region or even the future sensitivity in the small g aγγ region.There are additional cosmological constraints from cosmic microwave background and big bang nucleosynthesis on ALP coupling to photons [59,60], which are model-dependent, and hence, not shown here.
Discussion and Conclusions.-(i)This pertains to the relative utilities of mergers and supernovae in probing ALPs.If one ignores multimessenger and stacked analyses of mergers, supernovae enjoy an advantage that can be understood from the following considerations.Assuming that both merger and supernova emit equal ALP fluxes, and that the ALP decay length L (in the lab frame) is short compared to the distance D between the source and Earth, the photon flux from ALP de- Since D SN, merger ≫ L, the above ratio is dominated by the quantity (D merger /D SN ) 2 , and therefore, a close-by supernova will provide a higher photon flux than a far-away merger.This is the main reason why our GW170817 constraint in Fig. 4 (purple-shaded) is weaker than the SN1987A constraint (blue-shaded).
However, after incorporating temporal information of the photon signal, the situation might change.Mergers offer a natural arena for such temporal multimessenger studies, since the time t = 0 can be defined relatively cleanly as the time of arrival of the GW signal.This is unlike the supernova case, where the timing uncertainty can be really large, especially if the neutrino signal is not detected, as e.g. in the case of SN2023ixf [61].As shown in Fig. 4, a combination of temporal-spectral merger data can be a powerful future probe of hitherto unconstrained parts of ALP parameter space.
(ii) It is important to emphasize that since only flux upper limits from electromagnetic observations are used in our study, the conservative constraints derived here do not depend on the detailed modeling of the complex astrophysical background.However, the question of discriminating ALP-induced photons from astrophysical photons becomes very interesting in the early time window t ∼ (0−1) sec.Generally, photons emanating from BSM species like ALPs are expected to arrive before photons emanating from standard astrophysical processes.This is because BSM species couple feebly to material in the merger environment and escape promptly, thereafter decaying to photons in a more pristine environment away from the merger.As for the astrophysical MeV photon background from r-process nucleosynthesis [62][63][64], an analysis of panchromatic dataset from GW170817 [26] suggests a delayed X-ray/gamma-ray signal with respect to the GW signal, consistent with the Fermi-GBM detection of the electromagnetic signal 1.7 sec after the GW detection [41].Therefore, as shown by our temporal analysis, electromagnetic observations of the merger within the first second of the GW detection (possible with the early-warning system [65]) would be crucial to isolate the ALP-induced signal from astrophysical background.This is the subject of our upcoming work [66].
We also list here other possible follow-up studies.During a NS merger, the nuclear matter profiles change considerably over time [39,67], so ALP production should be calculated on top of these time-dependent profiles (as has been done for supernovae [68]) and ultimately included in the simulations themselves (also done in the case of supernovae for alternative ALP models [69,70]).Finally, the local ALP production rates can be refined by improving the treatment of the photon screening in dense matter [71] and by including other production channels like the electro-Primakoff effect [72,73].This section is devoted to the details of ALP production in NS mergers.We first describe in Section A.1 the properties of photons in dense matter, which must be accounted for in the calculations of ALP production rates in the dense NS merger remnant.The rates for the two ALP production channels -the Primakoff process and photon coalescence -are calculated in Sections A.2 and A.3, respectively.

A.1 Photons in dense matter
In dense matter the behavior of photons is modified, which has consequences for the production of bosons, like ALPs, that couple to photons [1].In medium, the photon maintains two transverse modes, albeit with modified dispersion relations, but also acquires a longitudinal mode which is often called the plasmon (though sometimes the term plasmon is applied to all three polarization states of the photon in medium).We shall neglect this longitudinal mode and consider only the two transverse photon modes.The dispersion relation of the transverse modes was calculated in Ref. [74] (see Ref. [1] for a pedagogical discussion).For simplicity, we take the transverse photon dispersion relation to be (k being the photon momentum) where the photon essentially has a mass equal to the plasma frequency ω pl of the medium.This approximation is common in the literature [55,[75][76][77][78][79][80][81][82].In NS merger or supernovae conditions, the plasma frequency is dominated by the most mobile species of particles possessing electric charge, which in this case are the strongly degenerate and ultrarelativistic electrons.The plasma frequency is given by [74] where α is the fine-structure constant and µ e is the electron chemical potential.No assumptions were made about electron degeneracy in this expression, and the positron contribution is included (but is negligible when µ e ≫ T .)In a plasma, the electromagnetic interaction is screened by the motion of charged particles.Screening in the hot, dense matter in a merger remnant is likely dominated by the nondegenerate protons [1], leading to a Debye screening scale with n p being the number density of protons.Momentum exchanges below this momentum scale are suppressed.This expression can be improved slightly to include the finite degeneracy of the protons [68], but we choose not to do so here, given the other uncertainties in the ALP production.Furthermore, in reality the screening length is set by the longitudinal component of the polarization tensor [1], which in a certain limit yields the Debye screening scale in Eq. (S3).However, the strong interactions in the high-density matter in an NS will alter the polarization tensor [71], and this should be considered in improved calculations of the screening scale in the future.ALPs that couple only to photons can be produced in dense nuclear matter by two mechanisms: the Primakoff process and photon coalescence.We are now ready to calculate the ALP production rates via these two mechanisms.In the next subsection, we derive the local ALP production rate, and in a subsequent subsection, the ALP spectrum at infinity including the gravitational redshift and trapping effect.

A.2 ALP production from the Primakoff process
The Primakoff process, depicted in the left panel of Fig. S1, involves the conversion of a photon to an ALP through electromagnetic scattering.We shall consider only scattering with protons γ + p → a + p and not with electrons γ + e − → a + e − , because the electron population is strongly degenerate, which suppresses that scattering rate.The matrix element for the Primakoff process γ(p 2 ) + p(p 1 ) → a(p 4 ) + p(p 3 ) reads Summing over spins and the two transverse polarizations of the photon pol ε α ε * γ → −g αγ , one ends up with the squared matrix element spin, pol The photon mass m γ will be taken to be the plasma frequency ω pl .In the calculations of the Primakoff process, we are interested in the limit of infinite proton mass.The momentum transfer q = p 1 − p 3 becomes q ≈ (0, q), as the proton rest mass dominates the proton energy and allows little energy transfer.Therefore, p i • q ≈ −p i • q and q 2 ≈ −q 2 .In addition, the limit of infinite proton mass simplifies the dot products These limits reduce the squared matrix element to spin, pol Now we formally take the limit m p → ∞, leaving just three terms in the matrix element.The dispersion relation E 2 2 = p 2 2 + m 2 γ simplifies the squared matrix element further to spin, pol Momentum conservation dictates that q = p 4 − p 2 , so the squared matrix element becomes spin,pol From this matrix element, one can derive the cross section given in Ref. [1].
To get the ALP production spectrum, we now integrate this squared matrix element over the phase space: where f represents the Fermi-Dirac distribution of the proton and g the Bose-Einstein distribution of the photon.The ALPs are assumed to have a sufficiently weak coupling to free-stream through the merger matter, so their Bose-Einstein distribution is not included.The protons are assumed to be non-relativistic, yielding the expression (S10) Next, we use the 3-dimensional δ-function to integrate over p 3 and then set up the coordinate system defined by where u and s are cosines of polar angles and ϕ is an azimuthal angle.Integrating over the angles which are made trivial by our choice of coordinate system gives a factor of 8π 2 .Then we integrate over ϕ, using the energy δ-function, and do the s integral.The procedure here is very similar to that detailed in the appendices in Refs.[83,84].At this point, the expression for the ALP production rate for the Primakoff process is with the condition The p 1 integral is now the easiest to do, and we find that the ALP energy spectrum (we relabel E 4 as E a for clarity) produced by the Primakoff process is given by where we define the functions This is almost our final expression, but we must take into account the screening of the Coulomb interaction discussed above.We do this through the propagator replacement (cf.Eq. ( 6.72) in Ref. [1]): where k S is the momentum scale of the Coulomb screening.In the language of the integration variables in Eq. (S14), taking into account screening corresponds to multiplying the existing integrand by (p

A.3 Photon coalescence
In medium, as well as in vacuum, ALPs can be produced through photon coalescence, as shown in the right panel of Fig. S1.We calculate the rate of ALP production from photon coalescence, taking into account the in-medium photon properties.The squared matrix element for γ + γ → a is The rate of photon coalescence γ + γ → a is given by the phase space integral where g 1 and g 2 are Bose-Einstein distributions for the two incoming photons.Again, we neglect the Bose enhancement factor for the ALPs because we only consider ALPs that couple weakly enough to free-stream through the system.We integrate over p 2 using the 3-momentum conserving δ-function and then set up the coordinate system where u is the cosine of the angle between p 1 and p a .Integrating over the three trivial angles and then converting the ALP momentum integral into an energy integral using the ALP dispersion relation, we obtain After doing the u integral using the energy δ-function, we convert the p 1 integral into an integral over E 1 which can be done analytically.We find the ALP production spectrum from γ + γ → a to be with the constraint m a > 2m γ .

B Axion emission from axisymmetric merger remnants
In this section, the NS merger remnant, in which the ALPs are produced, is described.The spacetime and thermodynamic properties of the remnant are detailed and the integration of the local ALP production rate over the remnant profile is discussed.S3.ALP production spectra for nine merger profiles presented in Ref. [31], for an ALP with ma = 1 MeV, gaγγ = 10 −10 GeV −1 and an ALP emission duration of 1 sec.The bold, orange curve represents the profile we consider in our analysis in the main text of the paper.Gravitational redshift and trapping effect is included.
The differentially rotating NS merger profiles considered in Refs.[31,85] have a spacetime metric where the metric fields η, A, B, ω are functions of r and θ only, as time-independence and axisymmetry are assumed.The coefficient η is known as the lapse, and it is also conventional to introduce the conformal factor ψ where A = B = ψ.The metric determinant is given by √ −g = A 2 Bηr 2 sin θ = ψ 3 ηr 2 sin θ.In Fig. S2, we show the merger profiles of various physical quantities for model (e) presented in the data accompanying Ref. [31].This is the model we use as a benchmark for the calculations in this paper.This particular benchmark profile reaches densities in excess of 2.3n 0 (n 0 = 0.16 fm −3 being the nuclear saturation density) and temperatures in excess of 50 MeV.The density decreases monotonically from core to edge, as can be seen from the top left panel of Fig. S2, while the temperature is peaked at about 5 km from the center (cf. the top right panel), in what is actually a ring-like surface in the interior of the merger remnant.This feature is seen in full numerical relativity simulations of NS mergers [39,[86][87][88].The plasma frequency reaches just above 10 MeV (cf. the middle left panel) and the momentum screening scale is of the order of a couple tens of MeV (cf. the middle right panel), both of which are consistent with the values estimated in supernovae cores [78,82,89,90].As shown in the bottom left panel of Fig. S2, the lapse factor η reaches a minimum value of about 0.5 in the core of the star, indicating that ALPs produced there would need energies of about twice their mass to escape to infinity from the remnant (cf.Eq. (S25)).As shown in the bottom right panel of Fig. S2, the nucleon effective masses decrease as the density increases, and in this case the protons have an effective mass of just over 300 MeV at the center of the merger remnant.Additional data from these merger models can be obtained on Zenodo [32].
To obtain the number of ALPs emitted at infinity per ALP energy per time d 2 N a / dω a dt, we integrate the ALP production rate dΓ a / dω a over the volume of the remnant: The time-integrated ALP production spectra are shown in Fig. S3 for nine merger profiles presented in Ref. [31], for an ALP with m a = 1 MeV, g aγγ = 10 −10 GeV −1 and an ALP emission duration of 1 sec.As we can see from this plot, the different merger profiles yield roughly similar ALP production rate within a factor of few.The bold, orange curve represents the profile we consider in our analysis in the main text.In Eq. (S24) we neglect the possibility that emitted ALPs may not have enough kinetic energy to escape the merger, and therefore would not decay to photons energetic enough to be seen by detectors near Earth.In order to include this "gravitational trapping" [13,37,78] into the ALP production rate, we further restrict the production integrals (S14 and S22) with the condition on the (local) ALP energy [37] where η is the gravitational lapse factor from Eq. (S23). 1 The gravitational trapping prevents those ALPs with low kinetic energy E a − m a from escaping the merger remnant, and thus the ALP production spectrum, counting only ALPs that have enough energy to escape the gravitational potential of the merger, "starts up" only at a value of E a greater than the ALP mass m a .In Fig. S4 we plot the ALP production spectra observed at infinity with (solid) and without (dashed) gravitational redshift and trapping, for the benchmark values of m a = 1 MeV, 10 MeV, 100 MeV, 250 MeV and 1 GeV.The effect of gravitational trapping is to cut off the part of the ALP spectrum that involves ALPs with low kinetic energy.The spectrum 'starts up" only for ALP energies E a ≳ 1.5 − 1.7m a .This figure can be compared to its supernova counterpart, Fig. 11 in Ref. [78], where the spectrum starts up at E a ≳ 1.12m a .The supernova evidently has less significant gravitational redshift and trapping because the dense core of the supernova where most of the ALPs are produced is less compact than a NS merger remnant (cf.Fig. S2).
The gravitational trapping effect becomes more important as the ALP mass m a increases.We find that for m a ≳ 10 MeV, the gravitational trapping cannot be neglected.In contrast, for the supernova case, the trapping effect only amounts to relaxation in the constraint [78] for m a ≳ 200 MeV, and that too, by merely 15% or so.

C ALP decay and geometry
In this section we derive the expected photon flux originating from ALP decays.Although the ALP flux is assumed isotropic in the frame of the NS merger, the photon flux must be obtained from geometry arguments [82], considering the possible ALP decay positions as depicted in Fig. S5.
The photon flux as measured by the detector around Earth is given by Eq. ( 2) of the main text.It is expressed in terms of several parameters including the ALP mass m a , the ALP-photon coupling constant g aγγ and the photon energy ω γ .The remaining quantities are derived from the aforementioned parameters.For example, the energy ω a , 1 In the Schwarzschild geometry present outside a spherical star of mass M and radius R, the lapse factor is η = 1 − 2GM/r for r > R. Inside the spherical star, the lapse factor as a function of r can be obtained once the Tolman-Oppenheimer-Volkoff (TOV) solutions for the enclosed mass M (r) and the pressure P (r) are known [91].In an axisymmetric geometry as in the merger case, η depends also on the coordinate θ. the speed β a and the boosted decay length ℓ a of the ALP are, respectively, given by while the distance traveled by the photon from the ALP decay point to Earth is (cf.Fig. S5) We note that in the triangle geometry of Fig. S5, the triangle is determined by two sides and one angle, but there exist two possible triangles associated to such a configuration.One of the triangles matches the one depicted in Fig. S5 with Eq. (S27) for the length of the third side, while the second possible triangle is discarded, since it has a minus sign in front of the square root in Eq. (S27) and leads to backward detection of the photon.Since D is the Earth-merger distance, we choose to measure time in the differential fluence above with respect to the arrival of the GW signal.On the RHS of Eq. (2) of the main text for the differential fluence, there is the differential ALP number assuming the merger occurs at t = τ and ALPs are emitted during a time period ∆t with a constant profile.In the fluence Eq. ( 2) of the main text, it is obviously retarded in time to take into account the time delay before the resulting photons are observed.The Jacobian in Eq. ( 2) of the main text for the change of variables from ω a to ω γ is subtle due to the existence of two possible values for ω a that can lead to a specific ω γ and z.Indeed, we are left with the following possibilities: for ω γ < m a /2 and z > 0 ⇒ ω + a , for ω γ < m a /2 and z < 0 ⇒ ω − a , for ω γ > m a /2 and z > 0 ⇒ ω ± a , for ω γ > m a /2 and z < 0 ⇒ impossible .gives the probability of ALP decay between L and L + dL.The Heaviside functions in Eq. ( 2) of the main text implement constraints.The first constraint forces the ALPs to decay outside the NS merger environment that would not otherwise allow the emitted photons to be observed.For mergers, we choose R ⋆ = 1000 km.The second constraint is necessary for the triangle to exist, as can be seen from the L γ expression in Eq. (S27).
Finally, the factor of 4πD(L γ + Lz) that appears in the denominator of Eq. ( 2) of the main text deserves a more detailed explanation.(i) Since ALP production is isotropic in the frame of the NS merger, dividing the ALP spectrum -which is the ALP number per unit ALP energy per unit time -by 4πL 2 leads to the ALP spectrum per unit ALP energy per unit time per unit area (ALP solid angle times ALP decay length square).Clearly, integrating over the infinitesimal area element L 2 dΩ θ , with dΩ θ the infinitesimal solid angle element for the angle θ as in Fig. S5, recovers the ALP spectrum.(ii) At the ALP decay point, only photons emitted with the appropriate angle α reach the detector.Considering every decaying ALP as a source thus implies an extra factor of 1/2πL 2 γ .Due to the α-angle distribution, this factor is the properly normalized photon solid angle times photon length square, since integrating over the infinitesimal area element L 2 γ dΩ α -where dΩ α is the infinitesimal solid angle element for the angle α -times the α-angle distribution m 2 a /2ω 2 a (1 − β a z) 2 , gives one.To take into account that for a given pair of ALP emission angle θ and ALP decay length L, only photons with the appropriate angle α reach the detector, a Dirac δ-function must be included: which can be integrated over the infinitesimal area element L 2 dΩ θ for ALPs, generating the appropriate factor 4πD(L γ + Lz) in the denominator, with the remaining integrals to be performed over the pair L and z.
With the time-dependent differential fluence Eq. ( 2) of the main text, it is natural to wonder as well about the spatial extent of the signal.To estimate the angular diameter of the observed photon signal, we assume ALPs emitted at 90 degrees (the y-axis in Fig. S5) decaying to 100 MeV photons after one ALP decay length.The resulting angular diameter is approximately 0.27 arcseconds, which is relatively small, of the order of the asteroid Vesta at its farthest from Earth.
Acknowledgments.-SPH acknowledges discussions with Gustavo Marques-Tavares and Georg Raffelt.BD acknowledges discussions with Jim Buckley and Sebastian Hoof.We thank Catie Newsom-Stewart for Figs. 1 and S5.BD is supported by the U.S. Department of Energy grant No. DE-SC 0017987 and by a URA VSP fellowship.The work of JFF is supported by NSERC.SPH is supported by the U.S. Department of Energy grant DE-FG02-00ER41132 as well as the National Science Foundation grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).KS is supported by the U.S. Department of Energy grant DE-SC0009956.YZ is supported by the National Natural Science Foundation of China under grant No. 12175039, the 2021 Jiangsu Shuangchuang (Mass Innovation and Entrepreneurship) Talent Program No. JSSCBS20210144, and the "Fundamental Research Funds for the Central Universities".

FIG. S2 .
FIG. S2.Merger profiles for the model (e) in the data of Ref. [31]: (a) baryon number density nB (top left panel), (b) temperature T (top right panel), (c) plasma frequency ω pl (middle left panel), (d) screening momentum scale kS (middle right panel), (e) lapse η (bottom left panel), and (f) proton (Dirac) effective mass m * p (bottom right panel).Only matter with a number density greater than 0.1n0 and a temperature greater than 1 MeV is shown.
FIG.S3.ALP production spectra for nine merger profiles presented in Ref.[31], for an ALP with ma = 1 MeV, gaγγ = 10 −10 GeV −1 and an ALP emission duration of 1 sec.The bold, orange curve represents the profile we consider in our analysis in the main text of the paper.Gravitational redshift and trapping effect is included.
FIG. S4.Comparison of ALP emission spectra with (solid) and without (dashed) the gravitational redshift and trapping effect for different ALP mass values.

FIG
FIG. S5.The ALP decay geometry.Detailed descriptions of the various symbols are given in the text.