Axion signatures from supernova explosions through the nucleon electric-dipole portal

We consider axions coupled to nucleons and photons only through the nucleon electric-dipole moment (EDM) portal. This coupling is a model-independent feature of QCD axions, which solve the strong CP problem, and might arise as well in more general axion-like particle setups. We revise the supernova (SN) axion emission induced by the nucleon EDM coupling and refine accordingly the SN 1987A bound. Furthermore, we calculate the axion flux from a future Galactic SN and show that it might produce a peculiar and potentially detectable gamma-ray signal in a large underground neutrino detector such as the proposed Hyper-Kamiokande. The possibility to detect such a signal offers a way to search for an oscillating nucleon EDM complementary to CASPERe, without relying on the assumption that axions are a sizeable component of the dark matter. Furthermore, if axions from SN produce an observable signal, they could also lead to an amount of cosmological extra-radiation observable in future cosmic surveys.


I. INTRODUCTION
Axions are (pseudo)-scalar fields predicted in many well-motivated extensions of the Standard Model (SM) [1,2]. The most notable example is the QCD axion [3,4], which emerges as an essential ingredient in the Peccei-Quinn (PQ) solution of the strong CP problem [5,6]. More generally, in the context of quantum field theory, axions emerge naturally as the Goldstone bosons of global symmetries that are broken at some high scale f a [3][4][5]7]. Ultra-light axions also appear in other frameworks such as supergravity or string theory [8][9][10]. Besides theoretical motivations, there is a huge attention towards axions since these are excellent candidates to account for some or all of the dark matter that we observe in the Universe [11][12][13].
Low-energy experimental tests depend on the axion effective couplings to photons and matter fields, notably * Electronic address: giuseppe.lucente@ba.infn.it † Electronic address: lmastrototaro@unisa.it ‡ Electronic address: pierluca.carenza@fysik.su.se § Electronic address: luca.diluzio@unipd.it ¶ Electronic address: mgiannotti@barry.edu * * Electronic address: alessandro.mirizzi@ba.infn.it L a = C aγ α 8π a f a F µνF µν + C aΨ ∂ µ a 2f aΨ γ µ γ 5 Ψ + . . . , (1) where F µν (F µν ) denotes the electromagnetic field strength (and its dual), ψ = e, p, n runs over low-energy matter fields, and C aγ,Ψ are naturally expected to be O(1) adimensional coefficients. From a phenomenological perspective, it is often assumed the presence of only one of the previous couplings and studied the possibility to constrain each of them separately. The axion-photon coupling (first one in Eq. (1)) is arguably the most used in experimental searches and phenomenological studies. Notably, in the presence of an external magnetic field, the axion-photon interaction leads to the phenomenon of axion-photon mixing [14]. This effect is exploited by several ongoing and upcoming axion search experiments (see [2,15,16] for recent reviews). The axion-photon coupling would also cause axions to be produced in stellar plasmas via the Primakoff process [17]. Therefore astrophysical observations of the Sun, globular cluster systems and supernovae (SNe) offer unique sensitivity to axion interactions (see [18,19] for reviews). The axion-fermion couplings in Eq. (1) also lead to axion production in different stellar systems, e.g. via electron bremsstrahlung in white dwarfs and red giants [20], or nucleon bremsstrahlung [21][22][23] and pion arXiv:2203.15812v2 [hep-ph] 26 Jun 2022 conversion [24,25] in SNe. Furthermore, experimental techniques sensitive to the fermion couplings have been recently conceived (see e.g. [26][27][28][29]).
Above the scale of QCD confinement, the axion interactions in Eq. (1) stem from an axion effective Lagrangian involving quarks and gluons (as well as other SM fields) where G a µν (G aµν ) denotes the gluon field strength (and its dual). The axion coupling with gluons is the most generic feature in the case of the QCD axion, introduced in order to solve the strong CP problem. The two-gluon coupling would allow for the gluonic Primakoff effect in analogy to what is expected for the photon coupling. This effect might be relevant for thermal axion production in the primordial hot quark-gluon plasma [30][31][32][33][34][35][36][37]. At the same time, the axion-gluon vertex is responsible for an irreducible contribution to the axion couplings to photons and nucleons in Eq. (1). These couplings, however, receive other equally important contributions dependent on the specific UV completion of the model. It is, thus, possible to conceive QCD axion models in which they are suppressed compared to their natural O(1) values, as we argue in Appendix A. In the following, we will assume that the couplings to photons and fermions are suppressed.
Finally, the gluonic vertex induces a model independent nucleon EDM portal interaction 1 (with N = p, n) which leads to an axion-dependent nucleon electric dipole moment (EDM), d N = g d,N a [38]. It should be noted that this interaction is in one-to-one correspondence with the axion-gluon coupling via the relation with C anγ = −C apγ = 0.0033 (15) [39] (i.e., g d,n = −g d,p ≡ g d ), and hence it is a model-independent feature of QCD axions. The nucleon-EDM interaction is particularly important for axion dark matter searches. 1 We observe that the axion-gluon coupling in Eq. (2) is not the only possible "microscopic" source for the nucleon EDM portal. In fact, the latter could arise as well from axion-like particle interactions with the basis of CP-violating SM quark and gluon effective operators, such as for example an interaction of the type aqiγ 5 σµν qF µν , with no relation to the solution of the strong CP problem. Hence, the axion-nucleon EDM portal in Eq. (3) describes a more general class of axion-like particle theories and could be decorrelated from the constraints stemming from the axion-gluon operator.
Indeed, the oscillating axion dark matter field would imprint the same oscillations into the EDM of protons and neutrons [38]. The detection of such a feature is the ambitious goal of the Cosmic Axion Spin Precession Ex-peRiment (CASPERe) experiment [29,40]. Oscillating electric-dipole moments of atoms and molecules can also be generated by the interaction in Eq. (2) [41][42][43].
As we shall see in detail in the discussion below, a strong indirect constraint on the axion interaction with the nucleon-EDM can be derived from the analysis of the SN 1987A observed neutrino signal. This possibility was originally presented in Ref. [38], which provided a backof-the-envelope (but, as it turns out, rather accurate) estimate of the axion emission rate from a SN through the N + γ → N + a process. In absence of a direct axion-nucleon coupling, this rate is proportional to g 2 d . It is well known that an overly efficient axion rate would reduce the duration of the observed SN 1987A neutrino signal [44,45], thus allowing to constraint the efficiency of the above process and, consequently, to set a bound on g d .
Given the relevance of the SN bound to constrain the nucleon-EDM portal, we devote our paper to investigate the bounds and signatures of such interactions from SN observations. The plan of our work is as follows. In Sec. II we present the calculation of the axion emissivity in a SN via the nucleon-EDM portal. In Sec. III, we characterize the bound from SN 1987A. In Sec. IV, we calculate the axion signal from a Galactic SN in a large underground neutrino detector, like Hyper-Kamiokande [46]. For completeness, in Sec. V we present an estimate of the cosmological bound on axion thermalization through the N + γ → N + a process. Finally, in Sec. VI, we discuss the complementary of our findings with other observables related to the nucleon-EDM portal and we conclude. There follow three Appendices. In Appendix A, we discuss model-building aspects of QCD axions with suppressed couplings to photons and matter fields. In Appendix B, we provide further details on the calculation of the SN axion emissivity while in Appendix C we discuss the non-degenerate nucleon limit.

II. SUPERNOVA AXION EMISSIVITY VIA NUCLEON DIPOLE PORTAL
Axions can be produced in SN through the nucleon dipole portal of Eq. (3). The two processes which contribute to the rate are the Compton scattering [38], N + γ → N + a, and the nucleon bremsstrahlung process, N + N → N + N + a. 2 Both processes are shown in contribution while the bremsstrahlung is suppressed by more than one order of magnitude, as further discussed in Appendix C. Therefore, in this section we present only a discussion of the Compton effect. Nevertheless, since the bremsstrahlung contribution has never been considered before in the literature, we provide a detailed derivation of the axion emission rate associated with it in Appendix C.
The matrix element for the Compton process involving a nucleon with a real photon in the initial state (see the left panel in Fig. 1) is where are the final and initial nucleon 4-momenta respectively and g d ≡ g d,n = −g d,p is the coupling to nucleon EDM, with the same magnitude but opposite sign for neutrons n and protons p.
In a SN, the photon acquires an effective mass m γ ≈ 16.3 MeV Y 1/3 e ρ 1/3 14 [48], where ρ 14 = ρ/(10 14 g/cm 3 ) and Y e is the electron fraction. Therefore, we evaluate the spin-averaged squared matrix element from Eq. (5) assuming the photon as a massive boson (with three degrees of freedom), getting where k is the photon 4-momentum and m N is the nucleon mass. The number of axions emitted per unit volume and per unit of time and energy is given by where |M| 2 is given by Eq. (6), and the distribution functions of the different interacting species are the usual Fermi-Dirac or Bose-Einstein distribution, where the + sign applies to fermions, the − is for bosons, and µ i are the chemical potentials for i = p, n, while photons have vanishing chemical potential. Corrections to the dispersion relations E i (p i ) of nucleons are incorporated through the equation [49,50] where U i is the non-relativistic mean-field potential and m * N is the effective nucleon mass in medium (see Ref. [23] for details). For definiteness, we take as benchmark for all the different input necessary to characterize the axion emission the SN model with 18 M progenitor simulated in spherical symmetry with the AGILE-BOLTZTRAN code [51,52].
The differential axion number luminosity, which is defined to be the total number of axions emitted in a specified energy range per unit time from the SN is obtained by integrating Eq. (7) over the SN volume and is given by The energy radiated in axions per unit volume and time, called the axion emissivity, can be calculated directly from Eq. (7) as [53] Phase space integration is performed following the procedure of Refs. [54,55], as documented in Appendix B.
In Fig. 2 we show the axion emissivity as a function of the SN radius r at different post-bounce times t pb (bottom right panel), together with the physical properties which determine it, namely the temperature T (upper left panel), the matter density ρ ∼ O(10 14 ) g cm −3 (upper right panel) and the effective photon mass m γ ∼ 15 MeV in the core (bottom left panel). In particular, at t pb = 1 s, the production zone is at r ∼ 10 km and it moves towards the star center at larger times, reflecting the behaviour of the peak temperature T max ∼ 30 − 40 MeV and showing the strong temperature dependence of the production rate. The axion energy luminosity, i.e. the energy emitted by axions per unit time, is obtained by integrating the emissivity over the stellar volume, i.e. L a = 4π dr r 2 Q a (r) .
We mention that redshift corrections need to be considered in order to evaluate the luminosity for a distant observer, as discussed in Refs. [56,57]. Indeed, after its emission, an axion will suffer a gravitational redshift before reaching an observer at infinity. This effect is encoded in the "lapse" factor α listed at each radius in the SN simulation data [52]. This means that the observed axion energy at infinity is E obs = E loc × α, where E loc is the axion energy in the local comoving frame of reference, in which SN-simulation data are provided. In addition, for the rate of emission another redshift correction is required, since the proper time lapse of a comoving observer is related to the distant observer time by the lapse function α [52]. Therefore, the contribution from local emission to the luminosity at infinity can be evaluated by including a factor α 2 . Moreover, since all physical properties of the star are given in the comoving reference frame of the emitting medium, a Doppler shift effect ∝ (1 + 2v r ) has to be considered, where v r is the radial velocity of the medium, which is always very small |v r | 1 [56,57]. For this reason, the observed axion luminosity at infinity is given by where Q a is the emission rate evaluated in the local comoving frame of reference. We stress again that since v r 1, the last term in brackets has a small impact on the observed luminosity, while the α factor reduces the luminosity by a factor ∼ 20 − 30%, in agreement with Ref. [56].

III. SN 1987A COOLING BOUND
The observation of the SN 1987A neutrino burst permits to constrain all the exotic energy losses that would significantly shorten its duration. For quantitative estimates, it is normally assumed that the luminosity associated with exotic processes, calculated at a representative time t pb = 1 s after the core-bounce, does not exceed the neutrino luminosity L ν 3 × 10 52 erg s −1 [56,58]. In order to place a bound on the axion coupling g d , we evaluate the axion luminosity adopting the "modified luminosity criterium", (see [56,59,60]) where the integral of the axion emissivity is performed on the emission region with R p = 40 km and the exponential suppression e −τ takes into account the possibility of axion absorption. In particular, e −τ is a directional average of the absorption factor where λ is the axion mean-free path calculated in Eq. (B10) in Appendix B, E a = E a α(r)/α r 2 + s 2 + 2 r s µ is the axion redshifted energy, µ = cos β and β is the angle between the outward radial direction and a given ray of propagation along which ds is integrated.
In Fig. 3, we show the expected bound on g d in the L a vs g d plane. The trend is a typical one often discussed in literature (see, e.g., Ref. [58]). The region for g d 10 −7 GeV −2 corresponds to the free-streaming case, where the axion production is dominated by a volume emission and L a ∝ g 2 d . Conversely, for g d 10 −6 GeV −2 axions enter the trapping regime, where the luminosity is dominated by a surface black-body emission from an "axion-sphere" with radius r a , where L a ∝ r 2 a T (r a ) 4 that is a rapidly decreasing function of r so that L a decreases when g d increases.
We exclude values of g d for which L a 3 × 10 52 erg s −1 , corresponding to the range 6.7 × 10 −9 GeV −2 g d 7.7 × 10 −6 GeV −2 . We notice that the bound on g d in the free-streaming regime is slightly weaker than the simple back-of-theenvelope estimation presented in Ref. [38], namely g d 4×10 −9 GeV −2 . Furthermore, as shown in Sec. V, values of g d larger than what excluded by the trapping limit are excluded by the extra radiation produced by the thermalization of axions in the early Universe. Therefore, in the next section we will focus on couplings below the free-streaming bound.

IV. AXION SIGNAL IN HYPER-KAMIOKANDE
Having calculated the SN axion spectrum produced through the nucleon dipole portal, our goal in this section is to discuss detection possibilities from a Galactic SN explosion with next generation neutrino detectors (see, e.g., Ref. [61] for a review). For definiteness, we focus on the neutrino underground water Cherenkov detector Hyper-Kamiokande, with a proposed fiducial mass of 374 kton [46]. In this case the detection channel is the scattering of the SN axions on free protons in the water producing a visible photon flux. In order to calculate the event rate in Hyper-Kamiokande, one has to consider the SN axion fluence from Eq. (10), including gravitational redshift. This is well-represented by the following quasithermal spectrum (see also [62]) 73 MeV and β = 3.09. This spectrum is shown in Fig. 4 for g d = 6 × 10 −9 GeV −1 and negligible axion mass.
The detection cross section, associated with the process in Eq. (16), is given by where in the last step we used the small axion mass limit and the non-relativistic approximation for nucleons, so The produced photon energy spectrum is given by where d is the SN distance from Earth, and N t is the number of targets in the detector, with N p = 2 the number of free protons per water molecule, N A the Avogadro number and m H2O = 18 g/mol the molar mass of water.
In Fig. 5, we show the event rate in Hyper-Kamiokande for the axion signal via a + p → p + γ (continuous curve) and for a Galactic SN at d = 0.2 kpc, representative of the distance of the red supergiant star Betelgeuse [63]. For comparison, we show the neutrino event rate associated with inverse beta decay process,ν e + p → n + e + (dashed curve) which is the dominant detection channel for SN neutrinos (see, e.g., Ref. [64]). It is interesting to realize that for E 100 MeV, the axion signal emerges over theν e background, offering a potential window of detection. The high-statistics SN neutrino detection can be used as an external trigger for the axion detection. Indeed, it selects a O(10) s time window to look at the coincidence of at least two photons from axions signal. Notably, the accidental background coincidence in a 10 s window is small, less than one event every three years. 3 The number of axion events for E > 100 MeV is given by Fig. 6, we show the number of expected events in Hyper-Kamiokande as a function of g d , for different values of the SN distance. It is apparent that a few hundreds of events would be detected for g d ≈ 6 × 10 −9 GeV −2 near the cooling bound and for a SN explosion at distance d = 0.2 kpc, such as Betelgeuse. Distances up to d 2 kpc would give a handful of events for the same value of the coupling. For a close-by SN at d 0.2 kpc, we expect to observe few events for couplings larger than g d ≈ 2 × 10 −9 GeV −2 . In order to quantify the sensitivity to g d as a function of the SN distance d, in Fig. 7 we show the Poisson probability to detect more than two photon events with E > 100 MeV in Hyper-Kamiokande as a function of the SN distance, for three different values of the coupling, evaluated as We see that for g d = 6 × 10 −9 GeV −2 , there is a nonnegligible probability (P ev 0.5) to detect an axionsignal up to 2.5 kpc. For g d = 4 × 10 −9 GeV −2 , the sensitivity radius is reduced to 1 kpc, and for g d = 2 × 10 −9 GeV −2 to 300 pc. There are ∼ 30 SN candidates in a radius d < 1 kpc [65]. According to our analysis, if one of these goes SN we might expect, together with a huge neutrino signal, a handful of high-energy events associated with the nucleon dipole portal to axions.

V. COSMOLOGICAL BOUNDS ON AXION EXTRA RADIATION
A complementary constraint on the axion nucleon dipole portal can be derived from measurements of extra radiation in the early universe. In fact, below the QCD phase transition, the process N + γ ↔ N + a becomes an effective process to produce a thermal population of axions which would contribute to extra radiation. The axion production rate in the early Universe is given by Γ = n N σ aN →γN where n N is the nucleon thermal number density and the production cross section σ aN →γN is given by Eq. (18). Axions decouple when Γ H, where H is the Universe Hubble expansion rate. Having determined the axion decoupling temperature T D , it is possible to calculate the effective number of relativistic degrees of freedom, appearing as extra radiation, as [2,66] ∆N eff 0.027 106.75 where g * ,s (T D ) are the entropic effective degrees of freedom (normalized to the total number of SM degrees of freedom). The sensitivity of the Planck 2018 data is enough to exclude ∆N eff 0.35 at 95% CL [67], which corresponds to g d 6 × 10 −6 GeV −2 . Therefore, the cosmological bound nicely connects with the exclusion given by the SN 1987A in the trapping regime. We remark that, for m a 1 eV, axions would be too heavy to be considered dark radiation and their constraint from contributing to dark matter is much weaker than the one from ∆N eff .
For values below the SN 1987A bound in the freestreaming regime, g d 7.7 × 10 −6 GeV −2 , axions would decouple before the QCD phase transition. In this case the processes relevant for their thermalization are the ones with gluons, rather than with nucleons. In this case the decoupling temperature can be estimated as T D 4×10 11 (f a /10 12 GeV) 2 [2] (see also [30][31][32][33][34][35][36]), which is typically well above the electroweak scale. Hence, from Eq. (24) it follows that ∆N eff 0.027, which is in the reach of future CMB-S4 surveys [68]. Requiring that the temperature of the Universe was high enough to bring the axion into thermal equilibrium, T RH > T D , CMB-S4 data would be able to probe [33]

VI. DISCUSSION AND CONCLUSIONS
In this work, we provided a careful quantitative investigation of the bounds and signatures of a nucleon dipole portal to axions from a SN explosion. First, we have revised the axion production channels in a SN. The most relevant channels are the Compton and the bremsstrahlung processes, the last of which had never been considered in the previous literature. We present a detailed calculation of the rates associated with both processes in Appendix C. We find that the SN 1987A cooling argument provides the limit g d 6.7×10 −9 GeV −2 . Furthermore, we have shown that for values of g d below this bound and larger than 10 −9 GeV −2 a future Galactic SN explosion within a radius d O(1) kpc would produce a handful of events through the process a + p → p + γ in the Hyper-Kamiokande detector.
In the case of QCD axion, Eq. (4) holds and the bound on g d can be translated into f a 5 × 10 5 GeV. However, we stress that in this case the SN bound on f a due to nucleon-EDM coupling would be weaker than the HB bound due to axion-photon coupling (f a 4×10 6 GeV) [69] and the SN bound due to axion-nucleon couping (f a 4 × 10 8 GeV) [23]. Therefore, the nucleon-EDM axion coupling would be the most important one for axion phenomenology only in the case in which both photon-and nucleon-axion couplings are suppressed. As further discussed in Appendix A, the required cancellation cannot be achieved with only a single tuning, but further non-trivial assumptions are needed.
It is interesting to compare the axion parameter space probed by SNe with sensitivities of other searches in the plane (g d , m a ), as shown in the upper panel of Fig. 8. There, the shaded green area is excluded by the nondetection of an oscillating nuclear dipole moment in experiments looking for a static one (nEDM, see [70]). The dashed lines represent future sensitivity estimates for different phases of the oscillating EDM experiment CASPERe [29]. Our SN 1987A bound is the blue shaded strip between 6.7×10 −9 GeV −2 g d 7.7×10 −6 GeV −2 . Higher values of the coupling are excluded by extra radiation ∆N eff produced after the QCD phase transition. Values 2×10 −9 GeV −2 g d 6.7×10 −9 GeV −2 (dashed curve) can be probed by axion events from a future closeby Galactic SN explosion.
In the lower panel of Fig. 8, we superimpose additional bounds and future sensitivities under the assumption that the origin of g d is the anomalous axion-gluon coupling, depending on f a in the case of the QCD axion (oblique yellow band) [see Eqs. (2)-(4)]. The region below the dotted grey line (g d 2.8 × 10 −22 GeV −2 ) corresponds to f a m Pl , while in the grey region (BBN) axions coupled to QCD are inconsistent with the production of the observed abundance of light elements during Big Bang Nucleosynthesis (BBN) [71]. Note that both CASPERe and the BBN bound rely on the assumption that the axion comprises the whole cold dark matter. Instead, the bounds denoted as "Earth" and "Sun" are due to finite density effects. In fact, in models where the axion mass is down-tuned, 4  mass can be spoiled in high-density stellar environments where the axion field relaxes to a = πf a , implying various experimental constraints (see Ref. [76] for more details). The shaded grey region around m a ∼ 10 −12 eV represents the most conservative Black-Hole superradiance bound [77] (but see also Ref. [78]). Finally, if axions decouple before the QCD phase transition, e.g. via the axion-gluon coupling, their contribution to extra radiation would be in the reach of future CMB-S4 observations [68], improving over existing constraints on g d . In the lower panel of Fig. 8 we show what the reach in g d would be, assuming a reheating temperature of T RH = 10 10 GeV.
We notice that the region probed by CMB-S4 is complementary to the direct search of axion dark matter by the CASPERe experiment [29]. Remarkably, also our SN signal is complementary with CASPERe. Indeed, for masses m a 10 −9 eV, if g d 10 −9 GeV −2 one would not observe any signal in CASPERe, but there is the possibility to get a few axion-induced events from a close-by SN and an excess of extra radiation in CMB-S4. This is a peculiar scenario where the nucleon dipole portal would be invisible to laboratory experiments and would show up only from the sky. happen via a tuning. More recently, exploiting the mechanism proposed in [73], Refs. [74,75] showed that the axion mass can be exponentially suppressed in terms of a Z N symmetry, while solving the strong CP problem and the axion being dark matter. These works motivate the region on the left of the "QCD axion band" in Fig. 8.

Appendix A: QCD axions with suppressed couplings to photons and matter fields
In this Appendix, we explore the question of whether it is possible to conceive a QCD axion model where the nucleon-EDM portal provides the leading axion interaction. This requires in turn that standard axion couplings to photons and matter fields are suppressed with respect to their natural O(1) values.
To formulate the problem in general terms, let us start from the axion effective Lagrangian below the electroweak scale where E/N is the ratio of the QED/QCD anomaly of the PQ current and f = u, d, e, . . . denotes SM Dirac fermions. The ellipses in Eq. (A1) stand for extra terms like off-diagonal fermion currents, including vector ones. This Lagrangian is matched with the axion effective Lagrangian below the scale of chiral symmetry breaking, Λ χ ≈ 1 GeV, which reads where we kept only terms that are relevant for axion phenomenology, namely nucleons (N = p, n), pions, electrons and photons. In fact, the axion-pion coupling is relevant for the axion hot dark matter bound through to axion-pion thermalization channel [79][80][81] while the other couplings are constrained by astrophysical considerations. The Wilson coefficients of the two effective Lagrangians in Eqs. (A1)-(A2) are related as follows (see e.g. [2]) where δ s = 0.038 (5) The condition that we want to impose corresponds to such that axion phenomenology is driven by the nucleon EDM couplings. From an effective field theory point of view the couplings in Eqs. (A4)-(A8) should be regarded as free parameters and hence it is conceivably possible that they are suppressed with respect to the O(1) values suggested by benchmark axion models.
Here, we want to provide a proof of existence of a UV completion that can realize the conditions in Eq. (A9). To this end, we start from the non-universal axion model of Ref. [83], which can realize the nucleo/pion/electrophobic conditions at the price of a single tuning. The model extends the scalar sector of the SM with three Higgs doublets H 1,2,3 ∼ (1, 2, −1/2) and a SM singlet φ ∼ (1, 1, 0). The scalar potential features the non-Hermitian SM invariant operators which imply the conditions (normalizing to the unity the PQ charge of φ, i.e. X φ = 1) where the latter condition arises from the orthogonality between the PQ and hypercharge currents, with H 1,2,3 = v 1,2,3 and v 2 = v 2 1 + v 2 4 + v 2 3 ≈ (174 GeV) 2 the square of the Higgs vacuum expectation value. The Yukawa sector features the following operators, with a non-universal assignment of the PQ charges in the quark sector with a 2+1 structure (i.e. first and second family, denoted by greek indices, are characterized by the same PQ charge) withH 1,2 = (iσ 2 )H * 1,2 . Differently from Ref. [83], that assumed a universal PQ charge assignment in the lepton sector, in order to obtain here a suppressed coupling to photons (see below) we assumē Neglecting flavour mixing, 5 the flavour diagonal axion couplings to the axial current read where we used the value of the QCD anomaly factor N given by and we also used Eqs. (A12)-(A13). Since, by construction, we have then Eq. (A4) implies C ap + C an ≈ 0 (up to O(5%) corrections from δ s ). On the other hand, the condition C ap − C an = 0 requires the tuning (see Eq. (A5)) that is where in the second step we used Eqs. (A12)-(A13). Note that this condition (satisfied for X 3 ≈ 0) also automatically guarantees C aπ = 0 (see Eq. (A6)) and C e ≈ 0 (see Eq. (A7) and Eq. (A20)). Defining the vacuum angles we can express and parametrize the couplings above in terms of the vacuum angles that are subject to perturbative unitarity constraints (see Ref. [83]). Finally, the QED anomaly factor can be split into a quark plus lepton contribution, i.e. E = E Q + E L , which read respectively and hence E/N = 2, corresponding to the photo-phobic coupling C aγ = 0.08(4).
The condition in Eq. (A9) is hence obtained at the prize of a single tuning, i.e. X 3 ≈ 0, where X 3 can be expressed in terms of vacuum angles of the extended Higgs sector, see Eq. (A28). In particular, X 3 = 0 is obtained for cos 2 β 1 = 1/3 which is fully within the perturbative domain of the model (see Ref. [83]). We remark that the advocated level of suppression of axion couplings can be kept also upon including running effects from f a to the QCD scale, albeit within fairly different parameter space regions than in the tree-level case [86] On the other hand, the level of cancellation that can be achieved in the model above (for X 3 ≈ 0) is not yet sufficient to make the nucleon-EDM axion coupling the most important one for axion phenomenology. Indeed, the two strongest astrophysical constraints on f a are due to C aγ and C aN , coming from HB stars [69] and SNe [23], respectively. An order of magnitude suppression in C aγ is sufficient to evade the former, while the latter is a factor ∼ 800 stronger than the SN bound due to nucleon-EDM coupling. Therefore the suppression proposed in the model in Eqs. (A31)-(A36) is enough to evade the HB bound but not sufficient for the SN bound due to C aN . Hence a further order of magnitude cancellation in C aN would be required. This can be achieved at the price of extra tunings. For instance, C aN can be further suppressed by taking into account flavour mixing effects for flavour-diagonal axion couplings (see Refs. [84,85] for details), while C aγ can be modified via an extra KSVZlike fermionic sector along the lines of [87,88] which contributes to the electromagnetic anomaly.
We conclude that although a percent level suppression of axion couplings seems perfectly feasible in explicit QCD axion models, going below that level requires further non-trivial assumptions.
Appendix B: Evaluation of the emissivity and absorption mean free path The axion emissivity due to Compton scattering N + γ → N + a is given by where |M| 2 is given by Eq. (6), where the + sign applies to fermions, the − is for bosons, and µ i are the chemical potentials for i = p, n, while photons have vanishing chemical potential. Following Ref. [54], the integration in Eq. (B1) can be simplified and the emissivity is given by where p a = |p a |, k = |k|, p f = |p f | and with θ the angle between the axion and the photon momenta, and Q, a, b, c given by where the constants κ, γ, are =p a k cos θ . (B6) The integration domain for cos θ is composed by α = sup[−1, cos θ min ] and β = inf[+1, cos θ max ], with α β and Finally, we stress that in the case of interest p a ∼ E a , since we are considering light axions (m a E a ). On the other hand, axions may be absorbed due to the inverse process a + N → γ + N . The absorption mean free path λ is where n N is the nucleon density and σ aN →γN (E a ) is the absorption cross section given by This expression differs from Eq. (B1) due to the absence of the integration over the axion energy, i.e. d 3 p a /(2π) 3 E a , and the interchange between the initial and final states. Therefore, following once more Ref. [54], one obtains

(B10)
Appendix C: Emissivities in the non-degenerate and non-relativistic limit

Compton effect
In order to obtain a simple expression for the Compton emissivity, let us evaluate it in the non-degenerate and non-relativistic limit for nucleons, ignoring also the effective photon mass. The general expression for the emissivity is where k) and P a = (E a , p a ) are the 4-momenta of the initial and final state nucleon N , the photon and the axion, respectively, with E i,f ≈ |p i | 2 /2m N + m N , E k = |k| ≡ k and E a = |p a | ≡ p a . In addition, all the f i 's are considered as Maxwell-Boltzmann distributions, i.e.
where µ i are the chemical potentials for nucleons, while photons and axions have vanishing chemical potential. In the non degenerate limit the Pauli blocking factors are negligible, (1 − f p f ) ≈ 1, and we can approximate To further simplify the expression, let us fix as reference system the center of momentum frame, in which p i = −k and p f = −p a . We can integrate over p i to eliminate δ 4 (P i + K − P f − P a ), obtaining a constraint on k, since and rewrite the transition matrix as being p a · k = |p a | |k| z. Since d 3 k = 4π dk k 2 , d 3 p a = 2π dE a E 2 a dz and d 3 p a = 4π dp f p 2 f , we obtain Integrating over z, and exploiting the relation in the nondegenerate and non-relativistic limit Now, using the relations dE a E 5 a e −Ea/T = 120 T 2 and dp f p 2 f e −p 2 f /2 m N T = π/2 (m N T ) 1.5 , we conclude that We are considering both the contributions of protons and neutron Y N = Y e + Y n ≈ 1, then the emissivity per unit mass ε a = Q a /ρ can be written in the following form For typical SN conditions (ρ = 3 × 10 14 g cm −3 , T = 30 MeV, Y e = 0.3), imposing the cooling bound ε < 10 19 erg g −1 s −1 , the coupling is constrained to be g d 3 × 10 −9 GeV −2 , in rough agreement with the estimation in Ref. [38].

Bremsstrahlung
In this Section, we evaluate the bremsstrahlung emissivity in the non-relativistic and non-degenerate approximation for nucleons, which is thought to be a good approximation in a SN core, and we show that this process is subdominant with respect to the Compton emission for typical SN conditions.
In the bremsstrahlung process, an axion is produced after the interaction between a nucleon (N 1 = N 3 ) and a virtual photon emitted by a proton (N 2 = N 4 = p). We classify the bremsstrahlung as np channel if N 1 = N 3 = n (see Fig. 9) and pp channel if N 1 = N 3 = p (see Fig. 10). Therefore the emissivity is given by the sum of the two channels Q a,B = Q a,np + Q a,pp , where Q a,np is due to np process and Q a,pp is due to the pp channel. Let us start from the np process. In this case the emissivity is given by with where we are averaging the nucleon spins of the initial and final states and the matrix element is (see Fig. 9) being e the electric charge, k = p 2 − p 4 the transferredphoton 4-momentum and g µν the metric tensor. Assuming that nucleons are non-relativistic and non-degenerate, where with i = p, n, and n p = ρY e /m N and n n = ρ(1−Y e )/m N . In addition, due to the non-relativistic approximation p a = (E a , p a ) , |p a | = E a .
Following Refs. [53,[89][90][91] let us introduce the center-ofmomentum variables In the non-relativistic limit, a typical nucleon with kinetic energy E kin has momentum |p i | = √ 2 m N E kin E kin ≈ E a . Then, the three dimensional delta function implies P = P . Due to the non relativistic approximation and using Eq. (C18) , we obtain p 1 · p 2 = m 2 N + 2|p i | 2 , p 3 · p 4 = m 2 N + 2|p f | 2 , Then, given |P| = P , |p i | = p i , |p f | = p f , p i ·p f = p i p f z, the squared matrix element can be written in terms of the new variables as and f 1 f 2 = n n n p 4 2π m N T 3 e −p 2 i /(m N T ) e −P 2 /(m N T ) .
(C22) One can introduce the variables then the delta function in Eq. (C15) can be rewritten as and we integrate over d 3 p 4 using the three-dimensional piece.
We now change the integration variables d 3 p 1 d 3 p 2 d 3 p 3 = 8 d 3 P d 3 p i d 3 p f where 8 is due to the Jacobian, and using where we fix u = v + x due to the δ-function and The pp-channel contribution (see Fig. 10) is given by where being with a the direct-diagram contribution (upper left panel in Fig. 10) and b the exchange-diagram, obtained interchanging the final fermion lines (upper right panel in Fig. 10). In Eq. (B1), S is the symmetry factor where 2 comes from the position where the axion can be attached (upper or lower vertex, see the lower panels of Fig. 10), and 1/4 comes from the identical particles in the initial and final states (for the np process S = 1). Assuming non-relativistic and non-degenerate nucleons, the matrix element in Eq. (C30) can be evaluated following a procedure analogous to the np process and the pp-contribution reads as Q a,pp = ρ 2 Y 2 e 64 π 7/2 which differs from Eq. (C26) due to the replacements (1− Y e )Y e → Y 2 e (since only protons are involved), 32 → 64 in the denominator (due to the symmetry factor S = 1/2) and |M np | 2 v,x,u=v+x → |M pp | 2 v,x,u=v+x , where |M pp | 2 v,x,u=v+x = g 2 d 4πα 2 T ((2v + x) 2 − 4 v z 2 (v + x)) 2 × m N (2v + x) 2 (4 m 2 N (4 v + x) + 7 m N T (4v 2 + 4 v x − x 2 ) + 4 T 2 (8 v 2 x + 4 v 3 + 7 v x 2 + x 3 )) + 16 v z 2 (v + x) (m 2 N (−(4v + x)) + 2 m N T v (v + x)) − 16 T v 2 z 4 (v + x) 2 (9m N + 4 T (3v + x) .