A viable QCD axion in the MeV mass range

The QCD axion is one of the most compelling solutions of the strong CP problem. There are major current efforts into searching for an ultralight, invisible axion, which is believed to be the only phenomenologically viable realization of the QCD axion. Visible axions with decay constants at or below the electroweak scale are believed to have been long excluded by laboratory searches. Considering the significance of the axion solution to the strong CP problem, we revisit experimental constraints on QCD axions in the O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(10 MeV) mass window. In particular, we find a variant axion model that remains compatible with existing constraints. This model predicts new states at the GeV scale coupled hadronically, and a variety of low-energy axion signatures, such as rare meson decays, nuclear de-excitations via axion emission, and production in e+e− annihilation and fixed target experiments. This reopens the possibility of solving the strong CP problem at the GeV scale.


Introduction
The Peccei-Quinn (PQ) mechanism is arguably the most compelling solution of the strong CP problem. Soon after it was originally proposed [1,2], it was realized that a light pseudoscalar would emerge in the infrared spectrum as a manifestation of the underlying PQ mechanism -the QCD axion [3,4]. As the parameters characterizing the QCD axion,
Given the significance of the QCD axion and the many resources dedicated to probing its existence, it is important to ensure that no caveats have been overlooked, and that no gaps in already probed regions of parameter space were left uncovered. Motivated by this, we revisit constraints on the MeV mass window for the QCD axion and its variants, and discuss a particular realization of the QCD axion that has not yet been definitively excluded. Our results reopen the possibility that the strong CP problem might be solved below the weak scale, and suggest new, hadronically-coupled degrees of freedom at the GeV scale.
In sections 2-5 of this paper, we refute previous, premature conclusions that the MeV mass range for the QCD axion has been completely ruled out. We then discuss a viable axion variant and its couplings to photons, nucleons, and electrons, and the relevant experimental implications (sections [6][7][8]. Finally, we comment on UV completions of such variants and associated phenomenology (sections 9-10). For an outline of this paper, vide table of contents.

Constraints on generic MeV axions
Generic constraints on QCD axions in the MeV mass window can be broadly classified into three categories: (i) amenable to model-building, (ii) plagued by large hadronic uncertainties, and (iii) evaded only by pion-phobia. The first category of constraints can be evaded by well-established model-building tools. Constraints in the second category are more difficult to avoid with model-building, but suffer from significant uncertainties which preclude them from fully and unambiguously ruling out the MeV mass range for the QCD axion. Finally, the third category encompasses the strongest constraints, which can only be avoided by a special class of axion variants which are pion-phobic, i.e., which have suppressed mixing with the neutral pion. While this can be achieved by model-building to some degree, the extreme pion-phobia needed to avoid exclusion also depends critically on the light quark mass ratio being close to a ratio of Peccei-Quinn charges (this will become 1 For fa 100 GeV, the QCD axion is also constrained by stellar evolution [13][14][15], and most recently by a combination of the CMB power spectrum and primordial 4 He and D/H abundances [16].

JHEP07(2018)092
clear in section 4). Indeed the most up-to-date determinations of m u /m d indicate that it is very close 1/2, making extreme pion-phobia a realistic possibility.

Constraints amenable to model-building
We start by discussing the main experimental observables that have excluded generic QCD axions in the MeV mass window. Despite being severe, these constraints can be "modelbuilt away" by deviating from generic models.
If the axion couples to heavy quark flavors, such as charm or bottom, it is strongly constrained by radiative decays of quarkonia, such as J/ψ → γa and Υ → γa. Wilczek [36] showed that such decay widths can be related to the leptonic widths via: where λ c and λ b are the axion couplings toc iγ 5 c andb iγ 5 b, respectively, and C J/ψ , C Υ ∼ O(1) encode QCD and relativistic corrections [37]. The MeV mass range for the QCD axion corresponds to decay constants f a in the O(1 − 10) GeV range, and consequently to large couplings of the axion to charm and/or bottom quarks, namely, λ c ∼ O(m c /f a ) and/or λ b ∼ O(m b /f a ). With such large couplings, radiative decays of quarkonia to γa would dominate over leptonic modes, in gross contradiction with observation. In fact, bounds from quarkonia decays alone [38][39][40][41][42] were sufficient to exclude the entire parameter space of the original PQWW axion (see, for instance, [43]). These and other bounds led to much activity during the 1980's in model building QCD axion variants, i.e., variations of the original PQWW axion that could evade existing constraints at the time and remain viable. A class of visible axion variants that trivially evades quarkonia bounds are those that couple the axion exclusively to first generation quarks. 2 Nevertheless, they still have to contend with other constraints. For instance, if such axions have suppressed couplings to leptons and decay dominantly to a pair of photons, they are sufficiently long lived and hence robustly excluded by a variety of beam dump experiments [44][45][46][47] in the range 100 keV m a 30 MeV. Limits from beam dumps are substantially degraded if axions couple to electrons with strength O(m e /f a ). In this case, axions heavier than 5-10 MeV become very shortly lived, τ a 10 −13 s, and their decay products stop in the earth shielding before reaching the detector, evading this class of bounds. Short-lived axions decaying to e + e − also evade severe constraints from K + → π + (a → invisible), whose branching ratio is bounded to be 4.5 × 10 −11 [48][49][50].
While coupling the axion to electrons is desirable, an analogous coupling to muons with strength O(m µ /f a ) would induce contributions to (g − 2) µ that would violate present bounds unless f a v EW . Therefore, in this work we shall restrict ourselves to variant axion models that couple exclusively the first generation fermions, namely, u, d, and e.

Pion-phobia
The mixing of the axion with the neutral pion poses a major challenge to the viability MeV axion models. Severe constraints on this mixing have been placed three decades ago, and the idea of avoiding them via pion-phobia is just as old [51].
Axion-pion mixing induces the rare decay process π + → e + ν e a → e + ν e e + e − . The width is given by [52]: where θ aπ is the axion-pion mixing angle, and θ c is the Cabibbo angle. The SINDRUM collaboration [53] searched for this specific decay, and put bounds on Br(π + → e + ν e (a → e + e − )) ranging from (0.5 − 1) × 10 −10 in the mass range m a ∼ (1 − 20) MeV. Using (2.3), this translates into a severe upper limit on θ aπ : As we shall see in section 4, typical MeV axion variants are in conflict with (2.4), since they predict θ aπ ∼ O(f π /f a ) ∼ (0.5 − 10) × 10 −2 . As a consequence, the experimental upper bound on θ aπ pushes viable models into special regions of parameter space where the axion is pion-phobic, and therefore has suppressed couplings to isovector currents.

Non-robust/uncertain constraints
Finally, some constraints are plagued by large hadronic uncertainties, and until those are under better control no robust exclusion claim can be made. That is the case of a − η and a − η mixings, which provide the dominant contribution to K + → π + a. Unlike K + → π + (a → invisible), the rare decay K + → π + (a → e + e − ) is much less constrained, Br (K + → π + (a → e + e − )) O(10 −6 ). In light of the allowed range (2.4) for θ aπ , the contribution to this decay from π 0 − a mixing easily satisfies the experimental bound. However, it has been claimed in the literature [54,55] that the contribution to this amplitude from a − η mixing is octet (∆I = 1/2) enhanced, and translates into a bound θ aη O(10 −4 ). Based on naïve estimates of θ aη from leading order in chiral pertubation theory, it was then concluded that axions in the MeV mass range were hopelessly ruled out by K + decay bounds. In section 5, we will revisit these statements and constraints, and argue that due to the uncertainties involved in these estimates, previous claims of exclusion were overstated, and bounds from K + decays alone cannot definitively rule out the class of axion variants we consider. Table 1 summarizes the main experimental constraints relevant for generic MeV axions, and also lists the reasons as to why the variant model we shall introduce avoids all present bounds. This is not a standalone table -we urge the reader to follow our arguments in the main body of the paper to become fully aware of all the assumptions and caveats implicit in the information contained in Br(a → χχ) ∼ (Q PQ χ mχ/Q PQ e me) 2 MB (assuming τa < 10 −12 s) χ = ν (e,µ,τ,s) , sub-MeV DM, . . . Table 1. Summary of the most relevant existing bounds on MeV axions, as well as the conditions for experimental viability discussed in this paper. The 3rd column indicates the typical range in generic axion models for the quantity being constrained, estimated at leading order (LO) in chiral perturbation theory (χPT). The 4th column contrasts it with the corresponding expectation for a short-lived, pion-phobic MeV axion which, as we argue in the text, remains compatible with experimental bounds. The 4th column takes into account large uncertainties from next-to-leading order (NLO) in χPT. The fifth column indicates the reason for the difference in expected properties of generic models (3rd column) and of the viable variant (4th column), where: MB indicates that the quantity in the 4th column is achievable via model building; NLO indicates that the error in the quantity in the 4th column comes from large corrections at NLO in χPT; and PP indicates that the quantity in the 4th column is suppressed due to accidental cancellations at LO in the axion-pion mixing angle θ aπ , i.e., pion-phobia.

Axion-meson mixings in χPT @ LO
We can trivially evade the constraints discussed in section 2.1 by restricting our investigation to QCD axion variants coupled only to the first generation of SM fermions. This assumption is implied for the remainder of this paper.
We are now in a position to introduce more concrete notation and review the extraction of axion-meson mixing angles to leading order in chiral perturbation theory (χPT).
We define the axion couplings slightly above the QCD scale, where all heavy fermions (c, b and t) have been integrated out. We also choose a specific basis (motivated by UV completions to be discussed in section 9) in which the operators (a/f a ) GG and ∂ µ (a/f a )J 5 µ SM are not present, or have been rotated away. With this choice of basis, the variant models we shall consider are unambiguously defined by the following couplings: where the superscript (0) on the mixing angles indicates that these expressions stem from leading order in χPT. Note, in particular, that when any of the quark masses is taken to zero, both the axion mass, m a , and its mixing with the SU(3) χ singlet, θ aη 0 , vanish. This is consistent with the fact that, in the limit of a massless quark, the axion becomes a bona fide Goldstone boson associated with a non-anomalous conserved current. Previous studies in the literature have assumed that these leading order estimations of the axion-meson mixing angles were good enough approximations. Based on those, they went on to extract rates for processes involving the axion, and infer bounds that widely excluded all axion variants with m a O(MeV). The core argument of this paper is that these assumptions, while correct for a wide class of axion variants, fail under special circumstances and invalidate broad exclusion claims. In these special circumstances, the mixing angles in (3.7)-(3.9) receive O(1) corrections from the next order in the chiral expansion, introducing large uncertainties to the axion couplings and dramatically affecting the inference of bounds that depend on these mixing angles.

Constraints from pion decays
As discussed in section 2.2, bounds from π + → e + ν e a severely restrict viable MeV axions to the extreme pion-phobic region. In this section we will discuss under what conditions pionphobia can be achieved, as well as implications for another relevant pion decay, namely, The leading order expression for the axion mixing with the neutral pion in (3.7) indicates that, for generic O(1) values of Q u,d , θ aπ ∼ O(f π /f a ). For this generic situation, the axion-pion mixing is about one to two orders of magnitude larger than the other axion-meson mixings, and therefore dominates the axion phenomenology involving hadronic couplings.
However, notice that |θ is not bounded from below. In particular, for Indeed, plugging in the PDG's 2017 weighted avarage for the light quark mass ratio [10], m u /m d = 0.483 ± 0.027, and setting Q d = 1 without loss of generality, (3.7) gives: which contains the pion-phobic range (2.4) still allowed by π + → e + ν e a. Therefore, extreme pion-phobia might be achievable in variant models with Q d /Q u = 1/2 ≈ m u /m d . Unfortunately, under these circumstances, the leading order expression (3.7) is not a reliable estimate of θ aπ . The reason is that the second order expansion in χPT can give comparable contributions to θ aπ , which should be included for a reliable estimate of this quantity. However, since the coefficients of many O(p 4 ) χPT operators that contribute to θ aπ are poorly known, a precise determination of θ aπ is not possible. Moreover, as we JHEP07(2018)092 will discuss in section 9, additional GeV states from the UV completion of such models also give model dependent contributions to θ aπ .
In summary, we have identified a particular class of axion variants, defined in (3.1) with Q u /Q d = 2, which is compatible with the condition of extreme pion-phobia, and hence remains a viable possibility. Due to present errors in the determination of m u /m d , uncertainties in the second order expansion in χPT, and model dependence of associated UV completions, the estimated range for the axion-pion mixing (4.2) is much broader than the experimentally allowed range in (2.4), but it is nonetheless consistent with those bounds given errors.

The KTeV anomaly
Were such pion-phobic axion to exist, θ aπ would ideally be determined from experiment. Another very sensitive probe of this mixing angle is the π 0 decay width to e + e − , which is highly suppressed in the Standard Model. In fact, there appears to be a discrepancy in the observed π 0 decay width to e + e − at the ∼ 3σ level. The most precise measurement of this branching ratio, from the KTeV-E778 collaboration [56], is: where the first and second errors are statistical and systematic, respectively. Theoretical estimates of this branching ratio [57][58][59] in the SM predict: corresponding to a 3.2σ discrepancy between theory and experiment. 4 In the SM, this decay is loop induced via the π 0 γγ coupling. The imaginary (absorptive) part of this amplitude is well understood and gives an irreducible contribution, the so called "unitarity bound", to this branching ratio [62]: where β ≡ 1 − 4m 2 e /m 2 π . By noting that the width for π 0 → γγ is given by: we can model the imaginary part of the SM amplitude ImM(π 0 → e + e − ) SM via an effective π 0 e + e − vertex given by [63]: (4.7)

JHEP07(2018)092
Parameterizing BSM corrections to this coupling as y BSM πee , we can finally write L πee = i y SM πee + y BSM πee π 0 eγ 5 e , (4.8) and predict the contribution to the width from the imaginary part of the amplitude: As a conservative bound on the BSM contribution to Γ(π 0 → e + e − ), we can demand that (4.9) does not exceed the observed value (4.3) by more than two standard deviations, leading to: (4.10) Note that when y BSM πee < 0, the BSM amplitude destructively interferes with the SM one, and in the range y BSM πee ∈ [−5.2, 0] × 10 −7 the contribution (4.9) to the width from the imaginary part of the amplitude is smaller than the SM unitarity bound (4.5). In that case, the real (dispersive) part of the amplitude would have to be unexpectedely large to account for the observed width.
If instead we attempt to explain the excess, 5 i.e., fit the BSM contribution to account for the discrepancy between (4.3) and (4.4), we obtain the following 2σ range: So far we have kept this discussion general, since for some UV completions of the axion variants we are considering, not only the axion, but also other GeV states, can contribute to y BSM πee . Nevertheless, we can consider the limit of these models in which the only important contribution to y BSM πee comes from the QCD axion. In this case, (4.12) The 2σ bound (4.10) in this instance then becomes: (4.13) and the range needed to fit the KTeV anomaly within 2 standard deviations is: (4.14) Remarkably, the θ aπ fit to the KTeV anomaly, (4.14), is compatible with the pion-phobic range imposed by π + → e + ν e a.

Constraints from charged kaon decays
It was widely claimed in the literature [54,55] that bounds from K + → π + a ruled out all QCD axion variants in the MeV mass range. In this section we critically re-examine the arguments which have led to this claim. We find that there are large uncertainties involved in obtaining K + → π + a, from the interpretation of experimental analyses, to assumptions in estimating the K + → π + (η * → a) amplitude, and finally to the derivation of the η − a mixing angle. We conclude that these uncertainties are significant enough to preclude a definitive exclusion of MeV axions from existing bounds on rare K + decays.

Experimental bounds on
Measurements of rare K + decays such as K + → π + νν, K + → π + γγ, and K + → π + e + e − have been improving over time for the past five decades or so. However, for the last 30 years, the low m e + e − region of K + → π + e + e − has been neglected from experimental scrutiny for BSM physics. This is due in part to the very large backgrounds from π 0 Dalitz decays, namely, K + → π + (π 0 → e + e − γ), which make the study of the low m e + e − region very difficult. To the best of our knowledge, the last time dedicated searches were performed in this region was in the 1980's, by two different experiments, one at KEK [65] and the other at BNL [66]. We shall now discuss these searches in some detail.
The KEK experiment E89 [65] was published in 1984. It consisted of a high resolution spectrograph that measured the momentum of charged pions from the decay of charged kaons at rest. They searched for a peak in the π + momentum distribution, which would be evidence of a two-body decay of the charged kaon, K + → π + X 0 , to a new pseudoscalar X 0 . Since they did not impose any vetoes nor requirements on the remaining decay products of K + , this search was sensitive to X 0 → anything. Unfortunately, we believe there are a couple of issues with this study.
First, it is not clear what was the lowest value of m X 0 to which this search was sensitive. The abstract states that "bounds are presented for the mass range of X 0 from 10 to 300 MeV/c 2 ." The concluding paragraph, on the other hand, quotes 5 − 300 MeV as the range of exclusion for m X 0 . The acceptance curve for the pion momentum, shown in figure 1 of [65], has a lower range of P π + ∼ 224 MeV, corresponding to a lower range of 50 MeV for m X 0 . Finally, the exclusion curve presented in figure 2, of Br(K + → π + X 0 ) versus m X 0 , has a lower range of m X 0 = 10 MeV.
Moreover, there are issues with figure 2. It misrepresents limits from previous analysis, namely, Asano et al. [67] and Abrams et al. [68]. The limits from Asano et al. on Br(K + → π + (X 0 → γγ) ) are depicted in figure 2 as a factor of ∼ 5 weaker than in [67], and shown only in the range m X 0 ∈ (50 − 90 MeV), when in fact Asano et al.'s limits apply from (0 − 100 MeV) as long as τ X 0 < 10 −9 s. It is unclear which limit from Abrams et al. [68] is being depicted in figure 2, but assuming it is the limit on Br(K + → π + γγ) in the region K π + < 70 MeV, corresponding to P π + > 156 MeV (M γγ > 237 MeV), then figure 2 of [65] misrepresents it by a factor of ∼ 5 as well.

JHEP07(2018)092
An earlier conference note [69] presenting preliminary results of the KEK E89 analysis might shed light on the origin of these discrepancies. 6 In figure 6 of conference note [69], limits were presented in terms of the relative branching ratio Br(K + → π + X 0 )/Br(K + → π + π 0 ) ≈ 5 × Br(K + → π + X 0 ), and represented more accurately the results from Asano et al. [67] and Abrams et al. [68]. This suggests that the factor of ∼ 5 discrepancy in the final publication [65] was a typo in recasting the plot from relative branching ratio to absolute branching ratio. Adding to the ambiguity on the lowest X 0 mass to which KEK E89 was sensitive, the limits presented on figure 6 of the conference note [69] did not extend below The second issue with the KEK E89 analysis was their statistical inference of exclusion limits in the low statistics region m X 0 80 MeV. The formula used for obtaining limits is appropriate in the case of a "bump-hunt" on top of a very large background, where the errors are gaussian and scale as N background . It fails, however, in the Poisson limit where the number of expected events from signal plus background is O(1), where it can overestimate the exclusion power by a factor of several.
All things considered, we conservatively choose to ignore the limits from KEK E89 analysis [65] in the region m X 0 50 MeV.
That leaves us with the other dedicated search for the rare decay K + → π + (a → e + e − ) in the low m e + e − region, published in 1987 by Baker et al. [66]. Their apparatus, located at BNL's AGS, was optimized for the decay of K + to three charged tracks, with efficient discrimination between electrons, pions and muons. They collected a large sample of approximately 2.8 × 10 4 K + → π + e + e − events, mostly Dalitz decays, K + → π + (π 0 → e + e − γ), and applied optimized cuts to remove this background while remaining sensitive to the K + → π + a signal. While they claim that their results exclude branching ratios Br(K + → π + a) 0.8 × 10 −7 for any mass in the range m a ∈ (1.8 − 100 MeV) at 90% confidence level, we believe their exclusion is too aggressive for two reasons.
Firstly, in their final signal region, they fit the observed m e + e − distribution in the range 1.8 − 100 MeV as a constant background not properly modeled by their Monte Carlo (MC), and then subtract it before extracting limits. We believe that a more conservative approach when deriving limits would be to treat any residual events in the signal region as potential signal. We can estimate what those more conservative limits would be based on the data shown in figure 2(d) of [66], together with the assumption that the signal efficiency as a function of m a does not vary substantially from the one quoted at m a = 1.8 MeV, and trusting their Monte Carlo signal yield of 5 expected events for Br(K + → π + a) = 10 −6 . Under these assumptions, we infer the limit Br( Our second concern regards this analysis' sole reliance on Monte Carlo to estimate all the acceptances of both background and signal. As figure 2(d) of [66] shows, their Monte Carlo grossly misestimates the Dalitz background contamination in the signal region by about one order of magnitude. They do not address this mismodeling and simply re-scale JHEP07(2018)092 their expected background to match the observed rate. More importantly, they do not address whether their Monte Carlo estimation of the signal acceptance could possibily be mismodeled as well. Were their signal acceptance off by a similar order of magnitude, their limits could be weakened to as much as Br(K + → π + a) 10 −5 in the low m e + e − region.
We therefore choose to remain agnostic about the precise experimental limit on MeV axions from kaon decays, and for the remainder of this paper we discuss implications from K + decays to axion-meson mixing parameters assuming the following range: We stress that new experimental studies of the low m ee kinematic region in K + → π + e + e − are relevant and warranted, considering the sensitivity of this final state to rare new physics processes, and the fact that measurements in this region could be greatly improved by modern experiments.
Since the π 0 is on-shell in the SM process K + → π + π 0 , it is trivial to relate the amplitude for this decay to the amplitude for K + → π + a via axion-pion mixing: If (5.2) were the dominant contribution to M(K + → π + a), we could plug in the upper bound (2.4) on θ aπ to obtain: which is at least three orders of magnitude below the existing bounds discussed in section 5.1. Antoniadis et al. [54] and Bardeen et al. [55] pointed out that the axion mixing with η and η also contribute to M(K + → π + a): (5.4) and that, due to ∆I = 1/2 enhancement (a.k.a. octet enhancement) of these amplitudes, these contributions dominate the branching ratio for K + → π + a and lead to tension with experimental constraints.
Let us briefly review the arguments of [55]. In order to obtain the amplitudes M(K + → π + η 8,0 ), [55] considers the leading order ∆S = 1 chiral Lagrangian describing non-leptonic kaon decays: where λ a are the Gell-Mann matrices, λ ds ≡ (λ 6 + iλ 7 )/2, and the Clebsch-Gordan coefficients C ab are given for instance in [70]. The first operator transforms as an octet under chiral SU(3) χ , while the second transforms as a 27-plet. In particular, while M(K 0 → ππ) receives contributions from both operators, only the 27-plet operator contributes to M(K + → π + π 0 ): Above, δ 0 , δ 2 are strong interaction S-wave ππ phase shifts. Within the effective framework of χPT, the coefficients g 8 and g 27 cannot be obtained from first principles. However, under the standard assumption that (5.6), (5.7), (5.8) provide the leading contributions nonleptonic kaon decays, g 8 and g 27 can be fit to match the observed K → ππ amplitudes [70] (see also [71]): Whereas naïvely one would expect g 8 ∼ g 27 , (5.9) shows a large enhancement of the octet coefficient relative to the 27-plet's, g 8 /g 27 31.2, tied to the fact that the hadronic width of K 0 S is much broader than that of K ± , e.g., This is known as the octet enhancement in non-leptonic kaon decays.
Ref. [55] noted that the octet operator in (5.5) also contributes to the (off-shell) amplitudes: Defining the quark flavor basis: and neglecting terms O(g 27 /g 8 ) and O(m 2 π /m 2 K ), we can then use (5.11), (5.12), (5.4) and (5.6) to obtain the following relation: 7 To be precise, the definition of axion mixing angles with states that are not mass eigenstates, such as

JHEP07(2018)092
Ref. [54] noted that expression (5.17) does not take into account the absence of strong final-state interactions between π + and a. Following [54], we correct this by introducing a fudge factor D ππ ∼ 1/ √ 3 multiplying the r.h.s. of (5.17). We then finally obtain: Finally, using the experimental upper bound (5.1), we infer the following constraint: We would like to stress that the arguments above, leading to (5.18) and (5.19), hinge on the assumption that the O(p 2 ) octet operator, shown in (5.5), is enhanced. While this assumption is not considered controversial, it is nonetheless possible that this is not the correct description of octet enhancement in non-leptonic kaon decays [72][73][74]. In particular, an equally good description is obtained by transfering the enhancement to the O(p 4 ) expansion of L ∆S=1 χ . Consider, for instance, where Λ ∼ 4πf π is a natural cut-off. If one assumes that the O(p 2 ) coefficients g 8 and g 27 are comparable, and instead the O(p 4 ) coefficient g 8 is enhanced, one obtains an equally good phenomenological fit of K → ππ and K → πππ data compared to the standard fit in (5.9). 8 However, unlike the octet operator in (5.5), (5.20) does not contribute to K + → π + η 8,0 . Therefore, the implication of the alternative fit in (5.21) would be that the amplitude for K + → π + a, which is proportional to g 8 , but unaffected by g 8 , would not be octet enhanced.
In that case, one would have to properly redo the fit to kaon data to obtain the precise η ud and ηs, goes as follows. Consider the interactions in the (canonically normalized) quark flavor basis: where J ud and Js are external sources. Integrating η ud and ηs out, one obtains:

JHEP07(2018)092
value of g 8 , but as a rough estimate, taking g 8 ∼ O(g 27 ), the constraint on a − η mixing would be weakened by a factor of ∼ 30, In short, the enhancement of the amplitude M(K + → π + a) depends on how octet enhancement is realized in χPT. If the amplitudes M(K + → π + η 8,0 ) are not octet enhanced, the resulting bound (5.22) on θ aη ud is relatively weak and does not presently pose a threat to the viability of MeV axions. Under the more conventional assumption of octet enhancement at O(p 2 ), however, it becomes important to establish how compatible is θ aη ud with the bound in (5.19), which we shall now address.

Axion-eta mixing
The last sources of uncertainty in predicting the rare decay K + → π + a, namely, the mixing angles θ aη 8 and θ aη 0 , have been the least critically examined in the axion literature.
It is well known that the leading order expansion in χPT does not adequately describe the η and η masses and mixing angles [77,78]. Indeed, the second order expansion of the Chiral Lagrangian provides important corrections to masses, decay constants and mixing angles of singlet and octet mesons, which typically scale as: Above, L i are the dimensionless coefficients of the O(p 4 ) operators in the chiral expansion, and are commonly known as Low Energy Constants (LECs) [79,80]. Many LECs are reasonably well-determined from experimental and/or lattice data, their typical size being L i ∼ O(10 −3 ). From (5.23) it is then evident that these encode O(1) effects in η-η mixing, and may very well have comparable importance in describing a-η and a-η mixing.
In order to illustrate the uncertainties involved in obtaining θ a η,η , we consider Leutwyler's study of η-η mixing in [81], which, based on 1/N c -expansion counting rules, retained only the following operators: 9 +OZI violating terms.

JHEP07(2018)092
Remembering that the axion is formally a phase of the light quark mass matrix M q (see (3.3)), and setting Q d = Q u /2 = 1, we can expand (5.24) to obtain: Above, we have expanded the coefficients in powers of m 2 π m 2 K and kept only the leading terms; we have also omitted all π and K dependent terms, since they are not relevant for the present discussion. Furthermore, we have simplified the notation by defining: Canonically normalizing the kinetic terms and integrating out η ud and η s , we obtain: Plugging in the numerical fit of the LECs obtained by [81] yields: Rather than claiming that we have calculated θ aη ud more precisely, our purpose with this exercise is to make the point that contributions from the second order chiral expansion change the leading order estimate, θ (0) aη ud , by O(1). While we have restricted this exercise to a couple of LECs at tree level, inclusion of the full O(p 4 ) expansion and loop corrections will most likely change the answer from (5.27). Had we derived such expression, even then we would not be able to evaluate it numerically because many of the LECs contributing to θ aη ud are still undetermined. For instance, operators such as [82] also give important contributions θ aη ud and θ aηs . While such operators also affect η-η mixing, and hence are subject to constraints, there are too many independent LECs to be fit by observations, and without further assumptions, the parameterization of the η-η phenomenology in χPT is underconstrained.

JHEP07(2018)092
We conclude, therefore, that currently it is impossible to reliably estimate θ aη ud and θ aηs , given that these mixings angles can receive O(1) corrections from the second order chiral expansion whose numerical values are undetermined by existing data. Therefore, one cannot claim with confidence that θ aη ud violates the K + → π + a bounds in (5.19) or (5.22), and therefore K + decay bounds by themselves do not provide a definitive exclusion of the QCD axion variant we are considering. Most likely, a reliable determination of the mixing angles θ aη ud and θ aηs would have to come from lattice calculations, or, were such axion to be realized in nature, from direct measurements.

The physical axion current
Because of much confusion in the literature, in this section we review the differences between the PQ current, the Bardeen-Tye current, and the current associated with the axion mass eigenstate. Identifying the latter, in particular, is a critical part of working out the proper phenomenology of axions, since it is the physical axion current that determines the coupling of the axion to the chiral U(3) χ currents that mediate meson decays, nuclear de-excitations, as well as the electromagnetic anomaly that induces the axion coupling to two photons.
In the UV, the PQ current for the class of axion variants defined in (3.1) is given by: 10 Bardeen and Tye [83] pointed out that the PQ current cannot be associated with the physical axion, i.e., the light mass eigenstate, since the divergence of this current receives a large contribution from the strong anomaly: which does not vanish in the limit m q → 0. The physical axion, on the other hand, becomes a bona fide Goldstone boson in the limit of a massless quark, and should therefore be associated with a current that is conserved as m q → 0. In order to identify the physical current, [83] subtracted the strong anomaly from J PQ , obtaining the Bardeen-Tye current: which is conserved in the limit m q → 0:

JHEP07(2018)092
We can rewrite J BT in order to show its explicit dependence on the chiral U(3) χ currents. First we define: Above, π 3 , η 8 and η 0 are not mass eigenstates.
The Bardeen-Tye current in (6.3) can then be written as: (6.7) The Bardeen-Tye current has been identified in numerous studies as the physical axion current. Under this assumption, it is straightforward to read the mixing angles θ a π/η from (6.5) and (6.6): θ BT aη 0 = 0 . (6.10) Note, however, that this assumption is only strictly correct in the limit where the η is decoupled, i.e., in the limit M 0 → ∞ (see (3.2)) and η → η 0 . Indeed, the Bardeen-Tye prescription to obtain the mixing angles θ a π/η is equivalent to the 1 st order χPT prescription with η 0 decoupled. This is corroborated by the fact that (3.7), (3.8), (3.9) and (6.8), (6.9), (6.10) agree in the limit M 0 → ∞.
We point out as a side remark that the condition that ∂ µ J BT µ vanishes in the limit m q → 0 is not sufficient to uniquely determine the Bardeen-Tye current. In particular, one could define: 5 µ + c 8 J which also satisfies the aforementioned condition with arbitrary finite coefficients c 3,8,0 . Again, in this context, this ambiguity can only be removed in the limit of decoupled η . In this limit, c 0 in (6.11) must necessarily vanish, since θ aη 0 ∝ c 0 → 0 as m η → ∞. 11

JHEP07(2018)092
Moreover, in this limit the mass eigenstates a phys , π and η are entirely determined by the non-anomalous currents J BT , J (3) and J (8) . Imposing then that J BT commutes with J π and J η , one can estimate c 3 and c 8 via current algebra methods, and conclude that these coefficients provide only small corrections to (6.8) and (6.9) of order m 2 a /m 2 π and m 2 a /m 2 η , respectively.
Unfortunately, once the η is included in the spectrum, the Bardeen-Tye prescription is no longer sufficient to determine the physical axion current. In particular, while π 3 in (6.5a) is a good approximation to the mass eigenstate π 0 , the physical states η and η have substantial components of both η 8 and η 0 , and hence are not amenable to current algebra methods, which are not applicable to anomalous currents.
Therefore, generically the physical axion current has components from all three neutral chiral U(3) χ currents: (6.12) The coefficients ξ a3 , ξ a8 and ξ a0 in (6.12) are related to the axion-meson mixing angles by: This discussion provides yet another perspective on the difficulties in estimating the hadronic mixing angles of the axion. This propagates into uncertainties in the axion's phenomenology, such as rates of rare meson decays, discussed in the previous two sections, as well as the axion couplings to photons and nucleons, which we will now consider.

Axion-photon coupling
As light quarks confine below the QCD scale, an effective operator coupling the axion to the electromagnetic dual field strength is generated. In this section we derive this operator, estimate the axion decay branching ratio to a pair of photons, and comment on bounds from K + → π + (a → γγ).
We start by rewriting the physical axion current (6.12) as: where

JHEP07(2018)092
With this notation, we can derive the electromagnetic anomaly of J a phys µ and the axion coupling to photons straightforwardly: where we assume that m a > 2 m e , so that the coupling of the axion to electrons does not contribute to this effective operator. The axion decay width to γγ is then given by: . (6.19) This expression neglects the loop induced contribution coming from the axion coupling to electrons, which becomes comparable to the hadronic contribution above if the linear combination of mixing angles in (6.19) is O(10 −4 ). The relative γγ to e + e − branching ratio is then: where β e = 1 − 4 m 2 e /m 2 a . Note the very weak dependence of this branching ratio on the axion mass when m a m e . From (6.20) we can estimate the branching ratio for the rare decay K + → π + (a → γγ) in terms Br(K + → π + (a → e + e − )): Br(K + → π + (a → γγ)) ≈ 10 −11 θ aπ + 5 3 θ aη ud + √ 2 3 θ aη s 10 −2 2 Br(K + → π + (a → e + e − )) 10 −6 , which is still below current experimental sensitivity by two or three orders of magnitude [84].

Axion-nucleon couplings and nuclear de-excitations
Another classic signature of visible axions is axion emission in nuclear de-excitations, first studied in [85][86][87]. We shall briefly review here the standard results from the literature. From (6.12) or (6.14) we can infer the axion coupling to protons and neutrons: A τ 3 + where ∆q is defined as: with s µ being the nucleon spin-vector. In terms of ∆q, we can alternatively recast g aN N as: The pion-nucleon coupling, set by g A , is well determined from nuclear β-decay, g The other axial coupling constants, g A and g (0) A , are much harder to extract from experiment, and various estimates based on data from semi-leptonic hyperon decays, proton deep inelastic scattering, and lattice calculations yield the following ranges [88][89][90][91][92][93][94][95][96][97][98][99][100] These uncertainties, compounded by the difficulties in extracting the axion hadronic mixing angles θ a π/η , are obstacles to making firm predictions for nuclear de-excitation rates via axion emission. 12 Nevertheless, the parametric dependence of such rates on nuclear and axion properties is relatively well understood. Donnelly et al. [86] showed that the axion acts as a "magnetic photon" in nuclear transitions, and using standard multipole techniques in the 12 Ref. [101] also considered properties of the axion in great detail, extracting the axion mass and coupling to photons and nucleons by combining results from Lattice QCD (LQCD) and χPT at NLO. We believe that the authors of [101] were overly optimistic regarding the smallness of the errors in some of the LECs used, as well as the errors in LQCD extractions of the nucleon spin content ∆q. But more importantly, their approach differs substantially from ours. For instance, they associate the axion with the Bardeen-Tye current in the two-flavor effective theory (i.e., chiral EFT with the strange quark integrated out). While the inclusion of NLO contributions in χPT should in principle correct for the use of the (non physical) Bardeen-Tye current, it is unclear to us whether the authors of [101] did it in a consistent manner such that the resulting axion couplings to photons and nucleons would indeed correspond to that of the mass eigenstate. For instance, applying their procedure to our model to extract the couplings of the axion to nucleons, one would be led to conclude that our axion variant does not couple to nucleons, which is incorrect. (In particular, in their notation, our model would correspond to: (1) of [101]; this would translate into cu = c d = c s,c,b,t = 0 in eq. (46) of [101]). A thorough comparison of our treatment with that of [101] is a non-trivial task which is beyond the scope of this paper.

JHEP07(2018)092
long-wavelength limit, derived the ratio of axion to photon de-excitation rates of a nuclear state. For an isoscalar (∆T = 0), p-wave magnetic (M1) transition, [86] derived: Above, ∆E is the energy splitting between the two nuclear levels. The parameters µ (0) ≈ 0.88 and η (0) ≈ 1/2 are related to the nuclear magnetic moment and the ratio of convection to magnetization currents, respectively. For an isovector (∆T = 1) M1 transition, the analogous expression is [86]: For mixed isospin transitions, the generalization is [102]: where c 0,1 are the probability amplitudes for the different isospin components. During the 80's, several experiments have searched for variants of the visible QCD axion in nuclear de-excitations [102][103][104][105][106][107][108][109][110]. Since most nuclear levels have splittings of a few MeV, the results of these experiments were only relevant for a range of m a already ruled-out by other measurements, in particular beam dump experiments (see section 7). In 1990, de Boer et al. [111] studied transitions of the 8 Be nucleus which were energetic enough to probe axion masses of up to ≈ 15 MeV. A few years later, however, the results of this experiment were revisited by the authors in [112,113] and the claimed limits were weakened substantially. We shall therefore not consider the results of [111][112][113] here.
Most recently, the ATOMKI collaboration [114] has measured several 8 Be nuclear transitions via emission of e + e − pairs. The two relevant transitions for our discussion are the M1 de-excitations of the J P = 1 + isospin doublet states, namely, 8  The ATOMKI collaboration claimed to have observed a deviation in the e + e − spectrum of the 8 Be * (18.15) → 8 Be(0) transition relative to the SM prediction of internal pair conversion (γ * → e + e − ). According to [114], this deviation was consistent with the onshell emission of a narrow resonance X of mass m X ≈ (16.6 ± 0.9) MeV promptly decaying

JHEP07(2018)092
to e + e − . The best fit for the relative de-excitation rate was Γ X /Γ γ ≈ 5.8 × 10 −6 , with a statistical significance of 6.8 σ. Moreover, in the original publication [114], no excess was observed in the e + e − spectrum of the 8 Be * (17.64) → 8 Be(0) transition. No error bars were quoted for either measurement, neither were upper bounds on emission rates of generic new particles, such as light vectors or pseudoscalars. Subsequent studies have attempted to understand these results via nuclear physics models [115], or via emission of a new light resonance [116][117][118][119][120][121][122][123].
The ATOMKI collaboration is revisiting these measurements with an improved experimental set-up [124,125]. At the time of writing of this paper, preliminary results by Krasznahorkay et al. have been released [125] indicating a possible excess in the 8 Be * (17.64) → 8 Be(0) transition as well, which was fit by a similar hypothetical particle X with mass m X = (17.0 ± 0.2) MeV, at a rate of Γ X /Γ γ ≈ 4.0 × 10 −6 .
At present, we believe the claims from the ATOMKI collaboration are inconclusive, and an independent measurement is warranted. Nonetheless, we can make order of magnitude predictions for axion emission in 8 Be M1 transitions using (6.35). While the excited states 8 Be * (17.64) and 8 Be * (18.15) are predominantly T = 1 and T = 0, respectively, their widths strongly indicate that they are isospin-mixed: Estimates on the level of isospin mixing vary between 0.18 sin θ 1 + 0.43 [118,122,126].
Taking the value of sin θ 1 + = 0.35 suggested in [122], we obtain the following range for the 8 Be * (18.15) axion de-excitation rate: Above, we also assume m a = 16.6 MeV and θ aπ = 0.5 × 10 −4 . Reiterating our point, while uncertainties in nuclear and axion parameters preclude us from making a precise prediction for the rate of axion emission in the de-excitation of 8 Be * (18.15), the range in (6.40) is compatible with the excess observed by the ATOMKI collaboration. Analogously, for the 8 Be * (17.64) de-excitation we obtain: Here, despite uncertainties in (6.41), we are able to make the more firm prediction that, should an anomaly persist in the 8 Be * (17.64) de-excitation rate at the level of Γ X /Γ γ ≈ O(10 −6 ), this would preclude our QCD axion variant from offering as a viable explanation of this excess.

JHEP07(2018)092
7 Axion-electron coupling Unlike its hadronic couplings, the axion coupling to electrons bears no consequence to the solution of the strong CP problem, and therefore it is a model dependent parameter that can be adjusted according to phenomenological constraints. Moreover, it is much less susceptible to calculational uncertainites, and can be probed in a variety of existing and upcoming experiments. In fact, it is conceivable that this MeV axion variant could be definitively excluded via its electron couplings before any substantial progress is made regarding its hadronic phenomenology. In this section, we will review existing bounds on the axion-electron coupling, and briefly discuss on-going experimental efforts to search for hidden photons that could test the MeV axion scenario as well.
For specificity, we reiterate the axion coupling to electrons: Q e f a m e aē iγ 5 e , (7.1) as well as the relation between the axion's mass and decay constant: and the axion lifetime, assuming that it is dominated by a → e + e − : 3) The constraints discussed in this section are summarized in figure 1.

Beam dump constraints
In the 80's, several beam dump experiments have specifically targeted the QCD axion.
Since the results of these searches require no re-interpretation, we refer the reader to the original papers for details on production and detection mechanisms. In this subsection, we compile the most significant contraints in the region of parameter space of interest. Constraints from beam dumps can be avoided if axions are sufficiently short-lived so as to decay in the earth shielding, before reaching the detectors. Moreover, in order to remain experimentally viable, invisible decay modes of the axion must be subdominant by at least O(10 −4 ) in order to avoid stringent constraints from K + → π + (a → invisible). In order to fulfill these requirements, the axion must be heavier than a few MeV and couple significantly to electrons.
In figure 1, we show contours of the axion lifetime as a function of the axion mass and coupling to electrons, Q e /f a . We also show the most relevant beam dump constraints as shaded gray regions, namely, FNAL E774 [127], SLAC E141 [128], and Orsay [129]. Other overlapping but less powerful constraints, such as Bechis et al. [132], SLAC E56 [133], FNAL E605/772 [134,135], and KEK [136] are omitted in order to make the plot more readable. Experiments with longer earth shieldings, such as SLAC E137 [45] and CHARM [44] lie outside the plot.  [127][128][129], and the shaded orange and green regions are excluded by BaBar [130] and KLOE [131] searches for dark photons, respectively. Constraints from the electron's anomalous magnetic moment are model dependent and therefore not shown in this plot (see section 7.2 for details).

Constraints from (g − 2) e
Next, we discuss the MeV axion contribution to the electron's anomalous magnetic moment. Constraints coming from (g−2) e are difficult to pin down due to uncertainties in the axion's two-loop contribution, as well as the model dependence associated with UV completions of MeV axions. As we shall discuss in section 9, generic UV completions of such models contain additional particles that can also contribute substantially to (g − 2) e .
State-of-the-art calculations of a e ≡ (g − 2) e /2 in the SM predict [137]: where the first three uncertainties are theoretical, and the last one stems from the error in the measurement of the fine structure constant. The most precise measurementof a e to date [138,139], on the other hand, gives:

JHEP07(2018)092
The dominant contribution to this quantity from the MeV axion comes at one-loop and is given by: where we have used (7.2) in the second equality. The two-loop contribution (which is a Barr-Zee type diagram) suffers from large uncertainties due to the poor determination of g aγγ , as well as the precise value of the cut-off scale Λ above which g aγγ can no longer be treated as a point-like coupling [140]. Setting Λ ∼ 4πf π , we have [141][142][143]: where, from (6.18), and A few observations are in order: the one-loop contribution is negative and exceeds the allowed range in (7.6) by more than two standard deviations if |Q e | 1/2. However, the two-loop contribution has the opposite sign if the product Q e g aγγ is positive, and could in principle partially cancel the one-loop contribution to a substantial degree. In particular, if the linear combination of hadronic mixing angles in (7.9) is ∼ 2 × 10 −2 , then Q e 1.3 is still consistent with (7.6) at 2σ. Note that this range of mixing angles, and therefore of g aγγ , is consistent with other experimental constraints, mainly because of the very poor determination of θ aη s . The rare decay K + → π + (a → γγ) is not yet sensitive to θ aη s O(10 −1 ) (see section 6.1), and the contribution of θ aη s to the 8 Be nuclear transitions is mild if ∆s ∼ −0.02, which is suggested by recent lattice results [94][95][96][97][98][99][100].
Given all these considerations, existing measurements of the electron anomalous magnetic moment cannot provide a robust and model independent exclusion of the MeV axion parameter space, and therefore we do not include it in figure 1.

Constraints from searches for dark photons
Over the past decade, there has been a spur in phenomelogical studies of light, weakly coupled vectors ("dark photons"), mostly motivated by dark matter phenomenology [144][145][146][147].

JHEP07(2018)092
As a result, many experimental proposals have been put forth to search for dark photons in the near future [148][149][150][151]. As it turns out, dark photon production and detection mechanisms are very similar to those of light pseudoscalars such as the axion, and several constraints on dark photons decaying visibly can be easily recast as limits on the MeV axion.
In this section, we briefly comment on two experimental strategies that are relevant for MeV axions, namely, production in e + e − collisions, and at fixed-target experiments. We also translate existing bounds as well as future projections of proposed experiments in the MeV axion parameter space.
The axion can be produced in association with a photon in e + e − annihilation via the aēiγ 5 e vertex. In the limit that the center-of-mass energy √ s m a , the cross-section computed at leading order is: where θ γ is the angle of the final state photon with respect to the beam axis, and we show the explicit dependence on m e that regulates the collinear divergence when |sinθ γ | → 0.
Recent results from BaBar [130] and KLOE [131] looking for e + e − → γ (A → e + e − ) can be easily recast in the axion parameter space using (7.11). We display the corresponding limits in figure 1.
The axion can also be produced in fixed target experiments via "axion-bremsstrahlung" when an energetic electron scatters off of a heavy nuclear target. The differential crosssection as a function of the axion emission angle with respect to the beam, θ a , and the fraction x a of the incoming electron energy carried by the radiated axion, x a ≡ E a /E 0 , has been worked out in [152]. Below, we just quote the result: and χ encodes nuclear form factors (in the notation of [145], χ = Z 2 Log). For further details, see [152][153][154].
Much like the case of dark photon bremsstrahlung, axion emission is almost collinear, and its characteristic emission angle is parametrically smaller than the angle of its decay products with respect to the incident beam [145]. Integrating (7.12) over emission angles, we obtain which is more peaked towards x a = 1 than the corresponding x A dependence of dark photon emission. Finally, we quote the total production cross-section integrated over angle and energy: m e m π f π and we have used (7.2) in the second equality. We see, therefore, that the total crosssection for axion emission increases logarithmically with m a while m a √ m e E 0 , whereas it decreases logarithmically with m a in the regime √ m e E 0 m a E 0 . In figure 2, we translate the projected reach of several planned dark photon experiments into the MeV axion parameter space, using the following rule of thumb relating the dark photon kinetic mixing parameter A and the axion-electron coupling that would yield comparable signal strengths: (7.17) We emphasize that (7.17) is only approximate, since the angular dependence of pseudoscalar cross-sections differs from that of vectors. Figure 2 includes ongoing and proposed dark photon searches via e + e − annihilation (VEPP-3 at BINP [155,156], PADME in Frascati [157], MMAPS at Cornell [151] and Belle-II at KEK [158]), as well as dark photon JHEP07(2018)092 bremsstrahlung from electrons scattering off of heavy fixed targets (HPS [159] and Dark-Light [160,161] at JLab, MAGIX [162] at Mainz). These projections were based on the Dark Sectors 2016 Workshop community report [151]. We refer the reader to this report for further details.

Other constraints
Finally, we comment on other potentially constraining observables that could probe the MeV axion parameter space. We discuss them in lesser detail because the resulting limits are either weaker than the ones previously discussed, suffer from similar hadronic uncertainties, or lack experimental information specific to the axion related signature. Bounds from the hyperfine splitting between the 1 3 S 1 and 1 1 S 0 positronium levels [163,164], for instance, are currently not competitive with beam dumps contraints shown in figure 1. Neither are bounds from neutron-nucleus and neutron-electron scattering [165,166], since the axion mediated contribution is spin-dependent [167,168].
We emphasize that axion production via its coupling to photons 14 should be subdominant to the processes discussed in section 7.3 induced by the axion's electronic couplings. This can be easily seen by noting that the effective axion-photon coupling, g aγγ a FF , is suppressed by an effective scale of O(10 − 100) TeV: A comparable effective scale suppresses the couplings of the axion to Zγ, thereby suppressing the rate for the decay Z → γa below current experimental bounds. This rare Z decay is typically a sensitive probe of generic "axion-like particles" (ALPs) with low decay constants (see recent studies in [192][193][194]). A generic ALP φ with decay constant f φ couples to γ Z with typical strength g φγZ ∼ α 4πf φ . This parametric dependence does not apply for an axion with f a ∼ O(GeV) for the following reason: in the case of a generic ALP, the φ-γ-Z coupling is generated by integrating out electroweakly charged fermions 13 Note that EW penguins with charm or top quarks in the loop do not contribute in this case, because the axion does not couple to heavy quarks.
14 See [192] for a recent study of ALP photo-production. . So far we have assumed that the MeV axion does not couple to neutrinos. This can be naturally realized if neither the lepton doublets nor right-handed Dirac and/or sterile neutrinos carry PQ charge. If we relax this assumption, the typical axion coupling to neutrinos would scale as m ν /f a ∼ O(10 −11 − 10 −10 ), which is likely too small to be phenomenologically interesting. On the other hand, PQ-charged sterile neutrinos could couple to the axion much more strongly, namely, with strength O(m νs /f a ). If m νs is heavy enough, this could lead to observable signatures. A particularly well motivated mass range for sterile neutrinos is m νs ∼ O(1 − 10) keV, where sterile neutrinos can constitute warm dark matter and potentially address structure formation problems such as the core vs cusp and missing satelites problems (see [195] and references therein). In this mass range, the branching ratio for axion decay into sterile neutrinos would be: An interesting consequence of (8.2) would be a contribution to the rare kaon decay K + → π + a → π + ν s ν s at the level of: which is compatible with the strongest bounds set by BNL's E787/E949 experiments [48][49][50], namely, Br(K + → π + (X 0 → invisible)) 0.45 × 10 −10 for m X 0 70 MeV. CERN's NA62 experiment will soon supersede these bounds by at least one order of magnitude, since it is expected to have sensitivity to the SM prediction of Br(K + → π + νν)| SM = (8.4 ± 1.0) × 10 −11 [196] with better than 10% accuracy [197]. Other interesting possibilities for neutrino phenomenology could come from offdiagonal couplings of the axion to active (ν ) and sterile neutrinos (ν s ), such as λ a ν ν s , especially if the sterile neutrinos are light enough to be produced in laboratory experiments and/or astrophysical processes. Effects such as ν s → ν a → ν e + e − (if m νs > m a ), non-standard neutrino interactions, MSW-type resonance effects in neutrino propagation in matter, and neutrino transport in core-collapse supernova are interesting directions to explore, but are beyond the scope of this paper.

GeV scale completions of MeV axion models
The Peccei-Quinn breaking scale f a suggests that there is new dynamics around 1−10 GeV. While there are many constraints on what type of new particles may be associated with this new dynamics, there is still a large degree of model dependence in UV completions of the MeV axion. A thorough exploration of the viable phenomelogy is beyond the scope of JHEP07(2018)092 this paper, and will be deferred to a future publication [198]. In this section, we briefly illustrate a few possibilities with a simple toy-model.
Consider introducing two new complex scalar degrees of freedom, Φ u and Φ d , with PQ charges Q PQ Φu = Q u = 2 and Q PQ Φ d = Q d = 1, respectively. We can then UV complete the PQ mechanism at the GeV scale by writing: (9.1) The scalar potential V (Φ u , Φ d ) enforces the PQ symmetry, and induces vacuum expectation values for Φ u and Φ d , hence breaking the PQ symmetry: We can then decompose Φ u and Φ d into real scalar and pseudoscalar components: Above, tan β PQ ≡ f u /f d , the pseudoscalar a is the MeV axion, and the axion decay constant, f a , is given by In this toy model, we have introduced 3 extra degrees of freedom besides the MeV axion: two real scalars, ϕ u and ϕ d , and one pseudoscalar, η PQ . The natural scale for their masses is set by f a . Such light states must, therefore, be electroweak singlets in order to be phenomenologically viable. As a consequence, the couplings of Φ u and Φ d to SM fermions in (9.1) descend from higher dimensional operators and must be generated after electroweak symmetry breaking. We shall return to this point in the next section.
Because of their couplings to light quarks, the new states ϕ u , ϕ d , and η PQ are produced in hadronic interactions, and decay dominantly to hadrons -though with smaller cross-sections and narrower widths than typical QCD hadronic resonances. Unfortunately, estimating their widths and mixings with QCD resonances is challeging, since they lie in a regime where neither perturbative QCD nor chiral perturbation methods are reliable. It is conceivable, at least in principle, that these states are not excluded. The scalar resonances ϕ u and ϕ d , for instance, might not have been identified if lying in the murky mass range below 2 GeV [199], where many broad scalar resonances overlap and might constitute a large and complicated background to disentangle. The same argument is less likely to apply for the pseudoscalar resonance η PQ , unless it lies in a mass range where the hadronic pseudoscalar spectrum is poorly understood. For instance, the 1300-1500 MeV range contains three 0 −+ states, namely, η(1295), η(1405) and η(1475). The exact nature of these states is still subject to debate [200], with interpretations ranging from those being radial excitations of lighter pseudoscalars [201][202][203] to pseudoscalar glueballs [204][205][206][207][208][209][210][211][212]. Some authors JHEP07(2018)092 dispute the existence of the η(1295) state and claim that there is a single pseudoscalar meson in this mass range, the η(1440) state, which would be the first radial excitation of the η [213,214]. There are also claims that the η(1405) and η(1475) structures might originate from a single pole [215], the splitting being due to nodes in the decay amplitudes [213], or amplitude mixing via a triangular singularity [216,217]. Whatever the nature of η(1295), η(1405) and η(1475) may be, a priori they might present a challenging background for η PQ , and possibly mix substantially with it. A less speculative and more careful investigation of these possibilities will be done in [198].
We conclude by commenting on the coupling of these states to electrons. One possibility is that either Φ u or Φ d have electron Yukawa couplings, in which case all real degrees of freedom (namely, ϕ u , ϕ d , η PQ and a) will couple to electrons (either directly or through mixing), with typical strength O(m e /f a ). Alternatively, we may introduce a third PQ-charged complex scalar, Φ e = ϕ e e i ηe/fe , which will be responsible for generating the electron mass and the axion-electron coupling. Mixings between the Φ e and Φ u,d degrees of freedom will likewise result in couplings of ϕ u,d,e and η PQ ,e to electrons, although in this case some of the couplings might be suppressed or enhanced relative to m e /f a . Depending on the details of a specific model, these states might give important contributions to (g − 2) e , Γ(π 0 → e + e − ), and might be searched for in e + e − annihilation and fixed target experiments as well. A phenomenological study will be considered in [198].

EW scale completions of MeV axion models
As discussed in the previous section, the dynamics that breaks the PQ symmetry at the GeV scale generically requires new light degrees of freedom coupling to SU(2) W -charged quark bilinears. However, these new light particles themselves must not be SU(2) W -charged, otherwise they would be excluded, for instance, by measurements of the Z 0 width. Therefore, the Yukawa couplings of the PQ sector to SM fermions violate SU(2) W × U(1) Y and can only arise from higher dimensional operators after EW symmetry breaking. These operators can be generated by integrating out new degrees of freedom at the EW scale, such as heavy PQ-charged scalar doublets, or heavy vector-like quarks, leading to interesting and distinct LHC signatures. Although a thorough exploration of all possibilities is beyond the scope of this paper, we will brielfy consider a simple EW completion of the toy-model discussed in section 9, and comment on the associated phenomenology.
For each EW singlet PQ scalar Φ f , we introduce a new SU(2) W doublet H f with the same corresponding PQ charge, Q PQ H f = Q PQ Φ f . We can then write EW preserving Yukawa couplings, H f f L f c , and tri-scalar "A-terms":

JHEP07(2018)092
Giving large masses m H f O(100 GeV) for the new PQ doublets H f , we can integrate them out and obtain the effective interactions below the EW symmetry breaking scale: L Λ=GeV PQ ⊃ − y u Φ u uu c + y d Φ d dd c + y e Φ e ee c + h.c. (10.2) where with H SM = 174 GeV. A parameter of particular relevance for phenomenology is the mixing angle between the light singlet Φ f and the neutral component of the heavy doublet H 0 f : Firstly, this parameter quantifies the degree of tree-level tuning required to maintain the hierarchy between the singlet states (at the GeV scale) and the doublet states (at the EW scale). A simple measure of this fine-tuning is given by: Such sizable Yukawa coupling would lead to a large production of H u at the LHC, with 13 TeV cross-sections ranging from O(1 − 10 3 ) pb for θ f ΦH ∼ (0.01 − 0.03) and m Hu ∼ (100 − 500) GeV. Existing LHC searches for leptophobic vectors (Z ) [218] can already place non-trivial upper bounds on processes such as pp → H u (→ jj) + j [219][220][221] and pp → H u (→ jj) + γ [221,222], ranging from Y u (0.1 − 0.4) depending on m Hu . Similar considerations hold for H d production.
Finally, θ f ΦH controls a variety of rare decays of the 125 GeV Higgs h, Z 0 and H f to final states with PQ scalars a, η PQ , ϕ u , ϕ d . We compile some of these possibilities in table 2. The reach of existing measurements and future search strategies for these rare decays will depend on how the boosted PQ scalars will be tagged by LHC detection algorithms. Typically, the GeV states ϕ u , ϕ d , η PQ will decay promptly to collimated final states of 2π, ηπ, KK, 3π, ηππ, etc, which may be tagged as jets, hadronic taus, photonjets [223][224][225] (for final states of 2 or 3 π 0 s), or even single photons (if the detector's granularity is poor enough, which might be the case for LEP detectors). The much lighter and longer lived MeV axion a, with cτ a ranging from 0.1 µm to 0.1 mm, will produce a JHEP07(2018)092  Table 2. Compilation of rare decays that could potentially probe the model in (10.1). The middle column shows potential signatures, assuming that the MeV axion a would be tagged either as a converted photon (γ * ) or a lepton-jet (LJ); and that ϕ u,d , η PQ would be tagged either as a jet (j), a photon-jet (P J), or a hadronic tau (τ h ). The third column shows typical ranges for signal cross-sections/branching ratios assuming θ f ΦH ∼ (0.01 − 0.03) and m Hu ∼ (100 − 300) GeV.
highly collimated e + e − pair, with a decay vertex displaced by a few mm to several cm. Here, too, the sensitivity of any particular analysis will depend on whether a is tagged as a converted photon, a prompt or displaced lepton-jet [226][227][228][229][230][231][232][233][234][235][236][237], or whether it will fail quality criteria for standard objects and simply be vetoed. A thorough study of these possibilities is deferred to [198]. 15 It suffices to say, nonetheless, that for θ f ΦH ∼ O(10 −2 ) all existing bounds are satisfied for the processes listed in table 2, regardless of assumptions on boosted-PQ-object-tagging.
If H e in (10.1) is independent of H u,d , additional leptonic signatures may arise at LEP (depending on the mass of H e ) and/or at the LHC (see [198]).
Finally, note that in order to write the flavor diagonal couplings of H u,d to first generation quarks, we implicitly assumed an MFV-type mechanism which generates the CKM flavor structure of the weak interactions of quarks without spoiling the flavor aligment of Yukawa couplings. Concrete realizations of such mechanism are possible, but are beyond the scope of this study.

Discussion
A short-lived, pion-phobic QCD axion with mass of several MeV might still offer a viable solution to the strong CP problem. Constraints that have excluded generic MeV axions can be evaded by coupling the axion exclusively to the first generation of SM fermions. Bounds from K + → π + a, previously believed to be severe, in fact suffer from large hadronic uncertainties and are currently sufficiently ambiguous to experimentally allow portions of the axion parameter space.

JHEP07(2018)092
The extreme pion-phobia needed to avoid exclusion is a realistic possibility in models with a special relation between the light quark masses and PQ charges, namely, In this study we have imposed this relation ad hoc, but it is easy to envision how it might arise dynamically. For instance, in supersymmetric models of flavor, quartics are often proportional to charges squared, and flavon VEVs (and thus fermions Yukawas) would then naturally be inversely proportional to the charges.
The associated phenomenology of these variants is rich and testable. Several axion signatures are similar to the those of visibly decaying dark photons, and can be searched for by ongoing and near-future dark photon experiments. They also offer alternative explanations to a few discrepancies in data usually attributed to dark photons, such as the Beryllium-8 and KTeV (π 0 → e + e − ) anomalies. At the GeV scale, these models predict new states coupled to light hadrons awaiting to be uncovered. Those, along with the MeV axion, may appear in rare decays of the SM Higgs, Z 0 and other BSM states, yielding exotic signatures with thin jets (i.e., "τ h -like"), prompt or displaced lepton jets, and photon jets.