Signals of the QCD axion with mass of 17 MeV/c^2: nuclear transitions and light meson decays

The QCD axion remains experimentally viable in the mass range of O(10 MeV) if (i) it couples predominantly to the first generation of SM fermions; (ii) it decays to $e^+ e^-$ with a short lifetime $\tau_a\lesssim 4\times 10^{-14}\,$s; and (iii) it has suppressed isovector couplings, i.e., if it is piophobic. Interestingly, these are precisely the properties required to explain recently observed anomalies in nuclear de-excitations, to wit: the $e^+e^-$ emission spectra of isoscalar magnetic transitions of $^{8\!}$Be and $^{4\!}$He nuclei showed a"bump-like"feature peaked at $m_{e^+e^-}\sim 17$ MeV. In this article, we argue that on-shell emission of the QCD axion (with the aforementioned properties) provides an extremely well-motivated, compatible explanation for the observed excesses in these nuclear de-excitations. The absence of anomalous features in other measured transitions is also naturally explained: piophobic axion emission is strongly suppressed in isovector magnetic transitions, and forbidden in electric transitions. This QCD axion hypothesis is further corroborated by an independent observation: a $2-3\,\sigma$ deviation in the measurement of $\Gamma(\pi^0\to e^+e^-)$ from the SM theoretical expectation. This article also includes detailed estimations of various axionic signatures in rare light meson decays, which take into account contributions from low-lying QCD resonance exchange, and, in the case of rare Kaon decays, the possible effective implementations of $\Delta S=1$ octet-enhancement in chiral perturbation theory. These inherent uncertainties of the effective description of the strong interactions at low energies result in large variations in the predictions for hadronic signals of the QCD axion; in spite of this, the estimated ranges for rare meson decay rates obtained here can be probed in the near future in $\eta/\eta^\prime$ and Kaon factories.


Introduction
The past decade has seen a resurgence of interest in the phenomenology of new light particles with feeble interactions with the Standard Model (SM) [1][2][3]. Motivations have been varied, spurred from the growing belief that dark matter might be part of a more complex dark sector with additional matter and force carriers [4][5][6][7][8], but also because light dark sectors could be parasitically explored in the broader U.S. and worldwide neutrino program [9][10][11][12]. Experimental signatures of dark sectors are being searched for by a diverse suite of experiments ranging from beam dumps/fixed targets to meson factories 1 . This effort drew on the legacy of an earlier, very active period of "intensity frontier" experiments initiated in the 1970's. This earlier period, however, was driven partly by studies of hadronic and neutrino physics, and partly by searches for the Higgs boson and the QCD axion. Indeed, special region of parameter space, the phenomenology of the QCD axion is no longer dominated by its isovector couplings; instead, it is largely determined by its isoscalar mixings with the η and η mesons. As such, it inherits the same strong dependence of the η and η on higher order terms in the chiral expansion, and its hadronic couplings suffer from O(1) uncertainties.
Despite these large uncertainties stemming from χPT, it is still possible to parameterize the dependence of a variety of hadronic signatures of the axion in terms of its isovector and isoscalar mixing angles, while remaining agnostic about their magnitudes. The usefulness of such parameterization is manifest when confronted with experimental data, not only in constraining the axion's hadronic mixing angles, but also in interpreting experimental anomalies as potential signals of the QCD axion. This will be the underlying philosophy of this study 4 .
Complementary, the underlying motivation for this study is a combination of the longstanding puzzle posed by the strong CP problem, and three independent experimental anomalies. The first two refer to bump-like excesses observed in specific magnetic transitions of 8 Be and 4 He nuclei via e + e − emission, with (naïve) significances of 6.8 σ [42] and 7.2 σ [43], respectively. The third anomaly is related to the persistently high central value observed for the width Γ(π 0 → e + e − ), whose most recent and precise measurement, performed by the KTeV collaboration in 2007 [44], showed a discrepancy from its theoretical expectation in the SM at the level of ∼ 2 − 3.2 σ [47][48][49][50]. In combination, these anomalies point to a common BSM origin: a new short-lived boson with mass of ∼ 16 − 17 MeV, coupled to light quarks and electrons, and decaying predominantly to e + e − . As an ad hoc explanation, there are only two possibilities for the spin-parity of this hypothetical new boson: it can either be a pseudoscalar (J P = 0 − ), or an axial-vector (J P = 1 + ), in order to simultaneously account for these three excesses 5 . Further constraints push these two possibilities into peculiar regions of parameter space, which may require contrived and/or baroque UV completions 6 . At face value neither of them is particularly compelling, leading many to believe that these anomalies are either the result of experimental systematics and/or poorly understood SM effects. In our opinion, this illustrates the paradoxical predicament of the light dark sector intensity frontier program: the generic models it seeks to discover or rule out are not strongly motivated, and, at least historically, it has been the case that experimental excesses without theoretically compelling interpretations tend to be the 8 Be anomaly cannot be emitted in the 0 − → 0 + transition of 4 He, nor does it contribute non-negligibly to Γ(π 0 → e + e − ). In [53], Feng et. al proposed an alternative explanation of the 4 He anomaly, whereby the e + e − excess stems from the de-excitation of the overlapping 0 + nuclear state. Recently, [54] argued that the protophobic vector boson hypothesis is excluded as an explanation of the 8 Be anomaly. 6 For instance, in the axial-vector case, the model building required to circumvent stringent bounds from electron-neutrino scattering restricts the axial-vector couplings of the 1 + state to light quarks to satisfy g A u = −2 g A d [55,56]; axial-vector models also typically require many ad hoc degrees-of-freedom to cancel gauge anomalies in the UV. In the axion case, in order to suppress a − π 0 mixing, the PQ charges of the up and down quarks must satisfy q u PQ = 2 q d PQ , with (nearly) vanishing PQ charges for the other quarks. Such flavor alignment, combined with the fact that fP Q ∼ O(GeV), requires nontrivial UV completion at the weak scale, see [45]. received with strong skepticism.
Fortunately, this predicament might not be warranted here. Nuclear transitions via axion emission and (modified) rare meson decays are smoking-gun signatures of the QCD axion which have been predicted over three decades ago [57][58][59][60][61][62][63][64]. The fact that some of these signatures have appeared in 8 Be, 4 He, and π 0 decays, and can be consistently explained by a QCD axion variant which remains experimentally viable (albeit with peculiar properties of electrophilia, muophobia and piophobia), should be taken with cautious optimism.
After a brief overview of the most relevant properties of the piophobic QCD axion in Sec. 2, we obtain the parameter space of axion isoscalar couplings favored by the 8 Be and 4 He anomalies, and, taking into account nuclear and hadronic uncertainties, show that they significantly overlap, favoring the QCD axion emission hypothesis as a single explanation of both anomalies (Sec. 3). We then turn to axion signals in rare meson decays. In Subsec. 4.1, we obtain the the parametric dependence of η/η di-electronic decays on the axion's isoscalar mixing angles. In Subsec. 4.2, we calculate the rate for axio-hadronic decays of the η and η mesons in the framework of Resonance Chiral Theory, an effective "UV completion" of χPT that incorporates low-lying QCD resonances and extents the principle of Vector Meson Dominance. Finally, in Sec. 5, we investigate various axionic decays of charged and neutral Kaons, considering distinct possible implementations of octet enhancement in χPT and their effect on axionic Kaon decay rates. We conclude in Sec. 6.

Brief overview of the piophobic QCD axion
Generic models of the QCD axion with mass of ∼ 16 − 17 MeV are largely excluded. However, as investigated in [45], all experimental constraints can be avoided in this mass range if the axion satisfies a few specific requirements: (i) it must be short-lived (τ a 0.4 × 10 −13 s), and decay predominantly to e + e − in order to avoid limits from beam dump and fixed target experiments, as well as constraints from charged Kaon decays such as K + → π + (a → γγ , invisible); (ii) the PQ charges of 2 nd and 3 rd generation SM fermions must vanish or be suppressed, in order to avoid limits from the muon anomalous magnetic dipole moment, (g − 2) µ , and from upper bounds on radiative quarkonium decays: J/Ψ, Υ → γ (a → e + e − ); (iii) the a − π 0 mixing must be suppressed, θ aπ O(10 −4 ), in order to respect upper bounds on Br π + → e + ν e (a → e + e − ) .
A simple phenomenological IR model realizing the requirements above can be easily incorporated in the post-EWSB SM Lagrangian by ascribing axionic phases to the masses of the up-quark, down-quark, and electron: Importantly, no additional operators should be present in this specific basis, such as derivative couplings of the axion to quark axial-currents, or the usual linear coupling of the axion to the gluon dual field strength operator.
In this IR model, requirement (ii) mentioned above has been imposed by fiat. Requirement (i) follows from the axion's coupling to e + e − , which dominates its decay width: For m a ∼ 16 − 17 MeV, existing bounds on the electron's PQ charge are very mild, limiting its range to 1/3 |q e PQ | 2. The upper bound is set by KLOE's 2015 search for visibly decaying dark photons [65], whereas the lower bound is set by the 2019 results from CERN's SPS NA64 fixed target experiment [66,67], constraining the axion lifetime to τ a 4 × 10 −14 s. The sensitivities of future experiments to the axion's electronic couplings (such as fixed targets and e + e − colliders) have been explored in [45].
From (2.1) and standard χPT at leading order, the axion mass is given by: It follows then that for q u PQ /2 = q d PQ = 1 and m a = 16.7 MeV, the axion decay constant is f a 1030 MeV. We will benchmark m a , f a , q u PQ , and q d PQ to these values for the remainder of this article.
For generic parameter space of QCD axion models, the quark mass hierarchy m u,d m s typically induces a hierarchy of axion-meson mixing angles, θ aπ θ aη , θ aη , resulting in the isovector couplings of the axion dominating its experimental signatures. This is not the case for the piophobic axion we are considering. Here, the a − π 0 mixing angle, to leading order in χPT, is given by: which, after taking q u PQ /2 = q d PQ = 1 and m u /m d = 0.485 ± 0.027 from [68], results in: It is clear from (2.6) and (2.7) that the axion's piophobia is the result of an accidental cancelation in χPT's leading order contribution to θ aπ . This cancelation stems from the near numerical coincidence between m u /m d and q d PQ /q u PQ = 1/2.
Unfortunately, χPT's prediction (2.7) alone is not precise enough to be useful. We instead have resort to observation to determine the allowed range for θ aπ with better precision. This can be achieved by requiring that the 3.2 σ excess in KTeV's measurement of Γ(π 0 → e + e − ) [44] be the result of π 0 − a mixing, which yields [45]: Given the suppressed value (2.8) for θ aπ , this model features an atypical hierarchy of mixing angles, θ aπ θ aη , θ aη , which results in the isoscalar couplings of the axion dominating its experimental signatures. This aggravates the loss of χPT's usual predictive power in axion phenomenology-given its state-of-the-art, χPT cannot numerically pin down the isoscalar mixing angles θ aη , θ aη with good accuracy. As argued in [45], θ aη , θ aη receive O(1) contributions from O(p 4 ) operators in the chiral expansion, many of which have poorly determined Wilson coefficients.
Any substantive theoretical progress in better determining the axion's hadronic couplings is unlikely to be accomplished anytime soon. Indeed, such efforts might be superseded by future experimental results which will be able to either exclude or narrow down the preferred ranges for the axion's isoscalar couplings. With this in mind, in this study we choose to remain agnostic about their magnitude, and instead simply parameterize the physical axion current as 7 : The axionic field a and the neutral meson degrees of freedom π 3 , η ud , and η s in (2.10) mix amongst themselves to yield the physical degrees of freedom (i.e., the mass eigenstates) a phys , π 0 , η, and η . In particular, the implication of (2.9) is that any strong or weak process involving the currents in (2.10) will lead to a corresponding axion signature for which one of the neutral mesons in the amplitude gets replaced by a phys properly weighted by the appropriate mixing angle. With the parameterization in (2.9), it is straightforward to obtain the axion's couplings to photons and nucleons. Specifically, below the QCD confinement scale, the electromagnetic anomaly of the physical axion current (2.9) leads to: 11) 7 We omit the dependence of (2.9) onēγµγ5e, which has no bearing on the axion-meson mixing angles.
which, combined with (2.2), yields the axion decay width and branching ratio to two photons: , (2.13) The axion's contribution to (g − 2) e stemming from its couplings to electrons and photons has been worked out in [45]. Finally, expressing the axion nuclear couplings generically as: the parameterization in (2.9) yields the following isovector and isoscalar axion-nucleon couplings, respectively: Above, N is the nucleon isospin doublet, m N is the nucleon mass, and ∆q quantifies the matrix elements of quark axial-currents in the nucleon via 2 s µ ∆q = N |qγ µ γ 5 q |N , with s µ the nucleon spin-vector. The combination in (2.15a) is well determined from neutron β-decay, ∆u − ∆d = g A 1.27 . In the following, we will use (2.15) to fit the recent 8 Be and 4 He anomalies, and (2.9) to obtain various rare meson decays.

Nuclear Transitions
One of the smoking-gun signatures of axions in the mass range O(keV − MeV) are magnetic nuclear de-excitations via axion emission [57][58][59]. Indeed, such signals have been extensively searched for during the 1980's [82][83][84][85][86][87][88][89][90][91]. However, since the energy of typical nuclear transitions ranges from a few keV to a few MeV, past searches did not place meaningful bounds on axions heavier than m a 2 MeV. Recently, the MTA Atomki collaboration led by A. Krasznahorkay reported on the observation of bump-like excesses in the invariant mass distribution of e + e − pairs emitted in the de-excitation of specific states of 8 Be and 4 He nuclei [42,43]. The energy difference ∆E between the nuclear levels involved in these particular transitions is atypically high, a priori allowing on-shell emission of particles as heavy as ∼ 17 − 18 MeV. Furthermore, consistent with the allowed values of angular momentum and parity carried away by the axion (J P = 0 − , 1 + , 2 − , 3 + , ...), these excesses appeared in magnetic (but not electric) transitions. Also, consistent with the emission of a piophobic axion, these excesses were observed only in predominantly isoscalar (but not isovector ) transitions. Axion emission rates for the magnetic dipole transitions of 8 Be have already been worked out in [45]; we briefly review the main results here to make this section self-contained. We then estimate the expected axion emission for the magnetic monopole transition of 4 He investigated by Krasznahorkay et al., and show that the reported excess rates for both nuclei favor the same range of axion isoscalar mixing angles.

Evidence for the QCD axion in 8 Be transitions
In [42], the MTA Atomki experiment selectively populated specific excited states of the 8 Be nucleus by impinging a beam of protons with finely-tuned energy on a 7 Li target. They then measured the energy and angular correlation of e + e − pairs emitted in de-excitations of these states to the ground state of 8 Be. From these measurements they were able to reconstruct final state kinematic variables, such as the invariant mass of the e + e − pair, m e + e − . The nuclear levels of interest de-excited to the ground state via magnetic dipole (M1) transitions: Above, 8 Be(0) is the J P = 0 + isospin-singlet ground state of the 8 Be nucleus, and 8 Be * (17.64) and 8 Be * (18.15) are J P = 1 + excited states, whose isospin quantum numbers are predominantly I = 1 and I = 0, respectively, but are nonetheless isospin-mixed: Their level of isospin mixing, quantified by θ 1 + , was estimated by ab initio quantum Monte Carlo techniques [52,92], and by χEFT many-body methods [56] to fall in the approximate range 0.18 sin θ 1 + 0.43. Following [52,56] we will consider a narrower range for sin θ 1 + which more accurately describes the width of the electromagnetic transition 8 Be * (18.15) → 8 Be(0) + γ: In the MTA Atomki experiment [42], a bump-like feature in the m e + e − distribution of the ∆I ≈ 0 transition (3.1b) was observed on top of the monotonically falling spectrum expected from SM internal pair conversion (IPC) [42]. A statistical significance of 6.8 σ was reported for this deviation relative to the IPC expectation. Additionally, ref. [42] claimed that the excess events were consistent with the emission of an on-shell resonance, generically labeled "X", with mass of m X = (16.7 ± 0.35 stat ± 0.5 syst ) MeV, promptly decaying to e + e − . This excess was later corroborated by the same collaboration with a modified experimental set-up [93], with a combined fit yielding a relative branching ratio of: with respect to the radiative γ-width of this transition, 8 Be * (18.15) → 8 Be(0) + γ, of Γ γ (18.15) ≈ (1.9 ± 0.4) eV [94]. As for the m e + e − spectrum of the ∆I ≈ 1 transition (3.1a), no statistically significant deviation from the IPC expectation was observed. Refs. [51,56] inferred a naïve upper bound of for the de-excitation rate of 8 Be * (17.64) via on-shell emission of this hypothetical "X(17)" resonance.
If it is confirmed that the observed excess originates from new, beyond the SM phenomena, as opposed to nuclear physics effects or experimental systematics, it could indeed be explained by the piophobic QCD axion. The prediction for axion emission rates from magnetic dipole nuclear transitions was first worked out by Treiman & Wilczek [57] and independently by Donnelly et al. [58] back in the late 1970's. For the two transitions in (3.1), the axion-to-photon emission rate is (see also [59,61,88]): where | 8 Be * denotes one of the states in (3.1), and its overlap with the isospin eigenstates |I = 0 and |I = 1 follows from (3.2). The quantities µ (0) = µ p + µ n = 0.88 and µ (1) = µ p − µ n = 4.71 are, respectively, the isoscalar and isovector nuclear magnetic moments, and η (0) , η (1) parameterize ratios of nuclear matrix elements of convection and magnetization currents [87]. In particular, η (0) = 1/2 due to total angular momentum conservation. The nuclear structure dependent parameter η (1) , to the best of our knowledge, has not been calculated for 8 Be ; we therefore conservatively vary η (1) in the range: Combining (3.6) with (2.15), (2.16), (3.1b), (3.2b), we can infer the axion isoscalar mixing angles that yield the observed excess rate (3.4). For concreteness, we vary θ aπ within the 1 σ range favored by the KTeV anomaly fit, (2.8), while also varying q e PQ in the range 1/2 ≤ q e PQ ≤ 2, and the nuclear structure parameters θ 1 + and η (1) in the ranges (3.3) and (3.7), respectively. We obtain: (not octe t-en han ced) (not octet -enha nced) Figure 1. Fits, constraints, and sensitivity projections in the parameter space of the axion isoscalar couplings. The upper (lower) plot assumes the same (opposite) relative sign between θ aη ud and θ aη s . The orange and yellow bands enclose the range of isoscalar mixing angles that can explain the 8 Be and 4 He anomalies, respectively, benchmarking ∆u + ∆d and ∆s to the values shown; cf. (3.8, 3.16). The shaded gray regions are excluded by the conservative upper bound Br(K + → π + (a → e + e − )) 10 −5 (under different scenarios for octet-enhancement in χPT) and by current upper bounds on η ( ) → e + e − , assuming q e PQ = 1/2; cf. (4.1a, 4.2a). The dashed gray (red) lines show the expected reach from measurements of (or bounds on) η (η) → e + e − , assuming a future experimental sensitivity to the branching ratios predicted in the SM, (4.1b) and (4.2b). opposite (lower plot) relative sign between θ aη ud and θ aη s . These bands shift non-negligibly as ∆u + ∆d and ∆s are varied within the ranges in (2.17).
Finally, we conclude this discussion by using (3.6) and (3.8) to predict the axion emission rate for transition (3.1a): Indeed, this rate can be down by as much as 2 orders of magnitude below the sensitivity of published results to date, but could potentially be detectable if sufficient statistics is accumulated in this channel.

Evidence for the QCD axion in 4 He transitions
More recently, the same collaboration led by A. Krasznahorkay investigated transitions of a different nucleus, 4 He [43]. With a 900 keV proton beam bombarding a 3 H fixed target, this experiment populated the first two excited states of 4 He: and similarly measured the emission of e + e − pairs from de-excitations of these states to the I = 0, J P = 0 + ground state, denoted here by 4 He(0). Such transitions are allowed via: but forbidden to occur via the following processes: Above, E0 and M0 refer, respectively, to the electric monopole (J P = 0 + ) and magnetic monopole (J P = 0 − ) multipolarities of these transitions. After cuts, background subtraction, and accounting for contributions to the m e + e − spectrum from (3.11a) and from external pair conversion originating from the radiative proton capture reaction 3 H(p, γ) 4 He, a suggestive bump-like excess was observed in the final m e + e − distribution, with a statistical significance of 7.2 σ. Under the assumption that this excess originated from on-shell emission of a narrow resonance from the M0 transition (3.11b), the fit to the data performed in [43] yielded a favored resonance mass of m a = (16.84 ± 0.16 stat ± 0.20 syst ) MeV, and de-excitation width 8 : It is encouraging that not only the same resonance mass (within error bars) is favored by fits to both the 4 He and 8 Be excesses, but also that they appear in magnetic and (dominantly) isoscalar transitions, compatible with the interpretation of piophobic axion emission. To further support this hypothesis, we must obtain the range of axion isoscalar mixing angles compatible with the observed rate. According to Donnelly et al. [58], the width of axionic emission in 0 − → 0 + nuclear transitions is estimated to be 9 : 14) where I N * is the isospin of the excited nuclear state, | p a | ≈ ∆E 2 − m 2 a is the magnitude of the axion's spatial-momentum in the rest frame of the decaying nucleus, Q is a typical nuclear momentum transfer (of order the nucleus Fermi momentum, M0 involve nuclear matrix elements of magnetization currents, and are of O(1), unless forbidden by isospin conservation. For an isoscalar transition such as (3.11b), (3.14) reduces to: Using (3.15), (2.15b), and varying a (0) we find that the axionic de-excitation width of the M0 transition (3.11b) in 4 He yields the observed rate (3.13) if: which is compatible with the range of axion isoscalar mixing angles favored by the 8 Be excess, (3.8). In Fig. 1 we likewise display the parameter space in θ aη ud vs. θ aη s favored by the 4 He anomaly (yellow bands) under the same assumptions for ∆u + ∆d, ∆s and relative sign between θ aη ud and θ aη s used in the computation of the 8 Be orange bands. The 4 He yellow bands also shift non-negligibly as ∆u + ∆d and ∆s are varied within the ranges in (2.17). It is remarkable that the piophobic QCD axion is able to simultaneously explain the reported rate of anomalous excesses in de-excitations of two very different nuclei, 8 Be and 4 He, as shown by the overlap between the fits to the axion isoscalar mixing angles (3.8) and (3.16), or, equivalently, by the overlap between the yellow and orange bands in Fig. 1. This weakens the case for a nuclear physics origin of the observed features in the m e + e − spectra of these transitions [97]. And the fact that "unexplained" features are absent in 9 In [53], the calculated rate for pseudoscalar emission in this 0 − → 0 + transition assumed a nonderivatively-coupled pseudoscalar "X" (see the effective operator in eq. (39) of [53]), resulting in an amplitude with no momentum dependence (eq. (49) of [53]) and an emission rate scaling as ΓX ∝ | pX |. Under this assumption, the authors of [53] concluded that the rate of pseudoscalar emission in this transition would be six orders of magnitude larger than the experimentally favored rate. We point out that their conclusion hinged on their assumption that the leading effective operator at the nuclear level mediating this transition was a relevant operator of dimension-3. In the case of the QCD axion, this assumption is not valid, since the axion only couples derivatively to nuclear axial currents. For the QCD axion, the leading effective nuclear operator is dimension-5, resulting in an amplitude scaling as ∝ | pa| 2 , and therefore an emission rate scaling as Γa ∝ | pa| 5 . Note that the axionic amplitude is still isotropic, as it should be for a monopole transition, despite its nontrivial momentum dependence. For details, see [58,59,96].
the m e + e − spectrum of several other measured transitions-(3.1a) being one example-also makes it less straightforward to "explain away" the observed excesses as poorly understood experimental systematics. We therefore reiterate our point, stated in the Introduction (Sec. 1), that the anomalies in 8 Be and 4 He transitions, and their quantitative compatibility with predicted signals from the QCD axion, should not be quickly dismissed. A cautiously optimistic attitude and support for an independent verification of these measurements are certainly warranted.

η and η decays
In light of the anomalies in nuclear de-excitations discussed in the previous section, a natural next step is to investigate other systems where the hadronic couplings of the piophobic QCD axion could be more precisely determined or more stringently constrained. In this section we consider rare decays of η and η mesons, which, with the prospect of future η/η factories, could become powerful future probes of axions and hadronically coupled ALPs more generally [98]. These include the second phase of the JLab Eta Factory (JEF) program [99], expected to improve existing bounds on rare η decays by two orders of magnitude, and the REDTOP experiment [100,101], a planned η/η factory projected to deliver as many as 10 13 η mesons and 10 11 η mesons. These will offer an unprecedented opportunity to study rare η/η decays and probe BSM physics.

Di-electronic η and η decays
Just as the precise KTeV measurement of π 0 → e + e − offered the best determination of θ aπ , future observations of η → e + e − and η → e + e − could narrow down the ranges for the axion isoscalar mixing angles θ aη ud and θ aη s . Present bounds on these dileptonic branching ratios [102,103] are still two orders of magnitude away from sensitivity to the predicted SM rate [104,105]: and Indeed, the highly suppressed SM contribution to these dileptonic channels makes them potentially sensitive to a variety of interesting new physics scenarios. In anticipation of a future discovery of these decay modes, we obtain the axionic contribution to the dileptonic decays η ( ) → e + e − due to a−η ( ) mixing. Assuming that this effect dominates these rates (i.e., that interference with the SM amplitudes can be neglected), we have: The mixing angles θ aη and θ aη can be re-expressed in terms of θ aη ud and θ aη s using the parameterization [106] |η = cos φ ud |η ud − sin φ s |η s , (4.4a) from which it follows that (4.1a) and (4.2a) translate into relatively weak bounds on the axion isoscalar mixing angles: (4.5b) Taking φ ud = 39.8 • and φ s = 41.2 • from [106] and conservatively assuming |q e PQ | = 1/2 for concreteness, we display the bounds (4.5) in Fig. 1, as well as the potential reach in θ aη ud and θ aη s assuming that future η/η -factories will achieve sensitivity to the branching ratios predicted by the SM, (4.1b) and (4.2b).

Axio-hadronic η and η decays
Hadronic decay channels of η and η mesons could in principle be hiding promising signals of the QCD axion and/or other hadronically coupled ALPs. Amongst the most obvious modes are the three-body final states η ( ) → π 0 π 0 a , π + π − a, which have only recently been explored in the literature [46,107]. Indeed, the amplitudes for these processes receive a direct contribution from the leading order potential term in the chiral Lagrangian 10 , and could in principle result in considerably large branching ratios. The difficulty with studying hadronic η and η decays lies in reliably predicting their rates. One of the earliest examples where this difficulty was encountered was in the calculation of η → 3π, which was significantly underestimated by χPT at leading order [108][109][110][111]. Indeed, it has long been understood that contributions from chiral logarithms and strong final state rescattering could not be neglected in the computation of η → 3π [112][113][114][115][116][117][118]. Similarly, neither the total width nor the Dalitz phase space of η → ηππ are properly described by χPT at O(p 2 ) [119,120].
Previous studies [121,122] have shown that such intermediate energy-processes can be satisfactorily described by extending χPT to include low-lying meson resonances carrying non-linear realizations of SU (3) χ -such as vectors (ρ, ω, K * , φ, ...), axial-vectors (a 1 , f 1 , K 1 , ...), scalars (a 0 , f 0 , σ, κ, ...), and pseudo-scalars (η , π(1300), ...)-and assuming the principle of "resonance dominance" (an extension of vector meson dominance), whereby the low energy constants (LECs) of the O(p 4 ) chiral Lagrangian are saturated by mesonic resonance exchange. This framework, dubbed Resonance Chiral Theory (RχT), has been quite successful phenomenologically as an interpolating effective theory between the short-distance QCD description and the low-energy χPT framework, by encoding the most prominent features of nonperturbative strong dynamics [123,124]. There does not appear to be consensus in the literature, however, on which low-lying resonances should be included as degrees-of-freedom in the RχT Lagrangian, and which resonances should be regarded as dynamically generated poles due to strong S-wave interactions [125,126]. Examples of such "ambiguous" poles include the σ(500) and the κ(700).
In this subsection, we estimate the rates for η → ππa and η → ππa using RχT. In both cases, we find that the leading order χPT predictions for these decay rates are significantly modified by inclusion of resonance exchange amplitudes. In particular, for η → ππa, there is substantial destructive interference between the leading order amplitude from the O(p 2 ) quartic term and the amplitudes generated by tree-level resonance exchange. This is corroborated by performing the same calculation in ordinary χPT at O(p 4 ), where one finds that the LECs, in particular L 4 , L 5 , and L 6 , provide O(1) contributions that destructively interfere with the O(p 2 ) amplitude. For η → ππa, the contributions from resonance exchange (alternatively, from χPT interactions at O(p 4 )) are the dominant effect in a significant portion of the parameter space, and may enhance this decay rate by an order of magnitude over the leading order χPT prediction. The justification for favoring the RχT framework over O(p 4 )-χPT for this calculation is that the former is expected to better capture the Dalitz phase-space of the final state, which is relevant when extracting the event acceptance due to momentum cuts in experimental analyses (in particular due the e ± selection criteria). Indeed, we will find that there is strong variation of the amplitude's momentum dependence as we vary the assumptions and parameters of the RχT description, which implies a strong variation in the estimated sensitivity of existing and future experimental analyses. Unfortunately, the variations in these assumptions cannot be narrowed down without further input from experiment. Our main conclusion, therefore, is that one cannot reliably predict neither the total branching ratio, nor the Dalitz phase-space, of the decays η ( ) → ππa. Under reasonable assumptions, our RχT-based estimates vary over two orders of magnitude in branching ratio, Br(η ( ) → ππa) ∼ O(10 −4 − 10 −2 ). Nonetheless, this motivates dedicated reanalyses of existing data in final states of η ( ) → π π e + e − , as well as dedicated searches for e + e − resonances in these final states in future η/η factories.
In order to motivate our use of RχT, and also to justify our later approximation of retaining only low-lying scalar resonances, we begin by obtaining the main contributions to the amplitude A(η ( ) → π 0 π 0 a) in ordinary χPT at O(p 4 ). The Lagrangian is: Above, f π = 92 MeV, B 0 = − qq /f 2 π , M 0 parameterizes the O(GeV) contribution to the mass of the chiral singlet η 0 from the strong axial anomaly, L i (i=1,...,10) are the O(p 4 ) χPT low energy constants (LECs), M q (a) is the axion-dependent quark mass matrix, transforming as an octet spurion of SU (3) χ : and U is the non-linear representation of the pseudo-Goldstone chiral nonet: Collecting the terms in (4.6) that provide the dominant contributions to A(η ( ) → π 0 π 0 a), we obtain: 11) and the kinetic terms, omitted in (4.9), have been canonically normalized. While there is large variation in the literature of the inferred values for L 4 and L 6 from fits to experimental data, depending on assumptions and chosen observables, it is well established from fits to f K /f π that L 5 is positive and relatively large, It is then easy to see from (4.9) and (4.11) that the contributions to A(η → π 0 π 0 a) from the first and second terms in (4.9) are comparable in magnitude and destructively interfere with each other. This leads to a suppressed rate for η → π 0 π 0 a relative to the naïve O(p 2 ) estimation in χPT, which, however, is quite sensitive to the value of L 5 . On the other hand, the O(p 2 ) contribution to A(η → π 0 π 0 a) from the first term in (4.9) may be subdominant to that of the second term, which is parametrically larger by a factor of O( L 5 m 2 η /m 2 η ). This may lead to an order-of-magnitude enhancement of the rate for η → π 0 π 0 a.
While it is now straightforward to extract χPT's prediction for Br(η ( ) → π 0 π 0 a) using (4.9), we will instead pivot to RχT, from which ordinary χPT can be recovered by integrating out the low-lying meson resonances. Under the assumption of resonance dominance, RχT predicts that the relevant LECs contributing to η ( ) → ππa (L 4 , L 5 , and L 6 ) are saturated by the exchange of scalar resonances. We will therefore omit the lowlying pseudoscalar, vector, and axial-vector resonances from our discussion. Following the notation in [123], we have: where S is the low-lying J P C = 0 ++ meson octet 11 , (4.13) Above, a 0 and f 0 are shorthand for a 0 (980) and f 0 (980), respectively [127], and we have not explicitly identified the scalar mesons with non-zero strangeness, since they do not contribute to η ( ) → ππa.
Accounting for the tadpole-induced non-zero vacuum expectation value of f 0 , (4.14) and canonically normalizing the kinetic terms, we can extract from (4.12) the RχT interactions contributing to η ( ) → π 0 π 0 a: where, following the notation for the dimensionless fields and derivatives in (4.9), we have additionally introduced: It is straightforward to recover (4.9) from (4.15) by integrating out the scalar resonances and making the following identifications: Early fits to a 0 (980) → πη [123], along with large-N c assumptions and imposition of short-distance constraints [128] (such as sum-rules between two-point correlators of two scalar vs two pseudoscalar currents [129], and vanishing of scalar form factors at q 2 → ∞ [130,131]) have been used to estimate the scalar octet couplings to be |c d | ∼ |c m | ∼ f π /2, with c d c m > 0. Here, we conservatively vary these values by ± 20% in our calculations of Γ(η ( ) → ππa), but retain, for simplicity, the assumption of |c d | = |c m |.
With the relevant interactions in (4.15), we can then obtain the tree-level amplitude 12 A(η ( ) → ππa) (see Fig. 2): 12 We ignore corrections to A(η ( ) → ππa) from ππ final state rescattering, based on the conclusions from [112,114] that these effects correct the η → 3π amplitude by modest amounts of O(10%), and on reference [113], which finds somewhat larger rescattering corrections, of ∼ 70%, which are still subdominant relative to other sources of uncertainties in our estimations.  where we have neglected subdominant contributions of O(m 2 π /m 2 η ( ) ). Above, the equality between amplitudes with neutral vs charged pions is due to isospin symmetry; p η ( ) , p π 1 , and p π 2 are relativistic 4-momenta; and For the η/η mixing angles and decay constants above, we will adopt the values from the unconstrained fit in [106], namely, θ 8 = −24 • , θ 0 = −2.5 • , f 8 = 1.51 f π , and f 0 = 1.29 f π . Finally, we can obtain the differential decay rate from (4.20): where S π 1 π 2 is the standard combinatorial factor (S π + π − = 1, S π 0 π 0 = 2! ), dΦ 3 is the 3-body phase-space differential element, and, when obtaining the total decay rate, the integration over invariant masses m 2 π 1 π 2 = (p π 1 + p π 2 ) 2 and m 2 π 2 a = (p π 2 + p a ) 2 = (p η ( ) − p π 1 ) 2 should be performed over the Dalitz phase-space (see, e.g., the PDF review on Kinematics [127] for explicit expressions for the Dalitz plot boundaries). In Fig. 3, we show the branching ratios for η ( ) → π 0 π 0 a, η ( ) → π + π − a (computed by integrating the differential rate in (4.22) over the final state phase space) as a function of the RχT couplings c d , c m of the low-lying scalar octet resonances to the pseudo-Goldstone mesons; see eqs. (4.12, 4.15, 4.17). As mentioned previously, we assume, for simplicity, that c d = c m , and vary their magnitudes by ±20% around their expected values of | c d | = | c m | = 1 in the large-N c limit. The range covered by the bands are due to the uncertainties in the masses and widths of the scalar resonances a 0 (980) and f 0 (980)-following the PDG [127], we varied these parameters independently within the following ranges: The lack of predictive power of our treatment, with an estimated range of branching ratios spanning two orders of magnitude, Br(η ( ) → ππa) ∼ O(10 −4 − 10 −2 ), is due to the χPT and RχT parameters falling on a special range of values that, within uncertainties, can lead to substantial destructive inference between the LO amplitude and the amplitudes originating from exchange of low-lying scalar resonances. This is perhaps unsurprising, considering that even the SM hadronic decays of the η and η could not be correctly predicted, but only "postdicted," and their experimentally determined branching ratios and Dalitz plot parameters have been used to verify the validity of various treatments and assumptions, such as RχT, QCD sum rules, large-N c limit, dispersive methods, etc [112][113][114][115][116][117][118][119][120][121][122].
In particular, the upper range of our estimations, Br(η ( ) → ππa) ∼ O(10 −2 ), is probably excluded or in tension with observations, though no dedicated searches for an e + e − resonance in η ( ) → π π e + e − final states have ever been performed, to the best of our knowledge. However, the lower range Br(η ( ) → ππa) ∼ O(10 −4 − 10 −3 ) likely remains experimentally allowed, and within the sensitivity of upcoming η/η factories, such as the JLab Eta Factory (JEF) and the REDTOP experiment.
The most recent and precise measurement of the SM decay η → π + π − (γ * → e + e − ), which shares the same final state of η → π + π − a, was performed by the KLOE collaboration at the Frascati φ-factory DAΦNE [132]. While their measurement yielded Br(η → π + π − e + e − ) = (2.68 ± 0.09 stat ± 0.07 syst ) × 10 −4 , it is non-trivial to infer any bounds from this analysis on Br(η → π + π − a). This is because, without proper Monte Carlo simulations, one cannot determine how the background rejection requirements would have affected the η → π + π − a signal efficiency. In particular, this search rejected events with m e + e − < 15 MeV whose reconstructed e + e − vertex was within a 2.5 cm distance from the beampipe. This cut could have significantly impacted the acceptance of the axion signal, depending on the m e + e − experimental resolution. Other event selection requirements on the momenta of the π ± and e ± charged tracks could in principle have rejected a large fraction of the axion signal as well.
An earlier measurement of Br(η → π + π − (γ * → e + e − )) by the CELSIUS/WASA collaboration observed, in hindsight, an upward fluctuation of the expected signal [133,134]. Indeed, considering KLOE's more precise measurement of this branching ratio, the CEL-SIUS/WASA analysis should have expected 10 SM signal events. It observed 24 events in the signal region, of which it determined that 7.7 were from background, and 16.3 were  Figure 4. The differential rate for η → π + π − a as a function of | p e + e − | ≡ | p e + + p e − | = p a , for three benchmark choices of RχT parameters specified in Table 1. For comparison, we also show the differential rate of the SM process η → π + π − e + e − , labeled "QED".
from the SM signal. Assuming, conservatively, that the 14 "excess events" were instead due to η → π + π − a decays, and taking into account the relative signal acceptance due to the minimum transverse momentum requirement of | p T | > 20 MeV for charged particles 13 , we estimate that branching ratios as large as Br(η → π + π − a) ∼ (1 − 3) × 10 −3 could be compatible with the CELSIUS/WASA measurement, although, without access to non-public information on details of the experimental analysis, this estimation is at the level of an educated guess.
Finally, the two existing measurements of Br(η → π + π − (γ * → e + e − )), performed independently by the CLEO [135] and BESIII [136] collaborations, were combined by the PDG [127] to give Br(η → π + π − e + e − ) = 2.4 +1.3 −0.9 × 10 −3 . However, both experimental analyses reported large external photon-conversion backgrounds in the signal region, peaked in the range m e + e − = (8 − 25) MeV (CLEO; see Fig. 2(d) of [135]) and m e + e − = (10 − 20) MeV (BESIII; see Fig. 2 of [136]). Events falling within these m e + e − windows were excluded from the analyses' inference of the SM branching ratio. Since events from η → π + π − a would have fallen precisely in this region where the photon conversion background peaked, it is difficult to estimate how strong a potential axion signal could have been. 13 We performed a simple MC event simulation to estimate the geometric acceptance resulting from the event selection requirement of | p e ± T | > 20 MeV, properly taking the momentum dependence of the amplitudes into account. This was done for both the SM signal using the amplitude in [133], as well as for the axion signal, assuming a few RχT benchmark parameters (see discussion below). We neglected contributions to the signal efficiency from other event selection requirements, and worked in the approximation of η mesons decaying at rest.
Simply requiring that the axion signal strength does not overpredict the number of events attributed to photon-conversion yields a conservative limit of Br(η → π + π − a) few×10 −2 , which is not particularly useful.
We end this section by remarking that an additional challenge with estimating the sensitivity of current and future experiments to η ( ) → ππa is the uncertainty in the final state Dalitz phase-space, which affects the signal acceptance resulting from event selection cuts. Consider, for instance, the differential decay rate dΓ(η → π + π − a)/d| p a | as a function of the axion's 3-momentum | p a |. The dependence of this rate on | p a | = | p e + + p e − | ≡ | p e + e − | varies dramatically depending of the numerical values chosen for the masses, widths, and couplings of the scalar resonances a 0 and f 0 . We illustrate this effect in Fig. 4, where we plot the differential decay rate dΓ(η → π + π − e + e − )/d| p e + e − | as a function of | p e + e − | for three different RχT benchmarks-corresponding to different choices of masses and widths for a 0 and f 0 within uncertainties (see Table 1)-as well as for the SM decay η → π + π − (γ * → e + e − ) [104] (labeled as "QED" in Fig. 4). While these different benchmark points yield close predictions for Br(η → π + π − a), their predictions for dΓ(η → π + π − a)/d| p a | differ dramatically, as shown in Fig. 4. In particular, for high enough cuts on the charged lepton momenta p e ± , the signal acceptance of benchmark B1 could be significantly lower than that of B3 (and of the SM decay η → π + π − γ * ). Indeed, this could lead to a variation of as much as an order of magnitude in the expected sensitivity of experimental searches.  Fig. 4, and the resulting prediction for the total decay rate of η → π + π − a .

Kaon decays
We conclude our study by exploring signals of the piophobic QCD axion in rare Kaon decays. Although the main focus of ongoing and near-future rare Kaon decay experimentssuch as NA62 at CERN [137] and KOTO at J-PARC [138]-has been on K → πνν, there is an under explored opportunity to search for BSM resonances in e + e − final states with low m e + e − , motivated not only by the piophobic QCD axion, but also by visibly decaying ALPs and dark-photons more generally [41]. Furthermore, the highly suppressed a → γγ decay mode might be a competitive final state in Kaon decay searches for which γγ backgrounds in the m γγ ∼ 17 MeV signal region are tamer than the e + e − backgrounds. In such cases, final states with K → (a → γγ) + SM can be obtained by combining the relevant branching ratios Br(K → a + SM) estimated in this section with Br(a → γγ) in (2.13).
The appearance of the axion in Kaon decay final states occurs via mixing with the neutral octet mesons, π 0 , η, and η . Therefore, the axionic amplitudes can be obtained from ordinary SM amplitudes properly reweighted by axion-meson mixing angles. While this prescription is straightforward for estimating "axio-leptonic" Kaon decays such as K + → µ + ν µ a (as we will show in Subsec. 5.1), it is ambiguous for "axio-hadronic" Kaon decays such as K + → π + a, K 0 S, L → π 0 a, and K 0 L → ππa. Firstly, the two-body hadronic width of the CP-even neutral Kaon, Γ(K 0 S → π 0 π 0 , π + π − ) ≈ 0.73×10 −5 eV, is enhanced by roughly three orders of magnitude relative to the two-body hadronic width of the charged Kaon, Γ(K + → π + π 0 ) ≈ 1.1 × 10 −8 eV. In χPT, this enhancement is parametrized as a large disparity in the magnitudes of the Wilson coefficients of the possible ∆S = 1 operators [108,[139][140][141]. Specifically, the coefficient of an SU (3) χ -octet (∆I = 1/2) operator is larger than the coefficient of the leading order 27-plet (∆I = 3/2) operator by a factor of ∼ 30. Secondly, phenomenologically, there are at least two choices of ∆S = 1 octet operators that could be responsible for this so-called "octet enhancement" (a.k.a. "∆I = 1/2 enhancement") in Kaon decays [142][143][144], namely, where λ ds ≡ (λ 6 +iλ 7 )/2, Λ ∼ 4πf π is a natural χPT cut-off, and the operators above occur at different orders in the chiral expansion: O 8 at O(p 2 ), and O 8 at O(p 4 ). Fitting existing data on K → ππ and K → πππ, while treating the coefficients g 8 and g 8 in (5.1) on equal footing, yields |g 8 + g 8 | 0.78 × 10 −7 [145,146]. In order to break this degeneracy in the fit, one must invoke the standard assumption under naïve power counting that ∆S = 1 octet enhancement should appear at lowest order in the chiral expansion, and therefore, g 8 g 8 ⇒ |g 8 | 0.78 × 10 −7 . However, there is no first principles derivation of this choice, and it could be incorrect. For example, in the Resonance Chiral Theory framework discussed in the previous section, it is easy to speculate that the origin of ∆S = 1 octet enhancement could be due to the weak interactions inducing a mixing between K 0 S and a broad J P C = 0 ++ resonance, such as the σ(500) 14 . Upon integration of the low-lying resonances, this effect would be captured by the operator O 8 in (5.1b), leading instead to This ambiguity directly affects predictions for rare Kaon decays to the piophobic axion, since the ∆S = 1 octet operators O 8 and O 8 contribute differently to the amplitudes A(K + → π + a), A(K 0 S, L → π 0 a), and A(K 0 L → ππa). In what follows, we will estimate the rates for various axionic Kaon decays in both scenarios, g 8 g 8 and g 8 g 8 . We will show that in the case of g 8 g 8 , all amplitudes A(K + → π + a), A(K 0 S, L → π 0 a), and A(K 0 L → ππa) are octet-enhanced, leading to higher axionic Kaon decay rates, and when relevant, more stringent constraints on the mixing angles θ aη ud and θ aη s . Conversely, in the scenario g 8 g 8 , the rates Γ(K + → π + a) and Γ(K 0 S, L → π 0 a) are significantly reduced, relaxing the otherwise strong constraints on θ aη ud and offering an exciting prospect for searching for these signals in near-future rare Kaon decay experiments.
14 Indeed, a naïve dimensional analysis estimation of g 8 ∼ |GF sin θc f 2 π Λ 2 /m 2 0 ++ | ∼ O(10 −7 ) does not immediately rule out this hypothesis. It is unclear whether a Dalitz plot analysis of K 0 L → 3 π data could distinguish it from the alternative description of octet-enhancement.
In upcoming subsections, we will normalize the calculated axio-hadronic rates to analogous Kaon decay rates in the SM. For later reference, we quote here the dependence of the relevant SM Kaon decay amplitudes on g 8 , g 8 , as well as g 27 , the coefficient of the 27-plet ∆S = 1 operator at O(p 2 ) in χPT, which is given by: Above, T ij kl are Clebsch-Gordan coefficients that project the 27-plet, ∆I = 3/2 part of the interaction [140]. The contributions of (5.1a), (5.1b), and (5.2) to the two-body hadronic Kaon decays of interest are [145,147]: As alluded to earlier, the hadronic K + decay amplitude is not octet-enhanced, and its consequent narrow width relative to that of K 0 S is parametrized by a hierarchy between the 27-plet and octet coefficients: Finally, for the relevant hadronic 3-body decay of K 0 L , we have [145,147]: where Y is one of the standard Dalitz plot variables, defined as: with p 1 and p 2 referring to the four-momenta of the charged pions, and p 3 the fourmomentum of the neutral daughter particle 15 , in this case π 0 .

Axio-leptonic K + decays
The amplitude for the axio-leptonic decay K + → + ν a can be easily related to the SM semi-leptonic amplitudes via the Ademollo-Gatto theorem [148], which states that the matrix elements of flavor-changing electroweak current operators can only deviate from their SU (3) χ -symmetric values to second order in chiral symmetry breaking [149,150]. This implies that, at zero momentum transfer, the following SU (3) χ relations hold:   where is a measure of SU (3) χ breaking. Then, since |a = θ aπ |π 0 + θ aη 8 |η 8 + θ aη 0 |η 0 , we have: Neglecting the difference in phase space, as well as finite momentum-transfer and SU (3) χ breaking corrections, which amount to O(10%) [147], it then follows from (5.9) that: In the specific case of a muonic final state, (5.10) yields: (5.11) In Fig. 5, we show the hypothetical reach in axion isoscalar mixing angles from K + → µ + ν µ a, assuming an experimental sensitivity to branching ratios Br(K + → µ + ν µ a) 10 −8 . Note that this branching ratio sensitivity figure has been chosen to facilitate comparison between different axionic Kaon decay modes (to be discussed in upcoming subsections), and is not informed by any experimental sensitivity projections.

Axio-hadronic K + decays
For axio-hadronic Kaon decays, we should first obtain the contributions to the amplitudes for K + → π + ϕ * (ϕ = π 0 , η ud , η s ) from operators (5.1a), (5.1b), and (5.2). Putting K + and π + on-shell, we have: 12c) The axionic decay K + → π + a is then induced by these amplitudes via axion-meson mixing: Note that (5.13) depends on g 8 but not on g 8 . This implies that A(K + → π + a) is only octet-enhanced in the scenario with g 8 g 8 , i.e., in the standard realization of octetenhancement in χPT via O 8 . In this case, using (5.3) and taking g 8 → 0 for simplicity, we can approximate (5.13) as: where K ππ ∼ 3 corrects for the fact that strong s-wave ππ final state interaction, present in K 0 S → π + π − , is absent in K + → π + a [60]. With (5.14) and (5.5) we finally obtain: In the scenario g 8 g 8 , A(K + → π + a) is not octet-enhanced, and, since g 8 is expected to be of the same magnitude as g 27 , there might be non-negligible interference between the contributions stemming from O 8 and O 27 . For simplicity, we will ignore this effect and consider the limiting case g 8 → 0, when O 27 provides the dominant contribution to axiohadronic K + decays. In this case, using (5.4), (5.13) can be approximated as: from which it follows that: (5.17) In Figs. 1 and 5, we show the excluded parameter space for the axion's isoscalar mixing angles-assuming a conservative experimental bound of Br(K + → π + a) 10 −5 (see discussion in [45])-for the two assumed scenarios for octet-enhancement in χPT which resulted in (5.15) and (5.17).
We also show in Fig. 5 the hypothetical reach in axion isoscalar mixing angles from K + → π + a, assuming an experimental sensitivity 16 to branching ratios Br(K + → π + a) 10 −8 . The K + → π + a contours in Fig. 5 make evident that this axionic Kaon decay channel is one of the most competitive in probing the piophobic QCD axion, and that updating the three-decades-old bounds on K + → π + e + e − with m e + e − 50 MeV [151][152][153] could cover presently unexplored and well motivated parameter space of light BSM sectors.

K 0
S decays The axio-hadronic decays of the CP-even neutral Kaon can be estimated via an analogous prescription as the one used in Subsec. 5.1.2. First, we obtain the contributions to A(K 0 S → π 0 ϕ * ), ϕ = π 0 , η ud , η s , from operators (5.1a), (5.1b), and (5.2). With K 0 S and π 0 on-shell, we have: The axionic decay K 0 S → π 0 a is then induced by these amplitudes via axion-meson mixing: Note that the only occurrence of g 8 in (5.19) stems from the axion-pion mixing contribution in (5.18a), and, as such, it is suppressed by the small θ aπ mixing angle. Because of this, A(K 0 S → π 0 a) parallels the behavior of A(K + → π + a) of only being octet-enhanced in the scenario with g 8 g 8 . In this case, using (5.4), we can approximate (5.19) as: 20) to then obtain: In the alternative scenario of g 8 g 8 , when g 8 and g 27 are expected to have comparable magnitudes, there will be non-negligible interference between the O 8 and O 27 contributions to the amplitudes (5.18b) and (5.18c). For simplicity we again consider the limiting case g 8 → 0 to arrive at the following approximation: from which it follows that: Note that in (5.23) the contribution from axion-pion mixing is non-negligible despite the suppression of θ aπ relative to θ aη s . As alluded to earlier, this is because (5.18a) is the only octet-enhanced amplitude contributing to K 0 S → π 0 a in the scenario with g 8 g 8 . In Fig. 5, we display the hypothetical reach in axion isoscalar mixing angles from K 0 S → π 0 a, assuming an experimental sensitivity 17 to branching ratios Br(K 0 S → π 0 a) 10 −8 , under both octet-enhancement scenarios in χPT. Since there is a non-negligible contribution from θ aπ in the scenario with g 8 g 8 , we have chosen the relative sign between θ aπ and θ aη s that yields the most conservative reach in the parameter space of Fig. 5 (cf. (5.23)).

K 0
L decays 5.3.1 CP-violating axio-hadronic K 0 L decays Direct CP-violation in the neutral Kaon system causes K 0 L to inherit the axio-hadronic decay modes of K 0 S . The resulting branching ratios can be trivially obtained by accounting for K 0 L − K 0 S mixing, parametrized by the parameter K 2.23 × 10 −3 : In particular, for the two octet-enhancement scenarios in χPT considered in Subsec. 5.2, we have: 5.3.2 CP-conserving axio-hadronic K 0 L decays The CP-conserving axio-hadronic decays of the CP-odd neutral Kaon can be estimated via an analogous prescription as the one used in Subsecs. 5.1.2 and 5.2. For specificity, we consider the final state with charged pions, and obtain the contributions to A(K 0 L → π + π − ϕ * ), ϕ = π 0 , η ud , η s , from operators (5.1a), (5.1b), and (5.2). Putting K 0 L , π + , and π − on-shell, we have: where the Dalitz plot variable Y has been defined in (5.7). Furthermore, an additional effect contributing to A(K 0 L → π π a) must be taken into account, namely, kinetic mixing between K 0 L and the neutral pseudoscalar mesons, induced by operators O 8 and O 27 : In particular, accounting for K 0 L − π 0 mixing is crucial in order to obtain the correct dependence of the SM amplitude A(K 0 L → π + π − π 0 ) on g 8 and g 27 , see (5.6). Similarly, the contribution to the axionic amplitude A(K 0 L → π π a) stemming from K 0 L − η ( ) mixing becomes important in the scenario with g 8 g 8 when |θ aη ud |, |θ aη s | O(10 −4 ). It can be straightforwardly obtained by re-weighting A η ( ) →ππa in (4.20): and we have used (4.4) and (4.21) in obtaining the approximate equality in (5.30). Unfortunately, the amplitude in (5.29) is subject to the same destructive interference effects, and therefore the same large uncertainties, as A(η ( ) → π π a) estimated in Subsec. 4.2. In particular, for the region of parameter space where (5.29) becomes important, these uncertainties make the estimation of Br(K 0 L → π π a) unreliable. Nonetheless, for the sake of illustration, we will include the effects of K 0 L − η ( ) mixing in our estimations below by benchmarking the RχT parameters entering in (5.29) The total amplitude A(K 0 L → π + π − a) is then given by the sum of (5.29) and (5.27a), (5.27b), and (5.27c) reweighted by axion-meson mixing: Note that A(K 0 L → π + π − a) is octet-enhanced in both scenarios, and therefore we can neglect the contributions from g 8 and g 27 when g 8 g 8 , and likewise neglect the contributions from g 8 and g 27 when g 8 g 8 . Despite this simplification, obtaining the dependence of Br(K 0 L → π + π − a) on the axion-meson mixing angles still involves non-trivial integration of the differential decay width over the three-body final state phase space. We performed this integration numerically for both octet enhancement scenarios under the assumptions stated above, and using (5.6) for normalization, obtained: aη s + θ aπ − 4.125 θ aη ud + 3.51 θ aη s − 6.14 θ aη ud θ aη s , and Br(K 0 L → π + π − a) g 8 g 8 ≈ 1.49 θ aπ + θ aη ud + √ 2 θ aη s 2 .
It is worth mentioning an exception to our estimations of A(K 0 L → π + π − a) presented in this section. Besides the two limiting cases we have been considering of g 8 g 8 and g 8 g 8 , there is also the possibility that g 8 and g 8 have comparable magnitudes. In this case, there would be non-negligible interference between the amplitudes for K 0 L → π + π − a originating from operators O 8 and O 8 , which could substantially modify the dependence of Br(K 0 L → π + π − a) on the axion-meson mixing angles. Despite the assumptions, simplifications, and uncertainties of our estimations, the K 0 L → π + π − a contours in Fig. 5 offer a compelling motivation for upcoming Kaon experiments to search for an m e + e − ∼ 17 MeV resonance in K 0 L → π π e + e − final states 19 . If these searches could achieve sensitivities down to branching ratios of O(10 −8 ), they could almost fully exclude the parameter space favored by (or verify the QCD axion explanation of) the 8 Be, 4 He, and KTeV anomalies.

Di-electronic K 0
L decays The last rare Kaon decay we shall consider is K 0 L → e + e − , whose amplitude can receive a potentially non-negligible contribution from K 0 L − η ( ) − a mixing. Indeed, using (5.28) and momentarily neglecting the SM contribution to the amplitude, it follows that the K 0 L → e + e − rate induced by K 0 L − η ( ) − a mixing would be: The estimate above should be contrasted with the observed K 0 L → e + e − rate [157], as well as the range of SM predictions [158][159][160] 18 This choice of branching ratio sensitivity benchmark of 10 −8 is intended to facilitate comparison between different axionic Kaon decay modes, and is not informed by any experimental sensitivity projections. 19 Dedicated reanalyses of existing data in K 0 L → π π e + e − final states could also be potentially sensitive to the axionic signal. See, e.g., [154][155][156].

(5.38)
For | q e PQ (θ aη ud + √ 2 θ aη s )| 10 −3 , however, the SM contribution to the amplitude cannot be neglected; in this case, the estimation in (5.35) is not accurate. Indeed, for a significant part of the parameter space favored by the 8 Be and 4 He anomalies, the axionic contribution to K 0 L → e + e − is subdominant to that of the SM, and in fact it is negligible in the octet-enhancement scenario with g 8 g 8 . A tantalizing possibility remains, however, that once measurements and theoretical predictions are improved over (5.36) and (5.37), a piophobic QCD axion signal will appear in this channel as an excess over the SM expectation.

Summary and Discussion
The PQ mechanism was conceived to address a problem intrinsic to the nonperturbative dynamics of QCD. Yet, presently, the prevalent view is that PQ symmetry breaking should take place at scales f PQ 10 10 Λ QCD . Why should these two scales be so widely separated? PQ cancellation of the strong CP phase would be much more robust against spoiling CP violation effects if f PQ ∼ Λ QCD . This possibility has long been dismissed due to stringent laboratory constraints on the visible QCD axion, in particular on its isovector couplings. More recently, however, we showed [45] that the O(10) MeV mass range for the QCD axion remains compatible with all experimental constraints if the QCD axion (i) couples dominantly to the first generation of SM fermions; (ii) is short-lived (decaying with lifetimes 4 × 10 −14 s to e + e − ); and (iii) is piophobic, i.e., has suppressed isovector couplings due to an accidental cancelation of its mixing with the neutral pion, θ aπ 10 −4 . These conditions require non-trivial UV completions, but so does any viable QCD axion model, whether "heavy" and short-lived or ultralight and cosmologically long-lived.
While this possibility forgoes the attractive feature of explaining the particle nature of dark matter, it offers a single, consistent explanation for a few persistent experimental anomalies: the observed rate for π 0 → e + e − , and the "bump-like" excesses in the e + e − spectra of (predominantly) isoscalar magnetic transitions of excited 8 Be and 4 He nuclei. Unsurprisingly, such signals have long been predicted as "smoking-gun" signatures of the QCD axion.
In this article we estimated the axionic emission rates of the relevant 8 Be and 4 He transitions, taking nuclear and χPT uncertainties into account, and showed that the piophobic QCD axion provides a natural and compelling explanation of the observed data for these two nuclei with quite distinct properties. We also considered in detail potential axionic signals in rare decays of the η, η , K ± , and K 0 S,L mesons. The (often ignored) hadronic and χPT uncertainties involved in estimations of these rare meson decays impede accurate predictions of axio-hadronic signals; nonetheless, the ranges we have obtained for several of the processes investigated can be probed in the near future by a variety of experimental programs, including η/η and Kaon factories.