On the interplay between astrophysical and laboratory probes of MeV-scale axion-like particles

Studies of axion-like particles (ALPs) commonly focus on a single type of interaction, for example couplings only to photons. Most ALP models however predict correlations between different couplings, which change the phenomenology in important ways. For example, an MeV-scale ALP coupled to Standard Model gauge bosons at high energies will in general interact with photons, $W^\pm$ and $Z$ bosons as well as mesons and nucleons at low energies. We study the implications of such scenarios and point out that astrophysical constraints, in particular from SN1987A, may be substantially relaxed, opening up new regions of parameter space that may be explored with laboratory experiments such as NA62.

As ALPs naturally emerge from the breaking of a global symmetry at some high scale they can be suitably studied within an effective field theory (EFT) approach [41][42][43]. In such a set-up many different interactions are expected to arise in a correlated way leading to a rich phenomenology. However, most existing ALP studies limit themselves to a single kind of interaction, most prominently the coupling to photons, while all other couplings are assumed to be zero or at least sufficiently small. In this case it is not possible to capture all of the phenomenological aspects of an ALP that has multiple interactions at low energies.
The present work therefore combines and improves various approaches for constraining ALPs from the literature: • For a given ALP interaction at high scales we consider the effect of all induced lowenergy interactions. For example, if the ALP couples to SM SU (2) L gauge bosons we do not only include the resulting W boson interaction but also the effective ALPphoton coupling and the Zγa vertex. In the case of ALP-gluon interactions we include the induced ALP-photon coupling together with the nucleon interactions.
• This approach makes it essential to go beyond the single-operator-study usually employed in the ALP literature (see also Refs. [44,45] for similar recent approaches). We consistently include all relevant effects by taking into account all required operators in both the ALP production and its decay or absorption. This treatment also differs from Refs. [44,46], where it is assumed that the ALPs decay invisibly into some additional light particles.
• We consistently include all relevant interactions to determine the bounds on ALPphoton and ALP-gluon couplings from the cooling of supernova SN1987A. A similar analysis has been performed in Ref. [22], but we substantially improve the calculation of the ALP optical depth following Ref. [23]. In particular, this approach leads to weaker bounds on the ALP-photon coupling than those obtained previously [22,47].
We apply this approach to the case of ALPs with MeV-scale masses that interact with the various SM gauge bosons. We are particularly interested in the question whether such ALPs may be discovered with laboratory experiments. The most promising search channel is the rare decay K + → π + + inv, which will be explored with unprecedented sensitivity at NA62 [48]. In the simplest models, for example in the case of an ALP coupling only to SU (2) L gauge bosons, the parameter region probed by NA62 is already excluded by astrophysical constraints, in particular by the bound from SN1987A. However, as soon as interactions with several different gauge bosons are considered simultaneously, supernova bounds can be strongly suppressed. The reason is that for sufficiently large couplings ALPs are trapped in the supernova core and can no longer contribute significantly to its cooling. We investigate this interplay of astrophysical and laboratory probes in detail and find that NA62 has great potential to explore this so-called trapping regime.
The paper is organized as follows. In section 2 we review the general aspects of the ALP EFT and derive the relevant interactions as well as meson-mixing contributions. We then discuss existing constraints and prospects on MeV-scale ALPs paying particular attention to supernova bounds in section 3. Our results for specific coupling scenarios of ALPs to SM gauge boson are presented in section 4, which is then followed by our conclusions in section 5. More details on the mixing between ALPs and pseudoscalar mesons as well as their contributions to K + → π + + inv. are given in appendices A and B. Finally, we give details on the computation of the ALP luminosity for SN1987A in appendix C and discuss the sensitivity of NA62 to ALPs with mass greater than the pion mass in appendix D.

ALP effective theory
In this work we employ an EFT approach to describe the interactions of ALPs with SM particles at energies below the new-physics scale Λ. For this we assume that the ALP a is a singlet under the SM gauge group, CP odd and that its mass is below the electroweak scale. We additionally assume no new sources of CP violation. Following the notation from Ref. [42], the Lagrangian then reads [41,43] where G a µν , W a µν and B µν describe the field strength tensors of the strong and electroweak sector with C GG , C W W and C BB being the corresponding Wilson coefficients and g s , g and g are the SM gauge couplings. The dual tensors are given byX µν = 1 2 µναβ X αβ ( 0123 = 1). Furthermore, Ψ F denote the SM fermion multiplets and C F are (hermitian) matrices characterizing their couplings strength to the ALP field. In the following we will focus on ALPs interacting predominantly with SM gauge bosons and therefore assume C F ≈ 0 at tree level.
The continuous shift symmetry of the ALP, a → a + const., used to construct the Lagrangian above is broken down to a discrete one by the anomalous G a µνG µν,a interactions and is furthermore softly broken by the explicit ALP mass term m a,0 . Combining the influence of these two terms on the ALP mass, we obtain [3,42] where F 0 ≈ 92 MeV is the pion decay constant. The parameter m a,0 therefore allows us to treat the ALP mass and the ALP-gluon coupling as independent parameters of our theory, in contrast to the case of the QCD axion [1][2][3][4]. The next step is to include the effects of electroweak symmetry breaking (EWSB) into the effective theory for ALPs. This induces several new interactions in the electroweak sector, of which the following ones are most relevant to our discussion [42] Here e is the electric charge, θ w is the weak angle with s w ≡ sin(θ w ) and c w ≡ cos(θ w ), F µν and Z µν are the photon and Z boson field strength tensors, respectively, and W ± are the fields for the W ± bosons.

Mixing effects of light ALPs
As we are interested in studying ALPs with masses below 1 GeV, the hadronic contributions to the processes of interest (discussed in more detail later) can be obtained from a chiral Lagrangian [41]. We start by computing the ALP-mixing with pseudoscalar mesons (π 0 , η, η ). For this we first have to determine the interactions resulting from eq. (2.1) at scales of about 1 GeV [42] where we have included the relevant QCD terms with q = (u, d, s) T and the quark mass matrix M q = diag(m u , m d , m s ). Before mapping this onto a chiral Lagrangian, it is convenient to rotate the gluonic interaction away by performing a chiral rotation on the light quark fields q → e iκq a 2fa γ 5 q where Tr[κ q ] = 1 and 1/f a = −32π 2 C GG /Λ. The interactions then read [41,42] (2.6) Here we introduced the rotated quark mass matrix M q (a) = e iκq a 2fa γ 5 M q e iκq a 2fa γ 5 and the charge matrix Q q = diag(2/3, −1/3, −1/3) containing the electric charges of the light quarks in units of e > 0.
These interactions can now be mapped onto a chiral Lagrangian at leading order (see appendix A for details). Here it is convenient to choose κ q = M −1 q /Tr[M −1 q ] in order to avoid tree-level mass mixing between the ALP and the octet mesons [41]. The mixing contributions of the ALP into the pseudoscalar mesons can then be determined to [49] where = F 0 /Λ and small mixing angles are assumed, i.e. the ALP mass should not be too close to any of the pseudoscalar masses. The kinetic mixing contributions K aP are given by while the mass mixing terms m aP are where B 0 = m 2 π /(m u + m d ) and c(θ) ≡ cos(θ) as well as s(θ) ≡ sin(θ) describe the η − η −mixing. 1

Interactions of light ALPs with photons and nucleons
The coupling studied most in terms of constraints and experimental prospects in the ALP literature is the one to photons, see e.g. Refs. [22,25,47,[51][52][53][54][55]. It is therefore crucial to capture all of the sizeable contributions to the photon coupling in our EFT set-up. Here it is important to notice that the tree-level coupling in eq. (2.4) generated after EWSB receives further hadronic contributions from the gluonic ALP coupling as is evident from eq. (2.6). We therefore define an effective ALP-photon coupling [42,56,57]  where the precise coefficient multiplying C GG can be obtained by including higher orders in the computation [58]. 2 We have also taken the mixing induced contributions into account, which are valid for ALP masses not too close to the respective pseudoscalar mass. We emphasise that the ALP-photon coupling is very generic in the sense that it receives a contribution from all the Wilson coefficients studied in this work and will therefore always be sizeable unless an accidental cancellation occurs. This simple observation demonstrates 1 We choose θ = −13 • [50] yielding mη ≈ 537 MeV and m η ≈ 1.15 GeV which is of sufficient accuracy for our purposes. 2 We neglect the small uncertainties of 0.04 in the determination of this coefficient for our discussion. We note however that recently Ref. [59] computed the coefficient to be 2.05, which suggests that uncertainties may in fact be larger. that studying one operator at a time, as usually done in the ALP literature, is insufficient to capture the full phenomenology of ALP models. Using this effective photon coupling the decay width to photons can then be written as where α is the electromagnetic fine-structure constant. We point out that, since we assume negligible tree-level interactions between ALPs and electrons, the ALP will decay dominantly into photons, i.e. Γ a tot ≈ Γ aγγ , and hence its lifetime will be very sensitive to its mass. Loop-induced decays into electrons/muons are negligible [44] and the hadronic decay channel a → 3π is relevant only for m a > 3m π , which is out of the mass range considered in this work.
For C GG = 0 the ALP will also couple to nucleons. These interactions can be determined by including protons and neutrons into the chiral EFT. Defining the nucleon coupling by L = ∂ µ a/(2Λ) C NN γ µ γ 5 N , where N is a nucleon field, we obtain 3 [61] C p + C n ≈ 16π 2 C GG , (2.12) Note that these relations imply that for an ALP coupling predominantly to gluons the effective proton coupling C p is about an order of magnitude larger than the effective neutron coupling C n .

Computation of K
One of the most sensitive laboratory probes of MeV-scale ALPs are searches for the rare decay of a charged kaon to a charged pion and an ALP, where the ALP escapes from the detector as an invisible particle [44,46]. Past experimental searches in this channel provide stringent constraints on ALP models [62,63], while upcoming experiments offer substantial discovery potential. For C W W = 0 the transition occurs via a loop-level FCNC, where the ALP is emitted from the internal W boson line of the penguin diagram. This transition generates an effective interaction [46] L = −g sd ∂ µ adγ µ P L s + h.c. , (2.13) with the effective coefficient (2.14) Here we have neglected the up-quark contribution and introduced the loop function With this effective interaction we can then determine the decay width to be where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2xy. We have neglected the hadronic form factor as it is close to unity [5,64,65].
In the case of gluonic ALP couplings, the ALP mixes with the pseudoscalar mesons as discussed in section 2.2 and will therefore participate in the K + decay through mixing effects. Taking into account that K + → π + a involves a change of isospin by 1/2, one expects the ∆I = 1/2 rule to apply [66]. By employing a chiral Lagrangian (see appendix B for calculational details), we can then approximately relate the ALP amplitude to the amplitude for K 0 → π + π − : Here A i and χ i describe the amplitudes and phases in the isospin decomposition of the kaon amplitudes according to the notation in Refs. [67,68]. The first term corresponds to the contribution of the K + → π + π 0 decay, which is usually omitted because of the A 2 /A 0 ≈ 1/22 suppression. However, for m a ≈ O(100 MeV) one finds θ aπ θ aη , θ aη , which can overcome this suppression. The resulting decay width reads [66,69] where τ are the kaon lifetimes, D 2 ππ ≈ 1/3 corrects for the absence of final-state interactions for π + a [70,71] and the three-momenta are | p X | = λ 1/2 (m 2 K , m 2 π , m 2 X )/(2m K ) in the lab frame. We emphasize that the computation presented above is only of approximate nature and is subject to uncertainties such as the η − η −mixing that can not be described accurately at leading order in a chiral EFT. We refer to Ref. [69] for a more detailed discussion.
We note that ALP-quark couplings would in general also give a contribution to the process K + → π + a. The resulting effect is however not easy to calculate in an EFT approach, because the loop-induced FCNC s → da is divergent [72]. In other words, it becomes difficult to reliably calculate experimental bounds and discovery prospects in the K + → π + a channel without a specific underlying model. To obtain a rough estimate for the case of flavour-universal couplings, one can replace the UV divergence by a leading logarithm [72]. Adopting this prescription we find that for an ALP that couples dominantly to quarks the most interesting regions of parameter space are already excluded by existing constraints. We will therefore not discuss this case in more detail and focus instead on the interactions with gauge bosons.  Figure 1: Differential distribution of m 2 miss for different K + decays taken from Refs. [48,73] and supplemented by the prediction for K + → π + a (purple) with C W W /Λ = 5·10 −4 TeV −1 and m a = 50 MeV, assuming a m 2 miss resolution of 0.001 GeV 2 . Here the regions labelled as "1" and "2" refer to the two signal regions of NA62, which are separated by the K + → π + π 0 background.

Experimental and Observational Constraints
In the following we will discuss the most important bounds from fixed-target experiments, collider searches and astrophysics. We will, however, not consider cosmology constraints as they are very model-dependent, see Ref. [30] for a detailed recent discussion. Instead we will focus on the general aspects of a light pseudoscalar particle interacting dominantly with SM gauge bosons independent of its precise cosmological history.

K
The fixed-target experiment E787 and its successor E949 have studied the decays of charged kaons at rest into a single charged pion and have placed constraints at 90% confidence level (C.L.) on the process K + → π + + X 0 [62,63], where X 0 denotes a new, long-lived particle.
To interpret the experimental bounds on the branching ratio in our context, we have to ensure that the ALP in fact escapes from the detector before decaying. We do this by including a probability factor that the ALP does not decay inside of the decay volume Here L det describes the characteristic size of the detector and l a = | p a |/(m a Γ a tot ) is the ALP decay length determined by its momentum | p a | in the lab frame. For E787 and E949 we take L det = 4 m and calculate | p a | given that the kaons decayed at rest. 4 We will also include prospects for NA62, a presently running experiment aiming to measure the ultra-rare K + → π + νν channel, which has an SM prediction of BR(K + → π + νν) = (8.4 ± 1.0) · 10 −11 [74][75][76], with 10% precision [48]. In this channel the νν pair escapes the detector without a trace, giving rise to a missing mass m 2 miss = (p K − p π ) 2 , where p K and p π are the kaon and pion four-momenta. When reinterpreting this search for K + → π + a we need to account for the fact that the SM process has a 3-particle final state. In this case the missing mass distribution is spread kinematically whereas for the 2-particle final state m 2 miss will be peaked. Note that, while NA62 cannot search for ALPs with mass m a ≈ m π due to prohibitively large backgrounds from K + → π + π 0 , it is sensitive to both the case m a < m π (signal region 1) and m a > m π (signal region 2).
As an example we plot the missing mass distribution m 2 miss in figure 1 assuming For this particular choice of mass this signal falls into the signal region 1 of NA62. One notices that despite K + → π + a having a much smaller branching ratio than K + → π + νν for these parameter points, the different kinematics imply that this ALP channel would have a visible impact on the measured distribution. However, as there are no official projections provided by NA62 concerning their sensitivity for the K + → π + a channel, we assume that NA62 will improve the experimental bounds by E787 and E949 by roughly an order of magnitude based on the information given in table 2 of Ref. [77]. 5 For eq. (3.1) we take | p a | = E a = 37 GeV and L det = 200 m in the case of NA62 [80].

Beam dump experiments
In beam dump facilities ALPs can be produced not only in flavour changing meson decays, but also via their direct couplings to photons and gluons, e.g. via the Primakoff process. Important constraints come from the E137 experiment at SLAC [81], where a search for ALPs with dominant photon couplings was performed. The original analysis was extended to the full ALP parameter space in Ref. [47], which uses the null result to obtain a bound at 95% C.L. . The same constraint also applies in the case of a dominant ALP-gluon coupling, for which the induced photon coupling (see eq. (2.10)) dominates the production in electron beam dump experiments.
Proton beam dump experiments also offer an additional promising search strategy for ALPs with gluon couplings. In Refs. [40,79] the constraints by CHARM have been revisited for the case that the production occurs through mixing with pseudoscalar mesons. The resulting production is found to be strongly enhanced compared to the case of couplings only to photons for ALPs with mass m a 100 MeV. Moreover, due to their higher beam energy, proton beam dump experiments are particularly sensitive to relatively large ALP masses and couplings, for which ALPs produced in an electron beam dump would decay before reaching the detector. However, the information provided in Refs. [40,79] is insufficient to allow for a reinterpretation of these bounds for ALPs with additional couplings to photons or under the assumption of different η and η mixing angles. We will therefore not include this bound in the following but stress that the sensitivity of proton beam dump experiments to ALPs with several couplings and large masses would be very interesting to explore in a proper simulation (see Refs. [16,55]).
In section 2.1 it was shown that the various ALP interactions in the electroweak sector are highly correlated and are therefore expected to arise simultaneously. In addition to the W ± boson and photon coupling, the coefficient C γZ can be relevant as it induces the Z → γ + a decay. The decay width for this process reads [42,82] A search in this channel was performed by L3 at LEP at the Z resonance. The null-results were converted into a bound on the branching ratio as a function of the photon energy, which for small m a corresponds to BR(Z → γa) 10 −6 at 95% C.L. [47,83].

Astrophysical Bounds
For MeV-scale ALPs that couple to photons and nucleons astrophysical probes also provide relevant constraints. In our discussion cooling constraints and the missing "ALP burst" from SN1987A as well as modifications of the stellar evolution of horizontal branch (HB) stars are relevant.

Cooling constraints from SN1987A
In the standard scenario of a core-collapse supernova most of its energy is released through neutrino emission, constituting the main cooling mechanism for the proto-neutron star, which forms as a result of the collapse. However, if additional, new weakly interacting particles exist with masses below O(100 MeV) these can also be produced in the SN core, which has a temperature of several tens of MeV. If this new energy loss mechanism is as efficient as (or even more efficient than) the energy loss through neutrinos, this would be in tension with the observed duration of the neutrino signal on Earth for the case of SN1987A. To make precise statements about the influence of a specific ALP on the supernova evolution one would in principle need to include all relevant effects in a detailed supernova simulation. Here we will instead use the approximate analytic criterion proposed by Raffelt [20] stating that the energy loss due to new particles should not exceed the neutrino luminosity L ν = 3 · 10 52 erg/s. 6 In order to determine the luminosity L a of ALPs produced in a proto-neutron star, we employ the ansatz from Refs. [23,85] and extend it to ALP-photon couplings. The luminosity is defined by a volume integral over the volume emission rate of ALPs within the neutrinosphere R ν , beyond which neutrinos stream freely, as ALPs produced in this region pose a new energy loss mechanism. However, for large ALP couplings the interaction will be so strong that ALPs produced in this region do not actually escape but get absorbed again, which is referred to as the "trapping regime". This is taken into account by including an optical depth factor in the integration characterizing the probability that an ALP produced in the core region reaches the radius R far , beyond which the ALP energy cannot get reprocessed efficiently into neutrino energy [85]. We provide more details on the luminosity formula and the ALP energy distribution contributing to the luminosity in appendix C.
For ALPs coupling dominantly to photons Primakoff conversion involving protons, γ + p → a + p, provides the main production mechanism [86]. In the case of ALP-gluon couplings bremsstrahlung in nucleon scattering N + N → N + N + a is most relevant [20,87]. Processes involving electrons are negligible as their phase space is Pauli-blocked in a supernova core [88]. Combining all of these points, we obtain for the ALP luminosity (3. 3) The integration variables p a and p γ describe the ALP and photon momenta, Γ γ→a and Γ a the production rates through Primakoff and bremsstrahlung, T the temperature as a function of radius r and τ ≡ τ (m a , ω, r, R far ) is the optical depth. Following Ref. [23], we set R far = 100 km. For the Primakoff process we have used that the photon energy ω γ is equivalent to the ALP energy ω a ≡ ω as the proton mass is much larger than all other energies involved [25,86]. In the integration we have also taken into account that the ALP/initial photon energy has to be at least m a to produce an ALP of this mass. The effects of sizeable ALP masses, i.e. m a > 1 MeV, are estimated by a phase space factor β = 1 − m 2 a /ω 2 [22,23]. We now go into more detail concerning the different quantities entering eq. (3.3): • We use the axion absorptive width Γ a for nucleon bremsstrahlung given in eqs. (4.2) and (4.3) in Ref. [23], where we perform the replacement f a → Λ. For the one-pion exchange correction factor γ h we use the parameter values for an electron fraction of Y e = 0.1 given in table 1 in Ref. [89]. We stress, however, that the different parameter sets are very similar at higher densities, where most of the ALP production and capturing will occur. For the factor γ p , which corrects for the finite pion mass and nucleon degeneracy, we take into account the fraction of nuclei that actually participate in the scattering. For example, an ALP with gluon couplings interacts much more strongly with protons than with neutrons. In this case we use a density of Y p ρ with the proton fraction Y p = 0.3 as input parameter.
• The transition rate of a (massless) photon with energy ω into an ALP of the same energy via the Primakoff process is given by [21] Γ γ→a = 16παC eff where κ is the screening scale as the Coulomb potential only has a finite range in a plasma [86]. Note that this transition rate is averaged over the photon polarizations, which is the origin of the factor 2 appearing in eq. (3.3). Taking into account that the protons are partially degenerate, the screening scale reads [88] with the effective number of proton targets n eff p . For simplicity we assume that n eff p = 0.5 n p = 0.5 Y p ρ within the core region, i.e. r ≤ 10 km, and n eff p = n p outside. This degeneracy has, however, only slight effects.
• The optical depth includes contributions from inverse Primakoff processes, ALP decays into photons and inverse nucleon bremsstrahlung. Combining these, it reads (see appendix C for the general definition) where we included a factor 2 due to the two photon polarizations and a Lorentz factor γ = ω/m a [22,90]. We also split the integration into r < R ν and r > R ν as the temperature and density profiles used do not go that far [85]. This will also give a conservative estimate of the effects in the outer region of the star as the temperature and density profile are kept constant beyond r = R ν instead of letting them decrease.
• The temperature and density profile used in this work is the fiducial one as given in Ref. [85]. We have also checked the exclusions contours for the profiles by Fischer (11 M and 18 M ) [91] as well as Nakazato (13 M ) [92] with the parameters as given in table 2 of Ref. [85] and have found the fiducial one to give the most conservative bound overall. For the remainder of this work we will therefore use the fiducial profile.
Let us now compare our results with those from the literature. It is important to stress that even when using the same cooling criterion, differences in the calculational approaches and the assumed temperature and density profiles typically lead to sizeable differences in the resulting bounds, which demonstrates the difficulty of extracting reliable supernova constraints. In the left panel of figure 2 we consider dominant ALP-photon couplings. We observe that the lower boundary, where ALP production becomes efficient enough to alter the duration of the neutrino signal (the free-streaming regime), is similar to the ones found in other recent calculations (Dolan et al. [47] and Lee [22]). In the trapping regime, however, for which we calculate a detailed optical depth factor, we obtain a weaker constraint than calculations based on the opacity.
In the right panel of figure 2, we focus on the case of a dominant ALP-gluon coupling, for which bremsstrahlung production dominates over Primakoff processes. For the comparison we identify g 2 aN N from Ref. [22] and 1/(2f a ) 2 from Ref. [23] with (C 2 p Y p + C 2 n Y n )/Λ 2 . While our determination of the trapping regime for small ALP masses gives similar results . We have taken the brown curve from Ref. [47], the red curve from Ref. [22], the yellow one from Ref. [42] (originally computed in Ref. [93]), the black as well as the green one from Ref. [23] and the purple coloured area corresponds to our result.
in this case, our bound differs for larger ones as well as in the free-streaming regime. Most notably, we find that the supernova bound is essentially independent of the ALP mass for m a 10 MeV, consistent with the fact that most ALPs produced in the supernova have an energy well above 10 MeV (see appendix C). We also point out that Ref. [94] obtains a stronger bound on ALP-nucleon interactions for small ALP masses (in particular in the trapping regime), but does not consider ALPs with masses m a > 1 MeV.

ALP burst from SN1987A
In the coupling regime below the supernova cooling constraint ALPs are still produced within the supernova core, but they do not lead to significant energy loss. In Ref. [51] it was pointed out that these ALPs are sufficiently long-lived to escape the supernova but would still decay before reaching Earth. This would have resulted in a (delayed) gamma ray burst that was, however, not observed. This missing ALP burst gives an important constraint for ALP-photon couplings below the cooling constraint.

Number count of events
In addition to the two constraints discussed above, there is another bound on ALP-nucleon interactions often quoted in literature that excludes much more parameter space above the SN1987A cooling bound [95]. In this strong coupling regime one has to consider that ALPs do not stream out of the core region but instead form a black body sphere, the axionsphere (analogous to the neutrinosphere) [96]. If the ALPs radiated from this sphere reached the Earth, they would induce more events in the Kamiokande experiment than observed, and hence such large couplings are excluded [95]. The ALPs that we consider, however, are too short-lived to reach the Earth so that this constraint does not apply.
What could play a role instead are ALPs produced on the axionsphere, leaving the supernova and then decaying to photons before reaching the Earth [97]. The absence of a gamma ray signal could then lead to additional constraints similar to the ALP burst discussed before. However, in the coupling range of interest and for ALP masses above a few MeV the decay length is too short to reach the effective radius of R eff ≈ 3 · 10 10 m needed to leave the supernova [51,98]. Of course, even in this case one would need to ensure that not too much energy is deposited in the mantle and envelope as most of the energy is liberated by neutrinos [97,99,100]. A careful study of these issues is beyond the scope of this work.

HB stars
The existence of light ALPs would also influence stellar evolution. Strongly affected by ALP-photon couplings are horizontal branch (HB) stars, i.e. stars that have entered the helium burning phase [101]. An additional energy loss through ALPs leads to a faster contraction of the core region and hence an increase in temperature [102]. This then results in a faster burning of the helium fuel and reduces the lifetime of horizontal branch stars. This influence has been studied in the context of globular clusters, where one can determine the relative number of horizontal branch stars and red giants, which are not as affected by the Primakoff production of ALPs due to their degenerate cores [103]. As the observed ratio is within 10% of the prediction, new energy loss contributions can be constrained [20,21]. The resulting bounds on the ALP-photon coupling have been computed in a simplified manner in Ref. [25], which we will use in the following. Note that for ALPs that couple to gluons the induced ALP-photon interaction dominates over the ALP-nucleon interaction, because the resulting Compton cross section is strongly momentum suppressed [20].

Results
In this section we combine the various constraints discussed above to determine the viable regions of parameter space and explore the prospects for probing ALPs with NA62. We consider first the case that ALPs couple only to electroweak gauge bosons, then explore the case of gluonic couplings and finally study the impact of varying all couplings simultaneously.

Couplings to electroweak gauge bosons only
Let us first consider the case of dominant couplings to SU (2) L gauge bosons, i.e. the case that all Wilson coefficients other than C W W can be neglected. The corresponding constraints are shown in the left panel of figure 3. One finds that the combined constraints from HB stars, SN1987A, L3, E137 and E787/E949 exclude essentially all ALPs with mass m a 100 MeV (unless C W W /Λ is much smaller than what is considered in this figure) and in particular the entire parameter region that can be probed with NA62. This conclusion relies crucially on the fact that couplings to SU (2) L gauge bosons induce both couplings to W ± bosons (which lead to constraints on FCNCs from E787/E949) and couplings to Z bosons and photons (which result in all the other constraints) after electroweak symmetry breaking. Previous studies have often found weaker constraints by either neglecting the ALP-photon coupling (as done e.g. in Ref. [44]) or assuming the presence of another light species that allows for ALPs to decay dominantly invisibly (see Ref. [46]).
However, it is also possible to relax the constraints without leaving the EFT framework introduced above. Since C γγ = C W W + C BB it is possible to shift the constraints that depend on the ALP-photon coupling relative to the constraints that depend on the ALP-W coupling by considering the case C BB = 0. As an example, we consider in the right panel of figure 3 the case that C BB = 10 C W W , such that the value of C γγ corresponding to a given value of C W W is substantially enhanced. The bounds that depend on the ALPphoton coupling are hence shifted downwards relative to the ones depending on FCNCs, so that the parameter space probed by NA62 falls into the trapping regime of SN1987A. By considering simultaneously the ALP couplings to SU (2) L and hypercharge gauge bosons, it is hence possible to open up parameter space where MeV-scale ALPs may be discovered in the laboratory.
We note that in the case that C BB = 10 C W W there are no constraints from E787/E949 (and no sensitivity of NA62) for m a 100 MeV, because the enhanced photon coupling implies that ALPs in this mass range would no longer escape from the detector before decaying. Instead, they would contribute to the channel K + → π + γγ, where the sensitivity is reduced by large SM backgrounds and branching ratios as large as O(10 −6 ) are allowed [80,104].
Finally, we point out that constraints could potentially be relaxed even further in the case of destructive interference, if C W W ≈ −C BB such that the ALP-photon coupling is suppressed. The constraints that depend on this coupling would then be shifted upwards relative to the constraints that depend on FCNCs. However, given the strength of the supernova constraints, a very precise cancellation would be required in order to evade these constraints entirely, making this solution less attractive.

Couplings to all Standard Model gauge bosons
The interplay of different couplings becomes even more interesting when including the ALPgluon coupling in the discussion. Figure 4 shows the resulting constraints as a function of m a and C GG for four different combinations of Wilson coefficients. The figure focuses on m a < m π (corresponding to signal region 1 in NA62), while the case m a > m π (signal region 2) is discussed in appendix D. The top-left panel corresponds to the case where only the gluon coupling is relevant and all other couplings can be neglected. As discussed in section 2, the low-energy phenomenology in this case is determined by the effective ALPphoton coupling, the effective ALP-nucleon coupling and the FCNCs induced by meson mixing. The ALP-photon coupling leads to relevant constraints from E137 and HB stars. The constraint from SN1987A depends on both the photon and the nucleon coupling, but the latter dominates the phenomenology because the high nuclear density in a supernova core leads to an enhancement of bremsstrahlung processes over Primakoff conversion. As a result, the cooling constraints are sensitive to much smaller Wilson coefficients than for the case of ALPs coupling to electroweak gauge bosons. At the same time, the trapping regime extends to much smaller couplings, opening up large regions of unconstrained parameter space.
For the case of ALP-gluon couplings, constraints on K + → π + a from E787/E949 are therefore highly relevant and the prospects of NA62 look very promising. However, it is important to keep in mind that these constraints, which are obtained from the ALP-meson mixing, are subject to hadronic uncertainties, see eq. (2.18) and the surrounding text. It is therefore interesting to consider the case where the effective coupling to W bosons gives a relevant contribution to the flavour-changing processes. One example is shown in the top-right panel, which corresponds to the case C GG = C W W (keeping C BB = 0). In this case the constraint from E787/E949 shifts downward and become largely independent of hadronic uncertainties. At the same time the constraint from E137 shifts upwards, because C GG and C W W contribute to the effective photon coupling with opposite sign and therefore cancel partially. The bound from SN1987A is unaffected by the ALP-photon coupling and therefore stays the same in all panels. Furthermore, we now also obtain constraints from L3, which rules out additional parameter space for m a ≈ 100 MeV.
In the bottom-left panel, where we take C GG = C W W = C BB , we observe an approximate cancellation in the effective photon coupling (see eq. (2.10)) leading to a strong suppression of the bound from E137 and moving the HB star bound to larger couplings. Moreover, the E137 bound now depends sensitively on the contribution to the effective photon coupling from ALP-meson mixing, which leads to a more complicated dependence on the ALP mass and larger hadronic uncertainties. For this specific combination of Wilson coefficients the constraints on the parameter space resulting from the ALP-photon interaction hence become less important. We point out that it is quite plausible that the Wilson coefficients C GG , C W W , and C BB would be comparable in size, for example if they are generated from the contribution of heavy new particles through triangle diagrams [45].
As an alternative, we consider in the bottom-right panel a more hierarchical coupling structure, with C GG = 0.1C W W and C BB = 0. As expected, the bounds from E787/E949 and L3 become significantly stronger in this case and the sensitivity of NA62 extends almost down to the parameter region excluded by SN1987A. For even larger hierarchy between C GG and C W W it may even be possible to explore the entire trapping regime with laboratory experiments. The key point, however, is that the individual bounds depends in different ways on the ALP couplings and can therefore shift relative to each other. A simple one-operator study is insufficient to capture all of these dependencies and may give a misleading impression of the status of ALP models. A more careful analysis instead reveals interesting regions of parameter space where ALPs may be discovered with ongoing or future experiments.

Conclusions
The interest in the study of hidden sectors, i.e. new light particles with extremely weak interactions, has increased substantially over recent years as a reaction to the absence of evidence for new physics at the electroweak scale. The growing theoretical and experimental efforts have paid particular attention to ALPs, which can arise as Pseudo-Nambu-Goldstone bosons from a broken global symmetry. In the present work we have studied MeV-scale ALPs in an EFT set-up with a particular focus on couplings to SM gauge bosons. We improve on the common approach to focus on a single effective operator at a time by taking into account all relevant low-energy interactions and their correlations. Relevant effects are found to stem from ALP interactions with photons, W ± bosons and nucleons, as well as ALP-meson mixing. We study the interplay of the resulting constraints and map out the phenomenology of ALPs with several types of interactions.
We are particularly interested in the prospects of searching for MeV-scale ALPs in the laboratory by performing measurements in the K + → π + + inv. channel, for which the presently running NA62 experiment promises to significantly improve experimental sensitivity. We compare the projected reach of NA62 with complementary constraints on the ALP parameter space, most importantly from SN1987A. We improve upon existing calculations by consistently including both ALP-photon and ALP-nucleon couplings. As shown in figure 2, our approach leads to a more conservative estimate of the trapping regime for ALP-photon couplings, while for ALP-gluon couplings our results are largely consistent with recent estimates. We furthermore discuss the details of ALP-meson mixing and estimate the resulting contribution to flavour-changing processes.
As a first application of our framework we have focussed on ALPs that interact dominantly with electroweak gauge bosons. The case that ALPs couple only to SU (2) L gauge bosons is already severely constrained and almost the entire parameter space relevant for NA62 is already excluded. However, as soon as a simultaneous coupling to hypercharge gauge bosons is considered, the different constraints can be shifted relative to each other. In particular, by enhancing the coupling to photons it is possible to extend the trapping regime and circumvent the cooling bound from SN1987A, see figure 3.
We then studied the entire range of interactions with gauge bosons by also considering the ALP-gluon coupling. The resulting ALP-nucleon interactions push the cooling constraint from SN1987A to much smaller couplings, opening up large regions of unconstrained parameter space that can be explored with laboratory experiments like NA62 (figure 4). By varying the relative magnitude of the different Wilson coefficients, we illustrate how the individual bounds and prospects (and hence the overall phenomenology) depend on the different couplings.
Our results demonstrate that it is essential to simultaneously consider multiple effective ALP interactions, which are generically expected to be generated simultaneously in many ultraviolet completions. Their correlations have an important influence on the interplay of different probes and the viable regions of parameter space. Focusing on single operators may introduce significant biases, for example by overestimating the strength of certain exclusion limits. Nevertheless, we still find that there are many strong constraints on the ALP parameter space, which can be improved with further theoretical and experimental efforts. In particular, it would be very interesting to include MeV-scale ALPs in a supernova simulation to study their impact in more detail and obtain more robust constraints from SN1987A. Likewise, we have left a detailed study of constraints and sensitivity projections for proton beam dump experiments to future work.

A Details on ALP interactions with hadrons
In this appendix we provide details on the computation of the ALP mixing with pseudoscalar mesons. We start from the second line of eq. (2.6) and map it onto a chiral Lagrangian at leading order. This yields [105][106][107][108] Here M 0 accounts for the mass term induced by the U (1) A breaking, B 0 = m 2 π /(m u + m d ) and U = e iφ/F 0 with the matrix φ containing the fields of the pseudoscalar mesons The mass matrix is given by M q (a) = e −iκqa/2fa M q e −iκqa/2fa with 1/f a = −32π 2 C GG /Λ. We choose κ q = M −1 q /Tr[M −1 q ] to ensure that there is no mass mixing between the ALP and the octet mesons [41]. The covariant derivative is defined by with the spurions reading v µ = −eQ q A µ , We first rotate η 8 and η 0 into η and η where we have introduced the short-hand notation c(θ) ≡ cos(θ) and s(θ) ≡ sin(θ). We adopt the value θ = −13 • from Ref. [50]. To diagonalize the η − η −sub-mass-matrix, we then have to take M 0 ≈ 1.05 GeV resulting in m η ≈ 537 MeV and m η ≈ 1.15 GeV sufficient for our intended accuracy. Expanding now eq. (A.1), one obtains both kinetic and mass mixing contributions. Following the notation of Ref. [49], we write these as with P = (a, π 0 , η, , η ). The matrices read where δ I = (m d − m u )/(m u + m d ) and For the remaining expressions we refer to eqs. (2.8) and (2.9). Diagonalisation of the kinetic as well as mass matrix results in the mixing angles given in eq. (2.7) at leading order.
B Mixing contributions to K + → π + a The hierarchy encountered in kaon decays, namely Γ(K 0 → π + π − ), Γ(K 0 → π 0 π 0 ) Γ(K + → π + π 0 ), is usually explained by considering the leading order chiral Lagrangian for ∆S = 1 transitions [67,68,109] L ∆S=1 ⊃ G 8 F 4 0 Tr Here λ i denote the Gell-Mann matrices, L µ = iU † D µ U and the indices attached to L µ denote specific matrix elements. From this Lagrangian, we obtain for the amplitude where A i and χ i describe the amplitudes and phases in the isospin decomposition of the kaon amplitudes [67,68]. We have expanded each amplitude in A 0 /A 2 ≈ 22 and additionally in m 2 K /m 2 π at leading order. As pointed out in Ref. [66] the enhanced term proportional to G 8 in eq. (B.1) also contributes to the off-shell amplitudes These amplitudes are enhanced compared to iM(K + → π + π 0 ) and are therefore also included. The K + → π + a amplitude can then be determined by taking all mixing contributions into account coinciding with eq. (2.17).

C Details on the computation of the ALP luminosity in SN1987A
The general expression for the ALP luminosity involves a volume integral over the energy loss rate Q, where also the probability e −τ of an ALP escaping the neutrinosphere is taken into account [23,85]: Here the energy loss rate Q (energy per volume and unit time) is defined by [87] where p j (ω j ) denote the three-momenta (energies) of the ALP (j = a), the initial states (j = i) and final states (j = f ). Moreover, f (ω j ) are the phase-space occupation numbers defined by n j = g j d 3 pf (ω j )/(2π) 3 with the number density n j and degeneracy factor g j of the particle species j. Note that a factor of 1 + f (ω j ) is taken for bosonic final states (stimulated emission), whereas 1−f (ω j ) applies to fermions (Pauli blocking). Furthermore, S denotes a symmetry factor for identical states in the initial or final state and |M | 2 is the squared matrix element which is summed over initial and final state spins and polarizations. The optical depth τ is defined by an integral over the inverse mean free path l = β/Γ abs [87,90]: with the absorption rate Γ abs for a process a + i → f reading [110]

(C.4)
We stress that the ALP energy integration in eq. (C.2) also affects the optical depth. In the case of the Primakoff process, where one integrates over the photon momentum in eq. (3.3), we exploit that the photon energy is equal to the ALP energy. The different rates have already been computed in the literature and simply need to be combined. For ALP production through bremsstrahlung in nucleon scattering and Primakoff processes we have [20,21] Q brems = d 3 p a (2π) 3 ωΓ a e −ω/T β , ωΓ γ→a e ω/T − 1 β , where ω a ≡ ω denotes the ALP/photon energy and the phase space factors β account for sizeable ALP masses. With these results, we can also determine the different absorption rate formulas 7 Γ brems abs = Γ a , Γ prim abs = 2Γ γ→a , Γ decay abs = γ −1 Γ aγγ , where an inverse Lorentz factor γ = ω/m a is included in the last line as Γ aγγ is computed in the ALP rest frame. Note that we neglected the Bose factors for the photon final states for simplicity as these barely influence the results. On the left we focus on a dominant ALP-photon coupling and on the right on an ALP-gluon coupling.
To illustrate the process of ALPs cooling the proto-neutron star, we provide in figure 5 the distribution of ALP energies contributing dominantly to the luminosity as a function of radius for ALP-photon couplings (left) and ALP-gluon couplings (right). Here we take m a = 1 MeV and consider two different coupling values in each case, corresponding to the upper and lower edge of the excluded parameter region. We make the following observations: For smaller coupling values most of the ALPs escaping the neutrinosphere get produced in the inner core region, where both temperature and density are largest, and therefore also have energies of up to O(100 MeV). On the other hand, those ALPs with rather large couplings are only able to leave the neutrinosphere when they are produced close to its edge as the trapping in the inner part is too strong. But even these ALPs have energies 10 MeV on average, making it clear why the supernova cooling bounds in figure 2 become independent of the ALP mass for m a 10 MeV.

D Prospects of NA62 for signal region 2
In this appendix we provide the constraints and prospects for the second signal region in NA62, which is sensitive to m a > m π . As one can see from figure 6, this signal region is not as promising as the first one since it is largely excluded by E137 and L3 combined. Moreover, the general structure of the bounds here is more complex, because the mixing contributions between the ALP and pseudoscalar mesons become crucial and lead to cancellations in the ALP-photon coupling.