Contents

The QCD axion, precisely

Giovanni Grilli di Cortona, Edward Hardy,

[6pt] Javier Pardo Vega and Giovanni Villadoro

SISSA International School for Advanced Studies and INFN Trieste,

Via Bonomea 265, 34136, Trieste, Italy

Abdus Salam International Centre for Theoretical Physics,

Strada Costiera 11, 34151, Trieste, Italy

We show how several properties of the QCD axion can be extracted at high precision using only first principle QCD computations. By combining NLO results obtained in chiral perturbation theory with recent Lattice QCD results the full axion potential, its mass and the coupling to photons can be reconstructed with percent precision. Axion couplings to nucleons can also be derived reliably, with uncertainties smaller than ten percent. The approach presented here allows the precision to be further improved as uncertainties on the light quark masses and the effective theory couplings are reduced. We also compute the finite temperature dependence of the axion potential and its mass up to the crossover region. For higher temperature we point out the unreliability of the conventional instanton approach and study its impact on the computation of the axion relic abundance.

## 1 Introduction

In the Standard Model, the sum of the QCD topological angle and the common quark mass phase, , is experimentally bounded to lie below (10) from the non-observation of the neutron electric dipole moment (EDM) [1, 2]. While would completely change the physics of nuclei, its effects rapidly decouple for smaller values, already becoming irrelevant for . Therefore, its extremely small value does not seem to be necessary to explain any known large-distance physics. This, together with the fact that other phases in the Yukawa matrices are and that can receive non-decoupling contributions from CP-violating new physics at arbitrarily high scales, begs for a dynamical explanation of its tiny value.

Among the known solutions, the QCD axion [3, 4, 5, 6, 7, 8, 9] is probably the most simple and robust: the SM is augmented with an extra pseudo-goldstone boson, whose only non-derivative coupling is to the QCD topological charge and suppressed by the scale . Such a coupling allows the effects of to be redefined away via a shift of the axion field, whose vacuum expectation value (VEV) is then guaranteed to vanish [10]. It also produces a mass for the axion . Extra model dependent derivative couplings may be present but they do not affect the solution of the strong-CP problem. Both the mass and the couplings of the QCD axion are thus controlled by a single scale .

Presently astrophysical constraints bound between few 10 GeV (see for e.g. [11]) and few  GeV [12, 13, 14]. It has been known for a long time [15, 16, 17] that in most of the available parameter space the axion may explain the observed dark matter of the universe. Indeed, non-thermal production from the misalignment mechanism can easily generate a suitable abundance of cold axions for values of large enough, compatible with those allowed by current bounds. Such a feature is quite model independent and, if confirmed, may give non-trivial constraints on early cosmology.

Finally axion-like particles seem to be a generic feature of string compactification. The simplicity and robustness of the axion solution to the strong-CP problem, the fact that it could easily explain the dark matter abundance of our Universe and the way it naturally fits within string theory make it one of the best motivated particle beyond the Standard Model.

Because of the extremely small couplings allowed by astrophysical bounds, the quest to discover the QCD axion is a very challenging endeavor. The ADMX experiment [18] is expected to become sensitive to a new region of parameter space unconstrained by indirect searches soon. Other experiments are also being planned and several new ideas have recently been proposed to directly probe the QCD axion [19, 20, 21, 22]. To enhance the tiny signal some of these experiments, including ADMX, exploit resonance effects and the fact that, if the axion is dark matter, the line width of the resonance is suppressed by ( being the virial velocity in our galaxy) [23, 24]. Should the axion be discovered by such experiments, its mass would be known with a comparably high precision, . Depending on the experiment different axion couplings may also be extracted with a different accuracy.

Can we exploit such high precision in the axion mass and maybe couplings? What can we learn from such measurements? Will we be able to infer the UV completion of the axion? and its cosmology?

In this paper we try to make a small step towards answering some of these questions. Naively, high precision in QCD axion physics seems hopeless. After all most of its properties, such as its mass, couplings to matter and relic abundance are dominated by non perturbative QCD dynamics. On the contrary, we will show that high precision is within reach. Given its extremely light mass, QCD chiral Lagrangians [25, 26, 27] can be used reliably. Performing a NLO computation we are able to extract the axion mass, self coupling and its full potential at the percent level. The coupling to photons can be extracted with similar precision, as well as the tension of domain walls. As a spin-off we provide estimates of the topological susceptibility and the quartic moment with similar precision and new estimates of some low energy constants.

We also describe a new strategy to extract the couplings to nucleons directly from first principle QCD. At the moment the precision is not yet at the percent level, but there is room for improvement as more lattice QCD results become available.

The computation of the axion potential can easily be extended to finite temperature. In particular, at temperatures below the crossover (170 MeV) chiral Lagrangians allow the temperature dependence of the axion potential and its mass to be computed. Around there is no known reliable perturbative expansion under control and non-perturbative methods, such as lattice QCD [28, 29], are required.

At higher temperatures, when QCD turns perturbative, one may be tempted to use the dilute instanton gas approximation, which is expected to hold at large enough temperatures. We point out however that the bad convergence of the perturbative QCD expansion at finite temperatures makes the standard instanton result completely unreliable for temperatures below , explaining the large discrepancy observed in recent lattice QCD simulations [30, 31]. We conclude with a study of the impact of such uncertainty in the computation of the axion relic abundance, providing updated plots for the allowed axion parameter space.

For convenience we report the main numerical results of the paper here, for the mass

 ma=5.70(6)(4)μeV(1012GeVfa),

the coupling to photons

 gaγγ=αem2πfa[EN−1.92(4)],

the couplings to nucleons (for the hadronic KSVZ model for definiteness)

 cKSVZp=−0.47(3),cKSVZn=−0.02(3),

and for the self quartic coupling and the tension of the domain wall respectively

 λa=−0.346(22)⋅m2af2a,σa=8.97(5)maf2a,

where for the axion mass the first error is from the uncertainties of quark masses while the second is from higher order corrections. As a by-product we also provide a high precision estimate of the topological susceptibility and the quartic moment

 χ1/4top=75.5(5) MeV,b2=−0.029(2).

More complete results, explicit analytic formulae and details about conventions can be found in the text. The impact on the axion abundance computation from different finite temperature behaviors of the axion mass is shown in figs. 5 and 6.

The rest of the paper is organized as follows. In section 2 we first briefly review known leading order results for the axion properties and then present our new computations and numerical estimates for the various properties at zero temperature. In section 3 we give results for the temperature dependence of the axion mass and potential at increasing temperatures and the implications for the axion dark matter abundance. We summarize our conclusions in section 4. Finally, we provide the details about the input parameters used and report extra formulae in the appendices.

## 2 The cool axion: T=0 properties

At energies below the Peccei Quinn (PQ) and the electroweak (EW) breaking scales the axion dependent part of the Lagrangian, at leading order in and the weak couplings can be written, without loss of generality, as

 La=12(∂μa)2+afaαs8πGμν~Gμν+14ag0aγγFμν~Fμν+∂μa2fajμa,0 , (1)

where the second term defines , the dual gluon field strength , color indices are implicit, and the coupling to the photon field strength is

 g0aγγ=αem2πfaEN, (2)

where is the ratio of the Electromagnetic (EM) and the color anomaly (=8/3 for complete SU(5) representations). Finally in the last term of eq. (1) is a model dependent axial current made of SM matter fields. The axionic pseudo shift-symmetry, , has been used to remove the QCD angle.

The only non-derivative coupling to QCD can be conveniently reshuffled by a quark field redefinition. In particular performing a change of field variables on the up and down quarks

 q=(ud)→eiγ5a2faQa([]cud),trQa=1, (3)

eq. (1) becomes

 La=12(∂μa)2+14agaγγFμν~Fμν+∂μa2fajμa−¯qLMaqR+h.c. , (4)

where

 gaγγ=αem2πfa[EN−6tr(QaQ2)], jμa=jμa,0−¯qγμγ5Qaq, (5) Ma=eia2faQaMqeia2faQa,Mq= (mu00md),Q=(2300−13).

The advantage of this basis of axion couplings is twofold. First the axion coupling to the axial current only renormalizes multiplicatively unlike the coupling to the gluon operator, which mixes with the axial current divergence at one-loop. Second the only non-derivative couplings of the axion appear through the quark mass terms.

At leading order in the axion can be treated as an external source, the effects from virtual axions being further suppressed by the tiny coupling. The non derivative couplings to QCD are encoded in the phase dependence of the dressed quark mass matrix , while in the derivative couplings the axion enters as an external axial current. The low energy behaviour of correlators involving such external sources is completely captured by chiral Lagrangians, whose raison d’être is exactly to provide a consistent perturbative expansion for such quantities.

Notice that the choice of field redefinition (3) allowed us to move the non-derivative couplings entirely into the lightest two quarks. In this way we can integrate out all the other quarks and directly work in the 2-flavor effective theory, with capturing the whole axion dependence, at least for observables that do not depend on the derivative couplings.

At the leading order in the chiral expansion all the non-derivative dependence on the axion is thus contained in the pion mass terms:

 Lp2⊃2B0f2π4⟨UM†a+MaU†⟩, (6)

where

 U=eiΠ/fπ,Π=(π0√2π+√2π−−π0), (7)

is the trace over flavor indices, is related to the chiral condensate and determined by the pion mass in term of the quark masses, and the pion decay constant is normalized such that  MeV.

In order to derive the leading order effective axion potential we need only consider the neutral pion sector. Choosing proportional to the identity we have

 V(a,π0) =−B0f2π[mucos(π0fπ−a2fa)+mdcos(π0fπ+a2fa)] =−m2πf2π√1−4mumd(mu+md)2sin2(a2fa) cos(π0fπ−ϕa) (8)

where

 tanϕa≡mu−mdmd+mutan(a2fa). (9)

On the vacuum gets a vacuum expectation value (VEV) proportional to to minimize the potential, the last cosine in eq. (2) is 1 on the vacuum, and can be trivially integrated out leaving the axion effective potential

 V(a)=−m2πf2π√1−4mumd(mu+md)2sin2(a2fa). (10)

As expected the minimum is at (thus solving the strong CP problem). Expanding to quadratic order we get the well-known [5] formula for the axion mass

 m2a=mumd(mu+md)2m2πf2πf2a. (11)

Although the expression for the potential (10) was derived long ago [32], we would like to stress some points often under-emphasized in the literature.

The axion potential (10) is nowhere close to the single cosine suggested by the instanton calculation (see fig. 1).

This is not surprising given that the latter relies on a semiclassical approximation, which is not under control in this regime. Indeed the shape of the potential is different from that of a single cosine, and its dependence on the quark masses is non-analytic, as a consequence of the presence of light Goldstone modes. The axion self coupling, which is extracted from the fourth derivative of the potential

 λa≡∂4V(a)∂a4∣∣∣a=0=−m2u−mumd+m2d(mu+md)2m2af2a, (12)

is roughly a factor of 3 smaller than , the one extracted from the single cosine potential . The six-axion couplings differ in sign as well.

The VEV for the neutral pion, can be shifted away by a non-singlet chiral rotation. Its presence is due to the - mass mixing induced by isospin breaking effects in eq. (6), but can be avoided by a different choice for , which is indeed fixed up to a non-singlet chiral rotation. As noticed in [33], expanding eq. (6) to quadratic order in the fields we find the term

 Lp2⊃2B0fπ4faa⟨Π{Qa,Mq}⟩, (13)

which is responsible for the mixing. It is then enough to choose

 Qa=M−1q⟨M−1q⟩, (14)

to avoid the tree-level mixing between the axion and pions and the VEV for the latter. Such a choice only works at tree level, the mixing reappears at the loop level, but this contribution is small and can be treated as a perturbation.

The non-trivial potential (10) allows for domain wall solutions. These have width and tension given by

 σ=8maf2a E[4mumd(mu+md)2],E[q]≡∫10dy√2(1−y)(1−qy). (15)

The function can be written in terms of elliptic functions but the integral form is more compact. Note that changing the quark masses over the whole possible range, , only varies between (cosine-like potential limit) and (for degenerate quarks). For physical quark masses , only 12% off the cosine potential prediction, and .

In a non vanishing axion field background, such as inside the domain wall or to a much lesser extent in the axion dark matter halo, QCD properties are different than in the vacuum. This can easily be seen expanding eq. (2) at the quadratic order in the pion field. For the pion mass becomes

 m2π(θ)=m2π√1−4mumd(mu+md)2sin2(θ2), (16)

and for the pion mass is reduced by a factor . Even more drastic effects are expected to occur in nuclear physics (see e.g. [34]).

The axion coupling to photons can also be reliably extracted from the chiral Lagrangian. Indeed at leading order it can simply be read out of eqs. (4), (5) and (14)111The result can also be obtained using a different choice of , but in this case the non-vanishing - mixing would require the inclusion of an extra contribution from the coupling.:

 gaγγ=αem2πfa[EN−234md+mumd+mu], (17)

where the first term is the model dependent contribution proportional to the EM anomaly of the PQ symmetry, while the second is the model independent one coming from the minimal coupling to QCD at the non-perturbative level.

The other axion couplings to matter are either more model dependent (as the derivative couplings) or theoretically more challenging to study (as the coupling to EDM operators), or both. In section 2.4, we present a new strategy to extract the axion couplings to nucleons using experimental data and lattice QCD simulations. Unlike previous studies our analysis is based only on first principle QCD computations. While the precision is not as good as for the coupling to photons, the uncertainties are already below and may improve as more lattice simulations are performed.

Results with the 3-flavor chiral Lagrangian are often found in the literature. In the 2-flavor Lagrangian the extra contributions from the strange quark are contained inside the low-energy couplings. Within the 2-flavor effective theory the difference between using 2 or 3 flavor formulae, is a higher order effect. Indeed the difference is which corresponds to the expansion parameter of the 2-flavor Lagrangian. As we will see in the next section these effects can only be consistently considered after including the full NLO correction.

At this point the natural question is, how good are the estimates obtained so far using leading order chiral Lagrangians? In the 3-flavor chiral Lagrangian NLO corrections are typically around 20-30%. The 2-flavor theory enjoys a much better perturbative expansion given the larger hierarchy between pions and the other mass thresholds. To get a quantitative answer the only option is to perform a complete NLO computation. Given the better behaviour of the 2-flavor expansion we perform all our computation with the strange quark integrated out. The price we pay is the reduced number of physical observables that can be used to extract the higher order couplings. When needed we will use the 3-flavor theory to extract the values of the 2-flavor ones. This will produce intrinsic uncertainties in the extraction of the 2-flavor couplings. Such uncertainties however will only have a small impact on the final result whose dependence on the higher order 2-flavor couplings is suppressed by the light quark masses.

### 2.1 The mass

The first quantity we compute is the axion mass. As mentioned before at leading order in the axion can be treated as an external source. Its mass is thus defined as

where is the QCD generating functional in the presence of a theta term and is the topological susceptibility.

A partial computation of the axion mass at one loop was first attempted in [35]. More recently the full NLO corrections to has been computed in [36]. We recomputed this quantity independently and present the result for the axion mass directly in terms of observable renormalized quantities222The results in [36] are instead presented in terms of the unphysical masses and couplings in the chiral limit. Retaining the full explicit dependence on the quark masses those formula are more suitable for lattice simulations..

The computation is very simple but the result has interesting properties:

 (19)

where , , and are the renormalized NLO couplings of [26] and and are the physical (neutral) pion mass and decay constant (which include NLO corrections). There is no contribution from loop diagrams at this order (this is true only after having reabsorbed the one loop corrections of the tree-level factor ). In particular and the combinations are separately scale invariant. Similar properties are also present in the 3-flavor computation, in particular there are no corrections (after renormalization of the tree-level result), as noticed already in [35].

To get a numerical estimate of the axion mass and the size of the corrections we need the values of the NLO couplings. In principle could be extracted from the QCD contribution to the - mass splitting. While lattice simulations have started to become sensitive to EM and isospin breaking effects, at the moment there are no reliable estimates of this quantity from first principle QCD. Even less is known about , which does not enter other measured observables. The only hope would be to use lattice QCD computation to extract such coupling by studying the quark mass dependence of observables such as the topological susceptibility. Since these studies are not yet available we employ a small trick: we use the relations in [27] between the 2- and 3-flavor couplings to circumvent the problem. In particular we have

 lr7 =mu+mdmsf2π8m2π−36L7−12Lr8+log(m2η/μ2)+164π2+3log(m2K/μ2)128π2=7(4)⋅10−3, hr1−hr3−lr4 =−8Lr8+log(m2η/μ2)96π2+log(m2K/μ2)+164π2=(4.8±1.4)⋅10−3. (20)

The first term in is due to the tree-level contribution to the - mass splitting due to the - mixing from isospin breaking effects. The rest of the contribution, formally NLO, includes the effect of the - mixing and numerically is as important as the tree-level piece [27]. We thus only need the values of the 3-flavor couplings and , which can be extracted from chiral fits [37] and lattice QCD [38], we refer to appendix A for more details on the values used. An important point is that by using 3-flavor couplings the precision of the estimates of the 2-flavor ones will be limited to the convergence of the 3-flavor Lagrangian. However, given the small size of such corrections even an uncertainty will still translate into a small overall error.

The final numerical ingredient needed is the actual up and down quark masses, in particular their ratio. Since this quantity already appears in the tree level formula of the axion mass we need a precise estimate for it, however, because of the Kaplan-Manohar (KM) ambiguity [39], it cannot be extracted within the meson Lagrangian. Fortunately recent lattice QCD simulations have dramatically improved our knowledge of this quantity. Considering the latest results we take

 z≡m¯¯¯¯¯¯¯MSu(2 GeV)m¯¯¯¯¯¯¯MSd(2 GeV)=0.48(3), (21)

where we have conservatively taken a larger error than the one coming from simply averaging the results in [40, 41, 42] (see the appendix A for more details). Note that is scale independent up to and Yukawa suppressed corrections. Note also that since lattice QCD simulations allow us to relate physical observables directly to the high-energy Yukawa couplings, in principle333Modulo well-known effects present when chiral non-preserving fermions are used., they do not suffer from the KM ambiguity, which is a feature of chiral Lagrangians. It is reasonable to expect that the precision on the ratio will increase further in the near future.

Combining everything together we get the following numerical estimate for the axion mass

 ma =5.70(6)(4) μeV(1012GeVfa)=5.70(7) μeV(1012GeVfa), (22)

where the first error comes from the up-down quark mass ratio uncertainties (21) while the second comes from the uncertainties in the low energy constants (2.1). The total error of 1% is much smaller than the relative errors in the quark mass ratio (6%) and in the NLO couplings (3060%) because of the weaker dependence of the axion mass on these quantities

 ma=[5.70+0.06z−0.480.03−0.04103lr7−74+0.017103(hr1−hr3−lr4)−4.81.4]μeV1012GeVfa. (23)

Note that the full NLO correction is numerically smaller than the quark mass error and its uncertainty is dominated by . The error on the latter is particularly large because of a partial cancellation between and in eq. (2.1). The numerical irrelevance of the other NLO couplings leaves a lot of room for improvement should be extracted directly from Lattice QCD.

The value of the pion decay constant we used ( MeV) [43] is extracted from  decays and includes the leading QED corrections, other corrections to are expected to be sub-percent. Further reduction of the error on the axion mass may require a dedicated study of this source of uncertainty as well.

As a by-product we also provide a comparably high precision estimate of the topological susceptibility itself

 χ1/4top=√mafa=75.5(5)MeV, (24)

against which lattice simulations can be calibrated.

### 2.2 The potential: self-coupling and domain-wall tension

Analogously to the mass, the full axion potential can be straightforwardly computed at NLO. There are three contributions: the pure Coleman-Weinberg 1-loop potential from pion loops, the tree-level contribution from the NLO Lagrangian, and the corrections from the renormalization of the tree-level result, when rewritten in terms of physical quantities ( and ). The full result is

 V(a)NLO= −m2π(afa)f2π{1−2m2πf2π[lr3+lr4−(md−mu)2(md+mu)2lr7−364π2log(m2πμ2)] + (25)

where is the function defined in eq. (16), and all quantities have been rewritten in terms of the physical NLO quantities444See also [44] for a related result computed in terms of the LO quantities.. In particular the first line comes from the NLO corrections of the tree-level potential while the second line is the pure NLO correction to the effective potential.

The dependence on the axion is highly non-trivial, however the NLO corrections account for only up to few percent change in the shape of the potential (for example the difference in vacuum energy between the minimum and the maximum of the potential changes by 3.5% when NLO corrections are included). The numerical values for the additional low-energy constants are reported in appendix A. We thus know the full QCD axion potential at the percent level!

It is now easy to extract the self-coupling of the axion at NLO by expanding the effective potential (2.2) around the origin

 V(a)=V0+12m2aa2+λa4!a4+… (26)

We find

 λa= −m2af2a{m2u−mumd+m2d(mu+md)2 (27) +6m2πf2πmumd(mu+md)2[hr1−hr3−lr4+4¯l4−¯l3−364π2−4m2u−mumd+m2d(mu+md)2lr7]}, (28)

where is the physical one-loop corrected axion mass of eq. (19). Numerically we have

 λa=−0.346(22)⋅m2af2a, (29)

the error on this quantity amounts to roughly 6% and is dominated by the uncertainty on .

Finally the NLO result for the domain wall tensions can be simply extracted from the definition

 σ=2fa∫π0dθ√2[V(θ)−V(0)], (30)

using the NLO expression (2.2) for the axion potential. The numerical result is

 σ=8.97(5)maf2a, (31)

the error is sub percent and it receives comparable contributions from the errors on and the quark masses.

As a by-product we also provide a precision estimate of the topological quartic moment of the topological charge

 b2≡−⟨Q4top⟩−3⟨Q2top⟩212⟨Q2top⟩=f2aV′′′′(0)12V′′(0)=λaf2a12m2a=−0.029(2), (32)

to be compared to the cosine-like potential .

### 2.3 Coupling to photons

Similarly to the axion potential, the coupling to photons (17) also gets QCD corrections at NLO, which are completely model independent. Indeed derivative couplings only produce suppressed corrections which are negligible, thus the only model dependence lies in the anomaly coefficient .

For physical quark masses the QCD contribution (the second term in eq. (17)) is accidentally close to . This implies that models with can have anomalously small coupling to photons, relaxing astrophysical bounds. The degree of this cancellation is very sensitive to the uncertainties from the quark mass and the higher order corrections, which we compute here for the first time.

At NLO new couplings appear from higher-dimensional operators correcting the WZW Lagrangian. Using the basis of [45], the result reads

 gaγγ=αem2πfa{EN−234md+mumd+mu+m2πf2π8mumd(mu+md)2[89(5~cW3+~cW7+2~cW8)−md−mumd+mulr7]}. (33)

The NLO corrections in the square brackets come from tree-level diagrams with insertions of NLO WZW operators (the terms proportional to the couplings555For simplicity we have rescaled the original couplings of [45] into .) and from - mixing diagrams (the term proportional to ). One loop diagrams exactly cancel similarly to what happens for and [46]. Notice that the term includes the contributions which one obtains from the 3-flavor tree-level computation.

Unlike the NLO couplings entering the axion mass and potential little is known about the couplings , so we describe the way to extract them here.

The first obvious observable we can use is the width. Calling the relative correction at NLO to the amplitude for the process, i.e.

 ΓNLOi≡Γtreei(1+δi)2, (34)

 Γtreeπγγ =α2em(4π)3m3πf2π, δπγγ=169m2πf2π[md−mumd+mu(5~cW3+~cW7+2~cW8)−3(~cW3+~cW7+~cW114)]. (35)

Once again the loop corrections are reabsorbed by the renormalization of the tree-level parameters and the only contributions come from the NLO WZW terms. While the isospin breaking correction involves exactly the same combination of couplings entering the axion width, the isospin preserving one does not. This means that we cannot extract the required NLO couplings from the pion width alone. However in the absence of large cancellations between the isospin breaking and the isospin preserving contributions we can use the experimental value for the pion decay rate to estimate the order of magnitude of the corresponding corrections to the axion case. Given the small difference between the experimental and the tree-level prediction for the NLO axion correction is expected of order few percent.

To obtain numerical values for the unknown couplings we can try to use the 3-flavor theory, in analogy with the axion mass computation. In fact at NLO in the 3-flavor theory the decay rates and only depend on two low-energy couplings that can thus be determined. Matching these couplings to the 2-flavor theory ones we are able to extract the required combination entering in the axion coupling. Because the couplings enter eq. (33) only at NLO in the light quark mass expansion we only need to determine them at LO in the expansion.

The decay rate at NLO is

 Γtreeη→γγ =α2em3(4π)3m3ηf2η, δ(3)ηγγ =329m2πf2π[2ms−4mu−mdmu+md~CW7+62ms−mu−mdmu+md~CW8]≃649m2Kf2π(~CW7+6~CW8), (36)

where in the last step we consistently neglected higher order corrections . The 3-flavor couplings are defined in [45]. The expression for the correction to the amplitude with 3 flavors also receives important corrections from the - mixing ,

 δ(3)πγγ =329m2πf2π[md−4mumu+md~CW7+6md−mumu+md~CW8]+fπfηϵ2√3(1+δηγγ), (37)

where the - mixing derived in [27] can be conveniently rewritten as

 ϵ2√3≃md−mu6ms[1+4m2Kf2π(lr7−164π2)], (38)

at leading order in . In both decay rates the loop corrections are reabsorbed in the renormalization of the tree-level amplitude666NLO corrections to and decay rates to photons including isospin breaking effects were also computed in [47]. For the rate we disagree in the expression of the terms , which are however subleading. For the rate we also included the mixed term coming from the product of the NLO corrections to and to . Formally this term is NNLO but given that the NLO corrections to both and are of the same size as the corresponding LO contributions such terms cannot be neglected. .

By comparing the light quark mass dependence in eqs. (35) and (37) we can match the 2 and 3 flavor couplings as follows

 ~cW3+~cW7+~cW114 =~CW7, 5~cW3+~cW7+2~cW8 =5~CW7+12~CW8+332f2πm2K[1+4m2Kfπfη(lr7−164π2)](1+δηγγ). (39)

Notice that the second combination of couplings is exactly the one needed for the axion-photon coupling. By using the experimental results for the decay rates (reported in appendix A), we can extract . The result is shown in fig. 2, the precision is low for two reasons: 1) are 3 flavor couplings so they suffer from an intrinsic (30%) uncertainty from higher order corrections777We implement these uncertainties by adding a 30% error on the experimental input values of and ., 2) for the experimental uncertainty is not smaller than the NLO corrections we want to fit.

For the combination we are interested in, the final result reads

 5~cW3+~cW7+2~cW8 =3f2π64m2Kmu+mdmu{[1+4m2Kf2π(lr7−164π2)]fπfη(1+δηγγ)+3δηγγ−6m2Km2πδπγγ} =0.033(6). (40)

When combined with eq. (33) we finally get

 gaγγ=αem2πfa[EN−1.92(4)]=[0.203(3)EN−0.39(1)]maGeV2. (41)

Note that despite the rather large uncertainties of the NLO couplings we are able to extract the model independent contribution to at the percent level. This is due to the fact that, analogously to the computation of the axion mass, the NLO corrections are suppressed by the light quark mass values. Modulo experimental uncertainties eq. (41) would allow the parameter to be extracted from a measurement of at the percent level.

For the three reference models with respectively (such as hadronic or KSVZ-like models [6, 7] with electrically neutral heavy fermions), (as in DFSZ models [8, 9] or KSVZ models with heavy fermions in complete SU(5) representations) and (as in some KSVZ “unificaxion” models [48]) the coupling reads

 gaγγ=⎧⎪ ⎪⎨⎪ ⎪⎩−2.227(44)⋅10−3/faE/N=0−0.870(44)⋅10−3/faE/N=8/3−0.095(44)⋅10−3/faE/N=2. (42)

Even after the inclusion of NLO corrections the coupling to photons in models is still suppressed. The current uncertainties are not yet small enough to completely rule out a higher degree of cancellation, but a suppression bigger than (20) with respect to models is highly disfavored. Therefore the result for of eq. (42) can now be taken as a lower bound to the axion coupling to photons, below which tuning is required. The result is shown in fig. 3.

### 2.4 Coupling to matter

Axion couplings to matter are more model dependent as they depend on all the UV couplings defining the effective axial current (the constants in the last term of eq. (1)). In particular, there is a model independent contribution coming from the axion coupling to gluons (and to a lesser extent to the other gauge bosons) and a model dependent part contained in the fermionic axial couplings.

The couplings to leptons can be read off directly from the UV Lagrangian up to the one loop effects coming from the coupling to the EW gauge bosons. The couplings to hadrons are more delicate because they involve matching hadronic to elementary quark physics. Phenomenologically the most interesting ones are the axion couplings to nucleons, which could in principle be tested from long range force experiments, or from dark-matter direct-detection like experiments.

In principle we could attempt to follow a similar procedure to the one used in the previous section, namely to employ chiral Lagrangians with baryons and use known experimental data to extract the necessary low energy couplings. Unfortunately effective Lagrangians involving baryons are on much less solid ground—there are no parametrically large energy gaps in the hadronic spectrum to justify the use of low energy expansions.

A much safer thing to do is to use an effective theory valid at energies much lower than the QCD mass gaps (100 MeV). In this regime nucleons are non-relativistic, their number is conserved and they can be treated as external fermionic currents. For exchanged momenta parametrically smaller than , heavier modes are not excited and the effective field theory is under control. The axion, as well as the electro-weak gauge bosons, enters as classical sources in the effective Lagrangian, which would otherwise be a free non-relativistic Lagrangian at leading order. At energies much smaller than the QCD mass gap the only active flavor symmetry we can use is isospin, which is explicitly broken only by the small quark masses (and QED effects). The leading order effective Lagrangian for the 1-nucleon sector reads

 LN=¯NvμDμN+2gAAiμ¯NSμσiN+2gq0^Aqμ¯NSμN+σ⟨Ma⟩¯NN+b¯NMaN+… (43)

where is the isospin doublet nucleon field, is the four-velocity of the non-relativistic nucleons, , is the vector external current, are the Pauli matrices, the index runs over isoscalar quark combinations, is the nucleon axial current, , and and are the axial isovector and isoscalar external currents respectively. Neglecting SM gauge bosons, the external currents only depend on the axion field as follows

 ^Aqμ=cq∂μa2fa,A3μ=c(u−d)/2 ∂μa2fa,A1,2μ=Vμ=0, (44)

where we used the short-hand notation . The couplings computed at the scale will in general differ from the high scale ones because of the running of the anomalous axial current [49]. In particular under RG evolution the couplings mix, so that in general they will all be different from zero at low energy. We explain the details of this effect in appendix B.

Note that the linear axion couplings to nucleons are all contained in the derivative interactions through while there are no linear interactions888This is no longer true in the presence of extra CP violating operators such as those coming from the CKM phase or new physics. The former are known to be very small, while the latter are more model dependent and we will not discuss them in the current work. coming from the non derivative terms contained in . In Eq. (43) dots stand for higher order terms involving higher powers of the external sources , , and . Among these the leading effects to the axion-nucleon coupling will come from isospin breaking terms .999Axion couplings to EDM operators also appear at this order. These corrections are small , below the uncertainties associated to our determination of the effective coupling , which are extracted from lattice simulations performed in the isospin limit.

Eq. (43) should not be confused with the usual heavy baryon chiral Lagrangian [50] because here pions have been integrated out. The advantage of using this Lagrangian is clear: for axion physics the relevant scale is of order , so higher order terms are negligibly small . The price to pay is that the couplings and can only be extracted from very low-energy experiments or lattice QCD simulations. Fortunately the combination of the two will be enough for our purposes.

In fact, at the leading order in the isospin breaking expansion, and can simply be extracted by matching single nucleon matrix elements computed with the QCD+axion Lagrangian (4) and with the effective axion-nucleon theory (43). The result is simply:

 gA=Δu −Δd,gq0=(Δu+Δd,Δs,Δc,Δb,Δt),sμΔq≡⟨p|¯qγμγ5q|p⟩, (45)

where is a proton state at rest, its spin and we used isospin symmetry to relate proton and neutron matrix elements. Note that the isoscalar matrix elements inside depend on the matching scale , such dependence is however canceled once the couplings are multiplied by the corresponding UV couplings inside the isoscalar currents . Non-singlet combinations such as are instead protected by non-anomalous Ward identities101010This is only true in renormalization schemes which preserve the Ward identities.. For future convenience we set the matching scale  GeV.

We can therefore write the EFT Lagrangian (43) directly in terms of the UV couplings as

 LN=¯NvμDμN+∂μafa{ cu−cd2(Δu−Δd)¯NSμσ3N +[cu+cd2(Δu+Δd)+∑q=s,c,b,tcqΔq]¯NSμN}. (46)

We are thus left to determine the matrix elements . The isovector combination can be obtained with high precision from -decays [43]

 Δu−Δd=gA=1.2723(23), (47)

where the tiny neutron-proton mass splitting  MeV guarantees that we are within the regime of our effective theory. The error quoted is experimental and does not include possible isospin breaking corrections.

Unfortunately we do not have other low energy experimental inputs to determine the remaining matrix elements. Until now such information has been extracted from a combination of deep-inelastic-scattering data and semi-leptonic hyperon decays: The former suffer from uncertainties coming from the integration over the low- kinematic region, which is known to give large contributions to the observable of interest; the latter are not really within the EFT regime, which does not allow a reliable estimate of the accuracy.

Fortunately lattice simulations have recently started producing direct reliable results for these matrix elements. From [51, 52, 53, 54, 55, 56] (see also [57, 58]) we extract111111Details in the way the numbers in eq. (48) are derived are given in appendix A. the following inputs computed at  GeV in

 gud0=Δu+Δd=0.521(53),Δs=−0.026(4),Δc=±0.004. (48)

Notice that the charm spin content is so small that its value has not been determined yet, only an upper bound exists. Similarly we can neglect the analogous contributions from bottom and top quarks which are expected to be even smaller. As mentioned before, lattice simulations do not include isospin breaking effects, these are however expected to be smaller than the current uncertainties. Combining eqs. (47) and (48) we thus get:

 Δu=0.897(27),Δd=−0.376(27),Δs=−0.026(4), (49)

computed at the scale  GeV.

We can now use these inputs in the EFT Lagrangian (2.4) to extract the corresponding axion-nucleon couplings:

 cp =−0.47(3)+0.88(3)c0u−0.39(2)c0d−0.038(5)c0s −0.012(5)c0c−0.009(2)c0b−0.0035(4)c0t,