Getting the most on supernova axions

.

The characterization of the ALP emission from the hot and dense nuclear medium typical of a young SN has revealed to be more complex than naively thought.
In particular, it has been recently realized that for typical conditions of a SN core the pionic Compton processes π + N → N + a [32,33] dominate the ALP production with respect to the nucleon-nucleon (N N ) bremsstrahlung N + N → N + N + a, accounted for in the original calculations [28,31,[36][37][38].
ALPs weakly-coupled to matter, g aN ≲ 10 −8 , are in the free-streaming regime in which, once produced, they leave the star unimpeded, contributing to energyloss and to the shortening of the observable SN neutrino burst [20].The most recent study shows that consistency with observations requires g aN ≲ 8 × 10 −10 for m a ≲ 10 MeV [34].However, this argument cannot be extended to arbitrary large values of g aN .For sufficiently large couplings, g aN ≳ 10 −8 , the SN environment becomes "optically thick" for the produced ALPs so that they cannot free-stream out of the star anymore and become trapped in the SN, analogously to neutrinos [20].This is known as the trapping regime, studied in a number of seminal papers after the SN 1987A event [20,23].More recently, ALPs in this regime were considered in Ref. [31], which treated the free-streaming and trapping limiting cases separately.Despite its relevance, a study connecting properly these two regimes has never been performed.As we show below, a recently formulated "modified luminosity criterion" [39][40][41][42] allows the extension of the energy-loss argument from the freestreaming to the trapping limit, smoothly connecting the two regimes.Essentially, the modified luminosity criterion requires the ALP luminosity L a to be lower than the neutrino luminosity L ν at post-bounce time t pb = 1 s, to avoid an excessive shortening of the neutrino burst.
FIG. 1. ALP luminosity at t pb = 1 s as a function of the ALP-proton coupling gap for three different masses, assuming gan = 0.The black line sets the value of the neutrino luminosity.Dashed lines depict the range of couplings where the ALP luminosity exceeds the neutrino one.
An example of how this method is applied is shown in Fig. 1, where L a is plotted as a function of the ALPproton coupling g ap for three values of the ALP mass and compared with L ν .These results have been obtained by assuming the ALP-neutron coupling g an = 0, in analogy with the Kim-Shifman-Vainshtein-Zakharov (KSVZ) axion model [43,44] which we will adopt as our benchmark in the following.Notably, the ALP luminosity grows quadratically with g ap in the free-streaming regime.After reaching the luminosity peak, ALPs enter the optically thick regime, in which their emission can be approximately described as a blackbody radiation from an axion-sphere [45], in close analogy with the case of neutrinos.In this regime, for increasing values of g ap the radius of the emission surface R a increases, while the local temperature of the axion-sphere T a becomes smaller.The combination of these effects leads the ALP luminosity L a ∼ R 2 a T 4 a to approach a plateau at very high couplings.We highlight that this behaviour has been already observed in other works employing the modified luminosity criterion [46].On the other hand, ALPs with m a ≳ 10 MeV suffer an additional Boltzmann suppression which exponentially reduces the luminosity when T a decreases.
Finally, as originally pointed out in Ref. [25] for sufficiently large couplings, g ap ≳ 10 −7 , the ALP burst might produce an observable signal in underground neutrino Cherenkov detectors.Specifically, ALPs would excite the oxygen nuclei in the detector, whose de-excitation would lead to the emission of observable photons.The nonobservation of this signal in coincidence with SN 1987A allowed the exclusion of an additional range of ALPnucleon couplings.
In this paper, we present an original and detailed study of the physics of trapped ALPs.In particular, in Sec.II we summarize the ALP production channels from nu-clear matter in SNe and discuss how the predicted ALP spectrum changes from the free-steaming to the trapping regime.In Sec.III we provide a smooth extension of the ALP luminosity to the trapping regime, with the goal of generalizing the SN cooling bound to the case of trapped ALPs.Our calculation considers generic massive ALPs, thus broadening the analysis beyond the classic QCD "hadronic axion" case [8,47].In Sec.IV we obtain a bound on strongly-coupled ALPs from the absence of a signal associated with the ALP burst in Kamiokande-II (KII) detector at the time of SN 1987A.This last result is based on a new calculation of the ALP-oxygen interaction cross section, presented in a companion paper [48].Our results allow us to exclude large regions of the ALP parameter space, including the entire region accessible to the next generation of cosmological probes sensitive to axion mass.The robustness of these results is discussed in Sec.V, where we analyze the possible uncertainties associated with these constraints.In Sec.VI we translate constraints on the axion-nucleon coupling on the QCD axion mass for canonical QCD axion models.Finally, in Sec.VII we summarize and conclude.An Appendix with details on the ALP emission calculation and the ALP-absorption mean free path follows.

II. ALP PRODUCTION AND ABSORPTION
ALPs can be produced in a SN core by N N bremmstrahlung N + N → N + N + a [28,31,[36][37][38] and by pionic Compton processes π + N → N + a [28,29,32].Recent N N bremsstrahlung calculations introduce corrections to the naive one-pion-exchange approximation, including many-body effects on the nucleon dispersion relations in the medium and its finite lifetime due to multiple scattering [31].Pionic processes have also been re-evaluated recently in Ref. [34] and now include also contributions from the contact interaction term [49], and from vertices associated to ALP-∆ coupling [50].We refer to Ref. [34] for the state-of-the-art calculation of the emissivities for these processes (see also Ref. [35] for a comprehensive review of the ALP production mechanisms).
For sufficiently high values of g aN couplings, ALPs could be reabsorbed in the nuclear medium inside the dense SN core, by means of reverse processes N + N + a → N + N and N + a → N + π.Starting from the expressions for the ALP production spectrum d 2 n a /dE a dt, introduced in Ref. [34], one can derive the mean free path (MFP) λ a associated to these two processes, as illustrated in Appendix A.
The integrated ALP spectrum over a spherically symmetric SN profile can be calculated as [39][40][41]45] where r is the radial position with respect to the center of the SN core and τ (E * a , r) is the optical depth at is the reduced absorption rate defined as in Ref. [41] and T (r) is the temperature profile.The expression in Eq. ( 2) is evaluated at E * a = E a × α(r)/α( r 2 + s 2 + 2rsµ), allowing us to take into account the gravitational redshift from the point of ALP production to the point of absorption, with α the lapse factor which embodies the gravitational effects.This procedure allows us to take into account ALP emissions in any direction, including backward, at each location.We highlight that Eq. (1) allows us to calculate the ALP spectrum in a unified framework which smoothly connects the free-streaming and trapping regimes.Such a recipe was missing in the previous literature.
We have calculated these quantities numerically, using as reference the 1D spherical symmetric GARCHING group's SN model SFHo-18.8 provided in [51], launched from a stellar progenitor with mass 18.8 M ⊙ [52] and based on the neutrino-hydrodynamics code PROMETHEUS-VERTEX [53].
The code takes into account all neutrino reactions identified as relevant for core-collapse SNe [54][55][56] and includes a 1D treatment of PNS convection via mixing-length theory convective fluxes [57] and muonic physics [56].Our benchmark model, already used in previous analyses (see, e.g., Refs.[41,[58][59][60]), leads to a neutron star (NS) with baryonic mass 1.351 M ⊙ and gravitational mass 1.241 M ⊙ .We highlight that GARCHING group's simulations do not account for the presence of pions in the SN core, since the pion properties in PNS matter are still under debate.Therefore, we estimated the pion chemical potential and a pion abundance by employing the procedure in Ref. [33], including the pion-nucleon interaction as described in Ref. [61].This post-processing addition of pions is justified as long as the impact of pionic matter on the PNS properties is not larger than the impact of muons.For our purpose, we employ a single snapshot of this SN model at t pb = 1 s.Fig. 2 displays the behaviour of the ALP spectrum in the massless case for different values of the ALP-proton coupling g ap at t pb = 1 s.We observe that in the free-streaming regime (g ap = 3 × 10 −10 , solid thick line) the spectrum clearly shows a bimodal shape with two peaks, one at E a ≃ 50 MeV associated to N N bremsstrahlung, and the other one at E a ≃ 200 MeV due to πN process [34].As the coupling grows both production and absorption processes become more efficient.
At any given g ap , inverse pion conversion (when kinematically allowed) is more efficient than bremsstrahlung, resulting in a dramatic reduction of the ALP MFP for E a ≥ m π (see Appendix A).Then the part of the spectrum due to πN process is suppressed because of pionic re-absorption that lowers the second peak of the spectrum for g ap = 3 × 10 −8 (dashed curve) till it washes it out completely for g ap = 10 −6 (dotted curve) and g ap = 3 × 10 −5 (solid thin curve).
On the other hand, since the N N bremsstrahlung is a thermal process, the first peak in the ALP spectrum reflects the temperature of the regions where escaping ALPs are produced.In fact, as the ALP-nucleon coupling increases, ALPs become more and more trapped inside the inner regions of the SN core and the only ones able to escape are those produced in the outer layers, where the temperature is lower.Thus, the peak of the spectrum associated to N N bremsstrahlung is shifted towards lower energies, from E a ≃ 50 MeV for g ap = 3 × 10 −10 in the free-streaming case down to E a ≃ 6 MeV for g ap = 3 × 10 −5 in the strongly-coupled case, where the emission takes place from a region at temperature T ≃ 4 MeV.

III. MODIFIED LUMINOSITY CRITERION
The ALP luminosity can be computed as [40][41][42] where the lower limit of integration m a /α cuts away the fraction of heavy ALPs gravitationally trapped in the interior of the core [34,42,59].Notice that in Eq. ( 4) we have set upper integration limit of the radial coordi-nate at the neutrino-sphere radius R ν in order to consider just the amount of energy carried away by axions which can impact the PNS cooling and the associated neutrino emission duration.An excess in energy-loss during the SN cooling phase would have shortened the duration of the SN 1987A neutrino burst.Therefore, the modified luminosity criterion requires that at t pb ∼ 1 s [38,40] the ALP luminosity L a computed on the unperturbed model must not exceed the total neutrino luminosity L ν provided by the same SN simulation.Namely, for our benchmark model, at t pb = 1 s we have This criterion allows us to exclude the blue region in Fig. 3.Note that in this case we exclude 6 × 10 −10 ≲ g ap ≲ 2.5 × 10 −6 for m a ≲ 10 MeV and 7 × 10 −10 ≲ g ap ≲ 3.5 × 10 −7 for m a ∼ O(100) MeV.
The behaviour of the lower bound corresponding to free-streaming ALPs has already been discussed in Ref. [34], so we will focus on the range g ap ≳ 10 −8 , where ALPs enter the trapping regime.As discussed in the previous Section, strongly-coupled ALPs with m a < m π can escape the SN core only if their energies are lower than the pion mass, where pionic production is not possible.Since the emission is mostly in the range E a ≤ m π the dominant absorption mechanism is inverse bremsstrahlung, which is less efficient than pionic absorption.Thus, high values of the ALP-proton coupling are required to saturate the condition in Eq. ( 5) and the bound settles at g ap ≃ 2.5 × 10 −6 .We observe that the bound updates the results in the trapping regime obtained in the previous literature [31,39].In Ref. [31] the ALP production and absorption took into account N N bremsstrahlung which is the only relevant process in the trapping limit.The origin of the differences with the results obtained in this work are due to different motivations.Most importantly, the opacity is underestimated due to a miscalculation of the total MFP [see Eq. (5.31) therein].In addition, the bound in Ref. [31] is obtained by assuming a different physics case (g an = g ap ) and by employing a less refined criterion based on a rough calculation of the axion-sphere radius.On the other hand, in Ref. [39] ALPs were constrained using the modified luminosity criterion, but employing different emission and absorption rates, not accounting for all the corrections considered in this work and in Ref. [34].On the other hand, for m a ≥ m π pionic processes are possible for all energies.Consequently, the dominant absorption process in the entire range of energies is inverse pion conversion, and the ALP absorption rate is significantly enhanced reducing the ALP luminosity.Then the constraint on g ap excludes a smaller region of the parameter space compared to the case of lower masses.
Strongly-coupled ALPs emitted during a SN explosion may lead to a detectable signal in large water Cherenkov neutrino detectors, as proposed in the seminal paper by Engel et al. [25], in which the authors proposed to look for axion-induced excitation of oxygen nuclei with the subsequent emission of a photon to relax the system A revised calculation of the cross-section for this process, using state-of-the-art nuclear models, is presented in Ref. [48].However, on a dimensional ground we can estimate the order of magnitude of the oxygen absorption cross section as σ ∼ g 2 ap /m 2 N ≃ 10 −46 cm 2 (g ap /10 −9 ) 2 , where m N = 938 MeV is the nucleon mass.This estimate is in remarkably good agreement with the numerical results [48].Then, for the calculation of the ALP-induced signal, we adopted a two-step approach similar to what was done in Refs.[25,63].First, we calculated the population of ALP-induced excited states in 16 O using the continuum random phase approximation approach of Ref. [64].Then, in a second step, we evaluated the de-excitation spectrum of 16 O * through γ-cascades and particle emission.The emission of nucleons and α-particles, as well as the γ-decays of their final nuclides, had also to be included because states above the particle separation energies in 16 O * are initially populated.The γ-and particletransitions were treated similarly to the calculation of the transmission coefficients in the statistical model reaction code SMARAGD [65], as also described in Ref. [66].Summed spectra comprising all emissions were obtained.We highlight that the numerical results for the cross section were derived by considering couplings inspired to the KSVZ axion model, which we adopted as our benchmark model.Finally, among the possible nucleon-nucleon interactions described in Ref. [48], in our estimates we employed the parametrization which gives the most conservative results.Starting from the cross-section σ and the ALP number flux F a , it is possible to calculate the number of events in the detector following the recipe described in Ref. [48].If ALPs were produced in the SN 1987A explosion, for sufficiently high values of g ap , this flux would have been detected at the KII water Cherenkov detector [15,16], as pointed out in Ref. [25].We remark that, our analysis, as well as the original one [25], focuses only on the KII detector, neglecting the SN neutrino detection at Irvine-Michigan-Brookhaven (IMB) observatory [17,18].This choice is due to the different energy thresholds of the two detectors, with KII extending down to E thr ≃ 4 MeV sensitive to all the possible oxygen radiative de-excitation channels [48].On the other hand the IMB threshold E thr ≃ 19 MeV is too high to have a significant axion detection [60] in a region of the parameter space not already ruled out by other experimental constraints [3].
FIG. 3. Summary plot of the bounds in the gap vs ma plane together with the QCD axion band (in yellow).The region labeled SNO is excluded by the search for p + d → 3 He + a (5.5 MeV) solar ALP flux in SNO data [3].The green and blue regions labeled SN1987A are ruled out from the non-observation of extra events inside the KII experiment and by the cooling argument.The orange line with the arrow within the QCD axion band shows the sensitivity of current and future cosmological experiments, ma ≳ 0.15 eV [14,62].See the text for more details.
In particular, the energies of the photons from oxygen de-excitation lie mostly in the range E γ < 15 MeV (see Ref. [48] for further details), so that they result to be out of the reach of IMB detection.In light of the most recent analysis of the SN 1987A neutrino burst, performed in Ref. [60], we have assumed that the first nine events detected in KII during the first second are actually in good agreement with neutrino events expected from SN simulations.Conversely, we do not include the latest three events at t pb ∼ 10 s in our further considerations since they are in tension with state-of-the-art simulations.During 2.7 days around the SN 1987A time the background at KII was n bkg ≃ 0.02 events/ s [15,16].Therefore, one can exclude all values of g ap leading to N ev ≳ 2 n bkg ∆t a in the time window ∆t a of the ALP signal.For ALPs with m a ≲ 20 eV [38], time-of-flight effects are negligible and one can assume that ALPs arrive in the same time window of neutrinos, so that ∆t a = ∆t ν .Using the results in Ref. [48], this leads to the constrain g ap ≳ 3 × 10 −7 .For higher values of the ALP mass, ALPs would travel more slowly than neutrinos and their arrival on the Earth would be delayed with respect to the first neutrino event by t = D m 2 a /2E 2 a [38], where D = 51.4kpc is the SN 1987A distance to the Earth.Furthermore, less energetic ALPs arrive later.Therefore, ∆t a is shifted at later times with respect to ∆t ν and it is also spread due to the energy-dependent delay.Since the oxygen energy levels lie in the range E ∈ [9.55,28] MeV, delays between the least and the most energetic ALPs detectable in KII can be estimated as where E min = 9.55 MeV, E max = 28 MeV and we have neglected the intrinsic dispersion of the ALP signal.Notice that for increasing ALP masses, the emission starts to be Boltzmann suppressed and the time delay rapidly increases.From our analysis we have excluded the green region in Fig. 3.For reference, we also show in red, the region excluded by the search for dissociation of deuterons induced by solar ALPs in the Sudbury Neutrino Observatory (SNO) data [3] (g ap ≳ 4 × 10 −5 for m a ≲ 1 MeV).

V. UNCERTAINTIES ON SN BOUNDS
The SN cooling argument employed in this work to constrain the ALP parameter space may be affected by many uncertainties related to the SN modelling employed in the calculations.Indeed, the SN temperature and density profiles significantly affect the ALP emission rates, since they are very sensitive to the stellar properties.Therefore, this Section is devoted to discuss how robust our results are under variations of the SN inputs.In particular, we will focus on three important aspects: the presence or absence of pions inside the SN core; the dependence of the bounds on the used SN model and the possible ALP absorption by heavy nuclei in the outer SN envelope, potentially relevant for strongly-coupled ALPs.

A. Dependence of the bounds on the SN models
The SN emissivity of ALPs coupled to nuclear matter is very sensitive to the SN temperature and the density profiles, which can be significantly different from one SN model to another.Therefore, it is important to test the robustness of the constraints placed in this work under the change of the SN profile employed.As discussed in Sec.II, the results reported in the previous sections are obtained by using the GARCHING SFHO-18.8 model.In the same way, we have also characterized ALP emission by employing a 1D spherically symmetric and general relativistic hydrodynamics SN model launched from a 18 M ⊙ progenitor [67] and based on the AGILE-BOLTZTRAN code [68,69].The latter code includes a full set of neutrino weak interaction processes [70,71] but, contrary to the GARCHING model, it does not account for convection effects.This SN model predicts a lower temperature and density in the SN core at t pb = 1 s compared to the GARCHING one, affecting the ALP emissivity computed in this work.The left panel in Fig. 4 shows the SN bounds obtained using the GARCHING model (green and blue regions delimited by solid lines) and the AGILE-BOLTZTRAN model (dashed lines).In particular, the lower part of the cooling bound is slightly relaxed by switching to the AGILE-BOLTZTRAN code.This is essentially due to the fact that the GARCHING simulation predicts higher temperatures in the inner regions of the SN core where most of the production of weakly-coupled ALPs takes place.Higher temperatures increase pion abundances in the SN core and the ALP production via pion conversion is significantly enhanced.Furthermore, ALP emission rates via N N bremsstrahlung and pionic processes are both strongly dependent on the temperature [see Eq. ( 4) and Eq. ( 5) of Ref. [34]].However, as well as for ALPs, also the neutrino emission is powered by higher temperatures and densities, so that the neutrino luminosity from the GARCHING simulation at t pb = 1 s results to be 70% higher than in the AGILE-BOLTZTRAN simulation.This effect balances the enhancement of the ALP emissivity in the cooling criterion reported in Eq. ( 5), leading to a difference of about 25%.On the other hand, ALP emission rates in the trapping regime are not significantly affected by the change of the SN model and the change in the upper part of the cooling bound is essentially due to the different neutrino luminosity predicted by the two models.In this regime, the ALP luminosity decreases as g ap becomes larger, as shown in Fig. 1.Thus, the lower neutrino luminosity predicted by the AGILE-BOLTZTRAN code implies that this model excludes slightly larger values of g ap in the trapping regime.Analogously, the KII event bound is not modified significantly, with the AGILE-BOLTZTRAN code excluding a slightly larger region in the small-mass limit.To summarize, the two different SN models induce a ∼ 25 − 30% uncertainty on the SN bounds, as shown by the hatched areas in the left panel of Fig. 4. The above discussion shows that different SN profiles (featuring different temperatures and densities) impact on the neutrino and the axion luminosity in a similar way.This implies that the SN cooling bound obtained by using the neutrino luminosity taken from the simulation (see Sec. III) is not significantly affected.Therefore, we do not expect that even taking SN models predicting very high neutrino luminosities, such as those described in Ref. [72], the SN bounds would change dramatically.

B. Presence of pions in the SN core
The abundance of negatively charged pions in the SN core depends on the Equation of State (EoS) employed in the simulation and a self-consistent treatment of pionic matter in SNe is a topic still under investigation (see Ref. [73] for recent developments).In order to estimate the maximum possible uncertainty coming from the abundance of pions in the SN core, we evaluate the relaxation of the SN bound in absence of pions, to obtain the most conservative estimate of the bound.In this scenario, the only possible production channel is N N bremsstahlung and the cooling bound relaxes.The hatched blue region in the right panel of Fig. 4 shows that the cooling bound in the free-streaming regime would be weakened by a factor ∼ 2, excluding g ap ≳ 1.6 × 10 −9 for m a < 10 MeV, and it would probe masses m a ≲ 150 MeV.In particular, the constrained region of the parameter space results to be smaller at higher ALP masses, since bremsstrahlung is more affected by Boltzmann suppression than pionic processes.On the other hand, as discussed in Sec.II, pion conversions play a marginal role in the emission of ALPs from the SN core when entering the trapping regime.Therefore, the absence of pions in the PNS would have no impact on the SN cooling bound in the trapping regime and on the bound from the non-observation of extra events in the KII experiment.

C. ALP absorption by heavy nuclei
After the onset of the explosion, a large amount of energy is released by neutrinos during the Kelvin-Helmoltz cooling phase (see Ref. [57] for an updated review).Meanwhile, the energy deposited by these neutrinos via capture and scattering events powers a baryonic outflow that expands with supersonic velocities, known as the neutrino-driven wind [74,75].The presence of these nuclei at radii larger than R H ∼ 100 km [76] may represent an additional source of absorption for ALPs escaping the SN.In particular, species with atomic numbers similar to 16 O, such as 9 Be and 12 C, and heavier nuclei, such as 56 Fe, may be excited by ALP resonant absorption similarly to oxygen.This effect may lead to a significant reduction of the ALP flux reaching Earth for sufficiently high couplings.By assuming that the cross section for ALP absorption is of the same order of magnitude of the ALP-oxygen one, evaluated in Ref. [48], it is possible to estimate the reduction of the ALP flux.The ALP absorption rate at each radius can be evaluated as where n H (r) = ρ(r) X H (r)/16 m N is the number density of heavy nuclei obtained by assuming that all the population is composed by species similar to oxygen1 , and ρ(r) and X H (r) are the matter density and the heavy nuclei fraction taken from the SN profile, respectively.Then, the flux of ALPs escaping the SN with a given energy E is reduced by a fraction The green hatched region in Fig. 4 displays the modification of the bound by taking into account this effect.It is remarkable to notice that at g ap = 3 × 10 −5 the ALP flux is reduced just by 20%.Therefore, the KII event bound is significantly affected by absorption effects only in a region of the parameter space which is already excluded by the SNO bound.Comparing the two panels in Fig. 4, one can conclude that the uncertainties induced by the SN models are subdominant compared to the ones related to the presence of pions in the SN core and the possible ALP absorption by heavy nuclei.Thus, the right panel of Fig. 4 is a good estimator of the total uncertainty affecting the SN bounds computed in this work.

VI. SN BOUNDS ON THE QCD AXION MASS
In the case of the canonical QCD axion, the SN bounds placed in this work can be read in terms of constraints on the QCD axion mass.In particular, given the axion coupling to nucleons g aN , the corresponding QCD axion mass is obtained as [77] m a ≃ 6.06 In the case of the KSVZ axion model, which we assumed as benchmark model for our analysis, the SN mass bound can be read directly from Fig. 4, where the solid dark-yellow line and the shadowed yellow band refer to the KSVZ and Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) axion models, respectively [78].One realizes that SN arguments place a bound on the KSVZ axion masses that varies in the range depending on the presence of pions in the SN core.These results, within the uncertainties discussed in the previous sections, are comparable with the bound placed in Ref. [79], which ruled out KSVZ axion masses m a ≳ 16 meV from isolated neutron star cooling.We highlight that in the KSVZ axion model also The red shadowed area shows the region of the parameter space ruled out by the RGB tips bound [87], while the grey and green shadowed region refer to the bound placed by isolated neutron star cooling [79] and the duration of the HB phase in Globular Clusters [81].
the axion-photon coupling g aγ is switched on.However, the SN bounds evaluated in this work are remarkably stronger than constraints from the horizontal branch (HB) stars in Globular Cluster [80,81], excluding g aγ ≳ 0.66 × 10 −10 GeV −1 , which corresponds to m a ≳ 440 meV [77].
Our results can be simply extended also to other canonical QCD axion models by substituting different values for the model dependent constants C aN .In particular, in the DFSZ axion model one has [82,83] C ap = −0.435sin 2 β − 0.182 ± 0.025 with tan β confined to 0.28 < tan β < 140 by perturbative unitarity constraints [78,84] (see Ref. [85] for a more recent evaluation of the DFSZ axion-proton coupling range).The solid blue line in Fig. 5 displays the constraint on the DFSZ axion mass placed in this work as a function of tan β.In particular, SN arguments exclude m a ≳ 5 meV for tan β ≳ 5 and m a ≳ 18 meV for tan β ≲ 1.In the absence of pions this bound is relaxed to m a ≳ 14 meV for tan β ≳ 5 and m a ≳ 28 meV for tan β ≲ 1 (see the dashed blue line).Moreover, since in DFSZ axion models also coupling to leptons are available, it is necessary to take into account also the red giant branch (RGB) tip bound introduced in Refs.[86,87].It is remarkable to notice that, in presence of pions, SN 1987A bounds are comparable to the RGB constraint, as shown by the red line in Fig. 5, representing the bound from Ref. [87].In this work, we have set stringent SN bounds on ALPs coupled to nucleons.In particular, we have extended the bounds on free-streaming ALPs to the case of stronglycoupled ALPs, updating the results in the recent literature, and constrained them from the non-observation of ALP-induced events in KII.The combination of these limits allows us to exclude values of the ALP-proton coupling g ap ≳ 6 × 10 −10 for ALP masses m a ≲ 1 MeV.Therefore, SN bounds strongly constrain the parameter space available for ALPs coupled to nucleons.This is in contrast with the original literature on the SN 1987A bound, which reported the existence of a window around a QCD axion mass m a ∼ O(1) eV which was not excluded by the energy-loss argument nor by the axion events in the KII (see, e.g., [38] for a detailed discussion).This region was dubbed the "hadronic axion window" [8,47].It was later shown that cosmological mass bounds on axions would close this window [9].Nevertheless, in non-standard cosmologies, one can relax the axion mass bounds [88], reopening this region of the parameter space.In this regard, it has been recently shown that in low-reheating scenarios multi-eV axions are allowed by cosmology [89].Therefore, the SN 1987A argument is extremely useful to give an independent constraint: with our new analysis we show that SN bounds are enough to exclude canonical QCD axion models (reported in the yellow band in Fig. 3) for masses m a ≳ O(10) meV.The absence of pions in the SN core, whose abundance is still debated, would relax the aforementioned constraint by a factor ∼ 2. This bound is comparable with the recent one placed from the cooling of young neutron stars [79] and it is stronger than the reach of current and future cosmological experiments, which would probe axion masses m a ≳ 150 meV [14,62], as depicted by the vertical orange line on the QCD axion band in Fig. 3.Note that for m a ∼ O(100) meV, QCD axions in SNe would be in an intermediate regime between the complete freestreaming case and the strongly-trapped one.As shown in Table I for the axion mass range probed by cosmology, the SN axion luminosity at t pb = 1 s would be significantly higher than the neutrino one, implying that the axion emissivity would dominate the SN cooling.This would be in tension even with the early time neutrino signal of SN 1987A.Therefore, contrarily to the neutrino case, it is unlikely that future cosmological probes would find signatures of the QCD axion mass as hot dark matter.Nevertheless, below the SN bound, QCD axions with masses m a ≲ 10 meV may contribute to subleading cosmological dark radiation [14], which is in the reach of future cosmological surveys (see Ref. [90]).• The relative sign between the Compton diagrams and the contact interaction term is different in a + p → n + π + with respect to a + n → p + π − .However, this is not relevant for us since the interference terms involving the contact diagram have been neglected in this work, being of higher order in 1/m N (see Refs. [36,49,50] for further details).
• In processes involving neutral pions only Compton-like diagrams are possible.Moreover, it is necessary to divide by a factor of 2 with respect to the same contribution estimated for charged pions [91].
• The integral on the nucleons kinetic energies in the formulas in Eq. ( 5) of Ref. [34] has to be generalized as where μin is the relativistic chemical potential of the incoming nucleon while μout refers to the outgoing nucleon.

Mean free path
Starting from the expressions for the ALP spectrum of production d 2 n a /dE a dt introduced in Ref. [34], the mean free path (MFP) λ a associated with the inverse pion conversion and inverse bremsstrahlung can be computed as [92] λ where E a and p a are the ALP energy and four-momentum, dΠ a = d 3 p a /(2π) 3 2E a is the ALP phase space and χ = ±1 for pionic processes and bremsstrahlung respectively.This change in sign of the ALP energy is due to the different role that the ALP plays in the two processes.In bremsstrahlung processes, depending on whether the ALP is absorbed or emitted, its energy is released to or soaked up from the nuclear medium and then it must change sign switching from the production to the absorption rate (see, e.g.[38]).On the other hand, in pionic processes the conservation of the energy for non-relativistic nucleons just requires that the ALP and the pion involved must have the same energy and no change in sign is necessary.However, note that in pionic absorption the pion is present in the final state implying the substitution f π → 1 + f π .Let us highlight that the expression of the MFP for inverse bremsstrahlung coincides with Eq. (4.27) in [38] extended to the case of massive ALPs.In particular we have used the expression of the structure functions S σ = Γ σ /(E 2 a + Γ 2 ) s σ in a generic degeneracy regime given in Ref. [31], where the integral functions in s σ have to be evaluated with the substitution x → −x and Γ has to be chosen in order to have properly normalized structure functions [93].The behaviour of the MFP in the massless case m a < 10 MeV at R ≃ 16.5 km is depicted in Fig. A2.In inverse bremsstrahlung the incoming ALP just yields a certain amount of kinetic energy to the system of nucleons.As a consequence, ALPs absorption is favored (shorter MFP) when they have energies of the order of the nucleons mean kinetic energy E kin ∼ 3T ∼ 20 MeV.On the other hand, larger energies cannot be efficiently absorbed by nucleons since their phase space becomes smaller.Therefore, the ALP MFP increases for energies larger than ∼ 20 MeV.On the contrary, the MFP associated to pionic processes decreases monotonically as the ALP energy increases.Since these processes consist in the conversion of an ALP into a pion (nucleons can be considered at rest), all the energy brought by the ALP is available in the center of mass to produce the outgoing pion.Consequently, once exceeded the threshold energy E a = m π , as the ALP energy increases the conversion process becomes more and more probable .

FIG. 4 .
FIG. 4. Summary of the bounds with their uncertainties in the gap vs ma plane together with the QCD axion band.The color code is the same as Fig. 3. Left panel: The hatched regions show the uncertainty on the bounds from the SN 1987A cooling (blue region) and the non-observation of extra events in the KII experiment in coincidence with SN 1987A (green region).In both the cases, solid thin lines delimit the region excluded by the GARCHING model while the dashed thick lines delineate the one excluded using the AGILE-BOLTZTRAN code.Right panel: The blue hatched region showing the uncertainty on the cooling bound is obtained by assuming the presence of pions in the SN core, or not.The green hatched region displays uncertainties regarding ALP absorption by heavy nuclei in the neutrino-driven wind.

FIG. 5 .
FIG.5.Summary plot of the bounds in the tan β v s $, ma plane for the DFSZ axion model.The blue shadowed region, delimited by the blue solid line, displays the region excluded by the SN bounds placed in this work.The hatched blue region between the solid and the dashed blue lines represents uncertainties related to the presence of pions.The red shadowed area shows the region of the parameter space ruled out by the RGB tips bound[87], while the grey and green shadowed region refer to the bound placed by isolated neutron star cooling[79] and the duration of the HB phase in Globular Clusters[81].
FIG. A1.Feynman diagrams of pionic Compton processes for ALP absorption, including also the contact and the delta mediated diagrams.The upper diagrams refer to a + n → π − + p, the middle ones depicts a + p → π + + n, while the lower ones refer to a + n → π 0 + n.

FIG. A2 .
FIG. A2.Behaviour of the ALP mean free path at the value of the upper bound for the case of massless ALPs at R ≃ 20 km.In particular the dashed blue line describes the absorption length by means of inverse bremsstrahlung, the red one by means of inverse pion conversion and the black solid curve shows the total MFP.

TABLE I .
[62]uminosity of QCD axion compared with neutrino luminosity at t pb = 1 s for different axion masses, where ma = 100 meV roughly corresponds to the reach of cosmological searches for hot dark matter axions[62].