The QCD axion, precisely

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.


Introduction
In the Standard Model, the sum of the QCD topological angle and the common quark mass phase, θ = θ 0 + arg detM q , is experimentally bounded to lie below O(10 −10 ) from the non-observation of the neutron electric dipole moment (EDM) [1,2]. While θ = O(1) would completely change the physics of nuclei, its effects rapidly decouple for smaller values, already becoming irrelevant for θ 10 −1 ÷ 10 −2 . 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 O(1) 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 f a . 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 O(m π f π /f a ). 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 f a .
Presently astrophysical constraints bound f a between few 10 8 GeV (see for e.g. [11]) and few 10 17 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, nonthermal production from the misalignment mechanism can easily generate a suitable abundance of cold axions for values of f a 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 v 2 ∼ 10 −6 (v 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, O(10 −6 ). 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 (T c ∼170 MeV) chiral Lagrangians allow the temperature dependence of the axion potential and its mass to be computed. Around T c 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 10 6 GeV, 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 m a = 5.70(6)(4) µeV 10 12 GeV f a , the coupling to photons 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 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 1/f a and the weak couplings can be written, without loss of generality, as where the second term defines f a , the dual gluon field strengthG µν = 1 2 µνρσ G ρσ , color indices are implicit, and the coupling to the photon field strength F µν is where E/N 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) j µ a,0 = c 0 qq γ µ γ 5 q is a model dependent axial current made of SM matter fields. The axionic pseudo shift-symmetry, a → a + δ, 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 eq. (1) becomes where 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 1/f a 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 M a , 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 M a 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: where · · · is the trace over flavor indices, B 0 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 f π 92 MeV. In order to derive the leading order effective axion potential we need only consider the neutral pion sector. Choosing Q a proportional to the identity we have On the vacuum π 0 gets a vacuum expectation value (VEV) proportional to φ a to minimize the potential, the last cosine in eq. (8) is 1 on the vacuum, and π 0 can be trivially integrated out leaving the axion effective potential As expected the minimum is at a = 0 (thus solving the strong CP problem). Expanding to quadratic order we get the well-known [5] formula for the axion mass 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 O(1) 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 is roughly a factor of 3 smaller than λ (inst) a = −m 2 a /f 2 a , the one extracted from the single cosine potential V inst (a) = −m 2 a f 2 a cos(a/f a ). The six-axion couplings differ in sign as well. The VEV for the neutral pion, π 0 = φ a f π can be shifted away by a non-singlet chiral rotation. Its presence is due to the π 0 -a mass mixing induced by isospin breaking effects in eq. (6), but can be avoided by a different choice for Q a , 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 which is responsible for the mixing. It is then enough to choose 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 O(m −1 a ) and tension given by The function E[q] 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, q ∈ [0, 1], only varies E[q] between E[0] = 1 (cosine-like potential limit) and E[1] = 4−2 √ 2 1.17 (for degenerate quarks). For physical quark masses E[q phys ] 1.12, only 12% off the cosine potential prediction, and σ 9m a f 2 a . 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. (8) at the quadratic order in the pion field. For a = θf a = 0 the pion mass becomes and for θ = π the pion mass is reduced by a factor (m d + m u )/(m d − m u ) √ 3. 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) 1 : 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 10% 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 O(m u /m s ) 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 O(30%) 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.

The mass
The first quantity we compute is the axion mass. As mentioned before at leading order in 1/f a the axion can be treated as an external source. Its mass is thus defined as where Z(θ) is the QCD generating functional in the presence of a theta term and χ top 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 χ top has been computed in [36]. We recomputed this quantity independently and present the result for the axion mass directly in terms of observable renormalized quantities 2 .
The computation is very simple but the result has interesting properties: where h r 1 , h r 3 , l r 4 and l r 7 are the renormalized NLO couplings of [26] and m π and f π 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 m 2 π f 2 π ). In particular l r 7 and the combinations h r 1 − h r 3 − l r 4 are separately scale invariant. Similar properties are also present in the 3-flavor computation, in particular there are no O(m s ) 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 l r 7 could be extracted from the QCD contribution to the π + -π 0 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 h r 1 −h r 3 , 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 The first term in l r 7 is due to the tree-level contribution to the π + -π 0 mass splitting due to the π 0 -η 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 L 7 and L r 8 , 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 O(1) 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 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 z is scale independent up to α em and Yukawa suppressed corrections. Note also that since lattice QCD simulations allow us to relate physical observables directly to the high-energy MS Yukawa couplings, in principle 3 , 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 z will increase further in the near future.
Combining everything together we get the following numerical estimate for the axion mass m a = 5.70(6)(4) µeV 10 12 GeV f a = 5.70 (7) µeV 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 (20). The total error of ∼1% is much smaller than the relative errors in the quark mass ratio (∼6%) and in the NLO couplings (∼30÷60%) because of the weaker dependence of the axion mass on these quantities 3 Modulo well-known effects present when chiral non-preserving fermions are used.
Note that the full NLO correction is numerically smaller than the quark mass error and its uncertainty is dominated by l r 7 . The error on the latter is particularly large because of a partial cancellation between L r 7 and L r 8 in eq. (20). The numerical irrelevance of the other NLO couplings leaves a lot of room for improvement should l r 7 be extracted directly from Lattice QCD. The value of the pion decay constant we used (f π = 92.21 (14) MeV) [43] is extracted from π + decays and includes the leading QED corrections, other O(α em ) corrections to m a 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 χ against which lattice simulations can be calibrated.

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 (m π and f π ). The full result is where m 2 π (θ) is the function defined in eq. (16), and all quantities have been rewritten in terms of the physical NLO quantities 4 . 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 l r 3,4 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 (25) around the origin We find where m a is the physical one-loop corrected axion mass of eq. (19). Numerically we have the error on this quantity amounts to roughly 6% and is dominated by the uncertainty on l r 7 . Finally the NLO result for the domain wall tensions can be simply extracted from the definition using the NLO expression (25) for the axion potential. The numerical result is the error is sub percent and it receives comparable contributions from the errors on l r 7 and the quark masses.
As a by-product we also provide a precision estimate of the topological quartic moment of the topological charge Q top to be compared to the cosine-like potential b inst 2 = −1/12 −0.083.

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 m a suppressed corrections which are negligible, thus the only model dependence lies in the anomaly coefficient E/N . For physical quark masses the QCD contribution (the second term in eq. (17)) is accidentally close to −2. This implies that models with E/N = 2 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 The NLO corrections in the square brackets come from tree-level diagrams with insertions of NLO WZW operators (the terms proportional to thec W i couplings 5 ) and from a-π 0 mixing diagrams (the term proportional to l r 7 ). One loop diagrams exactly cancel similarly to what happens for π → γγ and η → γγ [46]. Notice that the l r 7 term includes the m u /m s 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 couplingsc W i , so we describe the way to extract them here. The first obvious observable we can use is the π 0 → γγ width. Calling δ i the relative correction at NLO to the amplitude for the i process, i.e.
the expressions for Γ tree πγγ and δ πγγ read 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 thec W i couplings enter eq. (33) only at NLO in the light quark mass expansion we only need to determine them at LO in the m u,d expansion.
The η → γγ decay rate at NLO is 5 For simplicity we have rescaled the original couplings c W i of [45] where in the last step we consistently neglected higher order corrections O(m u,d /m s ). The 3-flavor couplingsC W i ≡ (4πf π ) 2 C W i are defined in [45]. The expression for the correction to the π → γγ amplitude with 3 flavors also receives important corrections from the π-η mixing 2 , where the π-η mixing derived in [27] can be conveniently rewritten as at leading order in m u,d . In both decay rates the loop corrections are reabsorbed in the renormalization of the tree-level amplitude 6 . By comparing the light quark mass dependence in eqs. (35) and (37) we can match the 2 and 3 flavor couplings as follows 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 C W 7,8 . The result is shown in fig. 2, the precision is low for two reasons: 1)C W 7,8 are 3 flavor couplings so they suffer from an intrinsic O(30%) uncertainty from higher order corrections 7 , 2) for π → γγ the experimental uncertainty is not smaller than the NLO corrections we want to fit.
For the combination 5c W 3 +c W 7 + 2c W 8 we are interested in, the final result reads When combined with eq. (33) we finally get 6 NLO 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 O(m u,d /m s ), which are however subleading. For the π → γγ rate we also included the mixed term coming from the product of the NLO corrections to 2 and to Γ ηγγ . Formally this term is NNLO but given that the NLO corrections to both 2 and Γ ηγγ are of the same size as the corresponding LO contributions such terms cannot be neglected. 7 We implement these uncertainties by adding a 30% error on the experimental input values of δ πγγ and δ ηγγ . Note that despite the rather large uncertainties of the NLO couplings we are able to extract the model independent contribution to a → γγ 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 E/N to be extracted from a measurement of g aγγ at the percent level.
For the three reference models with respectively E/N = 0 (such as hadronic or KSVZ-like models [6,7] with electrically neutral heavy fermions), E/N = 8/3 (as in DFSZ models [8,9] or KSVZ models with heavy fermions in complete SU(5) representations) and E/N = 2 (as in some KSVZ "unificaxion" models [48]) the coupling reads Even after the inclusion of NLO corrections the coupling to photons in E/N = 2 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 O(20) with respect to E/N = 0 models is highly disfavored. Therefore the result for g E/N =2 aγγ 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.

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 c 0 q 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 / N = 0 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 ∆ ∼ O(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 q 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 where N = (p, n) is the isospin doublet nucleon field, v µ is the four-velocity of the non-relativistic nucleons, D µ = ∂ µ − V µ , V µ is the vector external current, σ i are the Pauli matrices, the index q = ( u+d 2 , s, c, b, t) runs over isoscalar quark combinations, 2N S µ N =N γ µ γ 5 N is the nucleon axial current, M a = cos(Q a a/f a )diag(m u , m d ), and A i µ andÂ q µ are the axial isovector and isoscalar external currents respectively. Neglecting SM gauge bosons, the external currents only depend on the axion field as followŝ where we used the short-hand notation c (u±d)/2 ≡ cu±c d 2 . The couplings c q = c q (Q) computed at the scale Q 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 c q (Q) 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 A µ while there are no linear interactions 8 coming from the non derivative terms contained in M a . In Eq. (43) dots stand for higher order terms involving higher powers of the external sources V µ , A µ , and M a . Among these the leading effects to the axion-nucleon coupling will come from isospin breaking terms O(M a A µ ). 9 These corrections are small O( m d −mu ∆ ), below the uncertainties associated to our determination of the effective coupling g q 0 , 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 m a , so higher order terms are negligibly small O(m a /∆). The price to pay is that the couplings g A and g q 0 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, g A and g q 0 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: g A = ∆u − ∆d , g q 0 = (∆u + ∆d, ∆s, ∆c, ∆b, ∆t) , s µ ∆q ≡ p|qγ µ γ 5 q|p , where |p is a proton state at rest, s µ its spin and we used isospin symmetry to relate proton and neutron matrix elements. Note that the isoscalar matrix elements ∆q inside g q 0 depend on the matching scale Q, such dependence is however canceled once the couplings g q 0 (Q) are multiplied by the corresponding UV couplings c q (Q) inside the isoscalar currentsÂ q µ . Non-singlet combinations such as g A are instead protected by non-anomalous Ward identities 10 . For future convenience we set the matching scale Q = 2 GeV.
We can therefore write the EFT Lagrangian (43) directly in terms of the UV couplings as We are thus left to determine the matrix elements ∆q. The isovector combination can be obtained with high precision from β-decays [43] ∆u − ∆d = g A = 1.2723 (23) , where the tiny neutron-proton mass splitting m n − m p = 1.3 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-x 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.
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 ( computed at the scale Q = 2 GeV. 10 This is only true in renormalization schemes which preserve the Ward identities. 11 Details in the way the numbers in eq. (48) are derived are given in appendix A.
We can now use these inputs in the EFT Lagrangian (46) to extract the corresponding axionnucleon couplings: which are defined in analogy to the couplings to quarks as and are scale invariant (as they are defined in the effective theory below the QCD mass gap). The errors in eq. (50) include the uncertainties from the lattice data and those from higher order corrections in the perturbative RG evolution of the axial current (the latter is only important for the coefficients of c 0 s,c,b,t ). The couplings c 0 q are those appearing in eq. (1) computed at the high scale f a = 10 12 GeV. The effect of varying the matching scale to a different value of f a within the experimentally allowed range is smaller than the theoretical uncertainties.
A few considerations are in order. The theoretical errors quoted here are dominated by the lattice results, which for these matrix elements are still in an early phase and the systematic uncertainties are not fully explored yet. Still the error on the final result is already good (below ten percent), and there is room for a large improvement which is expected in the near future. Note that when the uncertainties decrease sufficiently for results to become sensitive to isospin breaking effects, new couplings will appear in eq. (43). These could in principle be extracted from lattice simulations by studying the explicit quark mass dependence of the matrix element. In this regime the experimental value of the isovector coupling g A cannot be used anymore because of different isospin breaking corrections to charged versus neutral currents.
The numerical values of the couplings we get are not too far off those already in the literature (see e.g. [43]). However, because of the caveats in the relation of the deep inelastic scattering and hyperon data to the relevant matrix elements the uncertainties in those approaches are not under control. On the other hand the lattice uncertainties are expected to improve in the near future, which would further improve the precision of the estimate performed with the technique presented here.
The numerical coefficients in eq. (50) include the effect of running from the high scale f a (here fixed to 10 12 GeV) to the matching scale Q = 2 GeV, which we performed at the NLLO order (more details in appendix B). The running effects are evident from the fact that the couplings to nucleons depend on all quark couplings including charm, bottom and top, even though we took the corresponding spin content to vanish. This effect has been neglected in previous analysis.
Finally it is interesting to observe that there is a cancellation in the model independent part of the axion coupling to the neutron in KSVZ-like models, where c 0 q = 0, the coupling to neutrons is suppressed with respect to the coupling to protons by a factor O(10) at least, in fact this coupling still is compatible with 0. The cancellation can be understood from the fact that, neglecting running and sea quark contributions and the down-quark spin content of the neutron ∆u is approximately ∆u ≈ −2∆d, i.e. the ratio m u /m d is accidentally close to the ratio between the number of up over down valence quarks in the neutron. This cancellation may have important implications on axion detection and astrophysical bounds.
In models with c 0 q = 0 both the couplings to proton and neutron can be large, for example for the DFSZ axion models, where c 0 u,c,t = 1 A cancellation in the coupling to neutrons is still possible for special values of tan β.

The hot axion: finite temperature results
We now turn to discuss the properties of the QCD axion at finite temperature. The temperature dependence of the axion potential and its mass are important in the early Universe because they control the relic abundance of axions today (for a review see e.g. [59]). The most model independent mechanism of axion production in the early universe, the misalignment mechanism [15][16][17], is almost completely determined by the shape of the axion potential at finite temperature and its zero temperature mass. Additionally, extra contributions, such as string and domain walls can also be present if the PQ preserving phase is restored after inflation, and might be the dominant source of dark matter [60][61][62][63][64][65][66]. Their contribution also depends on the finite temperature behavior of the axion potential, although there are larger uncertainties in this case coming from the details of their evolution (for a recent numerical study see e.g. [67]). 12 One may naively think that, as the temperature is raised, our knowledge of axion properties gets better and better-after all the higher the temperature the more perturbative QCD gets. The opposite is instead true. In this section we show that, at the moment, the precision with which we know the axion potential worsens as the temperature is increased! At low temperature this is simple to understand. Our high precision estimates at zero temperature rely on chiral Lagrangians whose convergence degrades as the temperature approaches the critical temperature T c 160-170 MeV where QCD starts deconfining. At T c the chiral approach is already out of control. Fortunately around the QCD cross-over region lattice computations are possible. The current precision is not yet competitive with our low temperature results but they are expected to improve soon. At higher temperatures there are no lattice results available. For T T c the dilute instanton gas approximation, being a perturbative computation, is believed to give a reliable estimate of the axion potential. It is known however that finite temperature QCD converges fast only for very large temperatures, above O(10 6 ) GeV (see e.g. [72]). The situation is particularly bad for the instanton computation. The screening of QCD charge causes an exponential sensitivity to quantum thermal loop effects. The resulting uncertainty on the axion mass and potential can easily be one order of magnitude or more! This is compatible with a recent lattice computation [31], performed without quarks, which found a high temperature axion mass differing from the instanton prediction at T = 1 GeV by a factor ∼ 10. More recent preliminary results from simulations with dynamical quarks [29] seem to show an even bigger disagreement, perhaps suggesting that at these temperatures even the form of the action is very different from the instanton prediction.

Low temperatures
For temperatures T below T c axion properties can reliably be computed within finite temperature chiral Lagrangians [73,74]. Given the QCD mass gap in this regime temperature effects are exponentially suppressed.
The computation of the axion mass is straightforward. Note that the temperature dependence can only come from the non local contributions that can feel the finite temperature. At one loop the axion mass only receives contribution from the local NLO couplings once rewritten in terms of the physical m π and f π [75]. This means that the leading temperature dependence is completely determined by the temperature dependence of m π and f π , and in particular is the same as that of the chiral condensate [73][74][75] where The function J 1 (ξ) asymptotes to ξ 1/4 e − √ ξ /(2π) 3/2 at large ξ and to 1/12 at small ξ. Note that in the ratio m 2 a (T )/m 2 a the dependence on the quark masses and the NLO couplings cancel out. This means that, at T T c , this ratio is known at a even better precision than the axion mass at zero temperature itself.
Higher order corrections are small for all values of T below T c . There are also contributions from the heavier states that are not captured by the low energy Lagrangian. In principle these are exponentially suppressed by e −m/T , where m is the mass of the heavy state. However, because the ratio m/T c is not very large and a large number of states appear above T c there is a large effect at around T c , where the chiral expansion ceases to reliably describe QCD physics. An in depth discussion of such effects appears in [76] for the similar case of the chiral condensate.
The bottom line is that for T T c eq. (55) is a very good approximation for the temperature dependence of the axion mass. At some temperature close to T c eq. (55) suddenly ceases to be a good approximation and full non-perturbative QCD computations are required.
The leading finite temperature dependence of the full potential can easily be derived as well, The temperature dependent axion mass, eq. (55), can also be derived from eq. (57) by taking the second derivative with respect to the axion. The fourth derivative provides the temperature correction to the self-coupling,

High temperatures
While the region around T c is clearly in the non-perturbative regime, for T T c QCD is expected to become perturbative. At large temperatures the axion potential can thus be computed in perturbation theory, around the dilute instanton gas background, as described in [77]. The point is that, at high temperatures large gauge configurations, which would dominate at zero temperature because of the larger gauge coupling, are exponentially suppressed because of Debye screening. This makes the instanton computation a sensible one.
The prediction for the axion potential is of the form V inst (a; the integral is over the instanton size ρ, n(ρ, 0) ∝ m u m d e −8π 2 /g 2 s is the zero temperature instanton density, m 2 D1 = g 2 s T 2 (1+n f /6) is the Debye mass squared at LO, n f is the number of flavor degrees of freedom active at the temperature T , and the dots stand for smaller corrections (see [77] for more details). The functional dependence of eq. (59) on temperature is approximately a power law T −α where α ≈ 7 + n f /3 + . . . is fixed by the QCD beta function.
There is however a serious problem with this type of computation. The dilute instanton gas approximation relies on finite temperature perturbative QCD. The latter really becomes perturbative only at very high temperatures T 10 6 GeV due to IR divergences of the thermal bath [78]. Further, due to the exponential dependence on quantum corrections, the axion mass convergence is even worse than many other observables. In fact the LO estimate of the Debye mass m 2
Both lattice [83] and NLO [79] results give a Debye mass m D 1.5 m D1 where m D1 is the leading perturbative result. Since the Debye mass enters the exponent of eq. (59) higher order effects can easily shift the axion mass at a given temperature by an order of magnitude or more.
Given the failure of perturbation theory in this regime of temperatures even the actual form of eq. (59) may be questioned and the full answer could differ from the semiclassical instanton computation even in the temperature dependence and in the shape of the potential. Because of this, direct computations from non-perturbative methods such as lattice QCD are highly welcome.  Figure 4: The temperature dependent axion mass normalized to the zero temperature value (corresponding to the light quark mass values in each computation). In blue the prediction from chiral Lagrangians. In different shades of red the lattice data from ref. [28] for different lattice volumes, and in shades of green the preliminary lattice data from [29] for different lattice spacings. The dotted grey curve shows the interacting instanton liquid model (IILM) result [84].
Recently several computations of the temperature dependence of the topological susceptibility for pure SU(3) Yang-Mills appeared [30,31]. While computations in this theory cannot be used for the QCD axion 13 , they are useful to test the instanton result. In particular in [31] an explicit comparison was made in the interval of temperatures T /T c ∈ [0.9, 4.0]. The results for the temperature dependence and the quartic derivative of the potential are compatible with those predicted by the instanton approximation, however the overall size of the topological susceptibility was found one order of magnitude bigger. While the size of the discrepancy seem to be compatible with a simple rescaling of the Debye mass, it goes in the opposite direction with respect to the one suggested by higher order effects, preferring a smaller value for m D 0.5m D1 . This fact betrays a deeper modification of eq. (59) than a simple renormalization of m D .
Unfortunately no full studies for real QCD are available yet in the same range of temperatures. Results across the crossover region, for T ∈ [140, 200] MeV, are available in [28], which used light quark masses corresponding to m π 200 MeV. Fig. 4 compares these results with the ChPT ones, with nice agreement around T ∼ 140 MeV. The plot is in terms of the ratio m a (T )/m a , which at low temperatures weakens the quark mass dependence, as manifest in the ChPT computation.
However, at high temperature this may not be true anymore. For example the dilute instanton computation suggests m 2 a (T )/m 2 a ∝ (m u + m d ) ∝ m 2 π , which implies that the slope across the crossover region may be very sensitive to the value of the light quark masses. In future lattice computations it is thus crucial to use physical quark masses, or at least to perform a reliable extrapolation to the physical point.
Additionally, while the volume dependence of the results in [28] seems to be under control, the lattice spacing used was rather coarse (a > 0.125 fm) and furthermore not constant with the temperature. Should the strong dependence on the lattice spacing observed in [31] be also present in full QCD lattice simulations a continuum limit extrapolation would become compulsory.
More recently, new preliminary lattice results appeared in [29], for a wider range of temperatures between 150 and 500 MeV. This analysis was performed with 4 dynamical flavors, including the charm quark, but with heavier light quark masses, corresponding to m π 370 MeV. These results are also shown in fig. 4, and suggest that χ(T ) decreases with temperature much more slowly than in the quarkless case, in clear contradiction to the instanton calculation. The analysis also includes different lattice spacing, showing strong discretization effects. Given the strong dependence on the lattice spacing observed, and the large pion mass employed, a proper analysis of the data is required before a direct comparison with the other results can be performed. In particular, the low temperature lattice points exceed the zero temperature chiral perturbation theory result (given their pion mass), which is presumably a consequence of the finite lattice spacing.
If the results for the temperature slope in [29] are confirmed in the continuum limit and for physical quark masses, it would imply a temperature dependence for the topological susceptibility (χ(T ) ∼ T −2 ) departing strongly from the one predicted by instanton computations. As we will see in the next section this could have dramatic consequences in the computation of the axion relic abundance.
For completeness in fig. 4 we also show the result of [84] obtained from an instanton-inspired model, which is sometimes used as input in the computation of the axion relic abundance. Although the dependence at low temperatures explicitly violates low-energy theorems the behaviour at higher temperature is similar to the lattice data by [28], although with a quite different T c .

Implications for dark matter
The amount of axion dark matter produced in the early Universe and its properties depend on whether PQ symmetry is broken or not after inflation. If the PQ symmetry is broken before inflation (H I f a ) and not restored during reheating (T max f a ), after the Big Bang the axion field is uniformly constant over the observable Universe, a(x) = θ 0 f a . The evolution of the axion field, in particular of its zero mode, is described by the equation of motion where we assumed that the shape of the axion potential is well described by the dilute instanton gas approximation, i.e. cosine like. As the Universe cools, the Hubble parameter decreases while the axion potential increases. When the pull from the latter becomes comparable to the Hubble friction, i.e. m a (T ) ∼ 3H, the axion field starts oscillating with frequency m a . This typically happens at temperatures above T c , around the GeV scale, depending on the value of f a and the temperature dependence of the axion mass. Soon after that the comoving number density n a = m a a 2 becomes an adiabatic invariant and the axion behaves as cold dark matter. Alternatively PQ symmetry may be broken after inflation. In this case, immediately after the breaking the axion field finds itself randomly distributed over the whole range [0, 2πf a ]. Such field configurations include strings which evolve with a complex dynamics, but are known to approach a scaling solution [64]. At temperatures close to T c , when the axion field starts rolling because of the QCD potential, domain walls also form. In phenomenologically viable models, the full field configuration, including strings and domain walls, eventually decays into axions, whose abundance is affected by large uncertainties associated with the evolution and decay of the topological defects. Independently of this evolution there is a misalignment contribution to the dark matter relic density from axion modes with very close to zero momentum. The calculation of this is the same as for the case where inflation happens after PQ breaking, except that the relic density must be averaged over all possible values of θ 0 . While the misalignment contribution gives only a part of the full abundance, it can still be used to give an upper bound to f a in this scenario. The current axion abundance from misalignment, assuming standard cosmological evolution, is given by where Ω γ and T γ are the current photon abundance and temperature respectively and s and n a are the entropy density and the average axion number density computed at any moment in time t sufficiently after the axion starts oscillating such that n a /s is constant. The latter quantity can be obtained by solving eq. (60) and depends on 1) the QCD energy and entropy density around T c , 2) the initial condition for the axion field θ 0 , and 3) the temperature dependence of the axion mass and potential. The first is reasonably well known from perturbative methods and lattice simulations (see e.g. [85,86]). The initial value θ 0 is a free parameter in the first scenario, where the PQ transition happen before inflation-since in this case θ 0 can be chosen in the whole interval [0, 2π] only an upper bound to Ω a can be obtained in this case. In the scenario where the PQ phase is instead restored after inflation n a is obtained by averaging over all θ 0 , which numerically corresponds to choosing 14 θ 0 2.1. Since θ 0 is fixed, Ω a is completely determined as a function of f a in this case. At the moment the biggest uncertainty on the misalignment contribution to Ω a comes from our knowledge of m a (T ). Assuming that m a (T ) can be approximated by the power law around the temperatures where the axion starts oscillating, eq. (60) can easily be integrated numerically. In fig. 5 we plot the values of f a that would reproduce the correct dark matter abundance for different choices of χ(T )/χ(0) and α in the scenario where θ 0 is integrated over. We also show two representative points with parameters (α ≈ 8, χ(1 GeV)/χ(0) ≈ few 10 −7 ) and (α ≈ 2, χ(1 GeV)/χ(0) ≈ 10 −2 ) corresponding respectively to the expected behavior from instanton computations and to the suggested one from the preliminary lattice data in [29]. The figure also shows the corresponding temperature at which the axion starts oscillating, here defined by the condition m a (T ) = 3H(T ). Notice that for large values of α, as predicted by instanton computations, the sensitivity to the overall size of the axion mass at fixed temperature (χ(1 GeV)/χ(0)) is weak. However if the slope of the axion mass with the temperature is much smaller, as suggested by the results in [29], then the corresponding value of f a required to give the correct relic abundance can even be larger by an order of magnitude (note also that in this case the temperature at which the axion starts oscillating would be higher, around 4÷5 GeV). The difference between the two cases could be taken as an estimate of the current uncertainty on this type of computation. More accurate lattice results would be very welcome to assess the actual temperature dependence of the axion mass and potential.
To show the impact of this uncertainty on the viable axion parameter space and the experiments probing it, in fig. 6 we plot the various constraints as a function of the Hubble scale during inflation and the axion decay constant. Limits that depend on the temperature dependence of the axion mass are shown for the instanton and lattice inspired forms (solid and dashed lines respectively), corresponding to the labeled points in fig. 5. On the right side of the plot we also show the values of f a that will be probed by ongoing experiments (solid) and those that could be probed by proposed experiments (dashed empty). Orange colors are used for experiments using the axion coupling to photons, blue for the others. Experiments in the last column (IAXO and ARIADNE) do not rely on the axion being dark matter. The boundary of the allowed axion parameter space is constrained by the CMB limits on tensor modes [87], supernova SN1985 and other astrophysical bounds including black-hole superradiance.
When the PQ preserving phase is not restored after inflation (i.e. when both the Hubble parameter during inflation H I and the maximum temperature after inflation T max are smaller than the PQ scale) the axion abundance can match the observed dark matter one for a large range of values of f a and H I by varying the initial axion value θ 0 . In this case isocurvature bounds [88] (see e.g. [89] for a recent discussion) constrain H I from above. At small f a obtaining the correct relic abundance requires θ 0 to be close to π, where the potential is flat, so the the axion begins oscillating at relatively late times. In the limit θ 0 → π the axion energy density diverges. Given the sensitivity of Ω a to θ 0 in this regime, isocurvatures are enhanced by 1/(π − θ 0 ) and the bound on H I is thus strengthened by a factor π − θ 0 . 15 Meanwhile, the axion decay constant is bounded from above by black-hole superradiance. For smaller values of f a axion misalignment can only explain part of the dark matter abundance. In fig. 6 we show the value of f a required to explain Ω DM when θ 0 = 1 and θ 0 = 0.01 for the two reference values of the axion mass temperature parameters.
If the PQ phase is instead restored after inflation, e.g. for high scale inflation models, θ 0 is not a free parameter anymore. In this case only one value of f a will reproduce the correct dark matter abundance. Given our ignorance about the contributions from topological defect we can use the misalignment computation to give an upper bound on f a . This is shown on the bottom-right side of the plot, again for the two reference models, as before. Contributions from higher-modes and topological defects are likely to make such bound stronger by shifting the forbidden region downwards. Note that while the instanton behavior for the temperature dependence of the axion mass would point to axion masses outside the range which will be probed by ADMX (at least in the current version of the experiment), if the lattice behavior will be confirmed the mass window which will be probed would look much more promising.

Conclusions
We showed that several QCD axion properties, despite being determined by non-perturbative QCD dynamics, can be computed reliably with high accuracy. In particular we computed higher order corrections to the axion mass, its self-coupling, the coupling to photons, the full potential and the domain-wall tension, providing estimates for these quantities with percent accuracy. We also showed how lattice data can be used to extract the axion coupling to matter (nucleons) reliably providing estimates with better than 10% precision. These results are important both experimentally, to assess the actual axion parameter space probed and to design new experiments, and theoretically, since in the case of a discovery they would help determining the underlying theory behind the PQ breaking scale.
We also study the dependence of the axion mass and potential on the temperature, which affects the axion relic abundance today. While at low temperature such information can be extracted accurately using chiral Lagrangians at temperatures close to the QCD crossover and above perturbative methods fail. We also point out that instanton computations, which are believed to become reliable at least when QCD becomes perturbative have serious convergence problems, making them unreliable in the whole region of interest. Recent lattice results seem indeed to suggest large deviations from the instanton estimates. We studied the impact that this uncertainty has on the computation of the axion relic abundance and the constraints on the axion parameter space. More dedicated non-perturbative computations are therefore required to reliably determine the axion relic abundance.

Acknowledgments
This work is supported in part by the ERC Advanced Grant no.267985 (DaMeSyFla).

A Input parameters and conventions
For convenience in table 1 we report the values of the parameters used in this work. When uncertainties are not quoted it means that their effect was negligible and they have not been used.
In the following we discuss in more in details the origin of some of these values.

Quark masses
The value of z = m u /m d has been extracted from the following lattice estimates: which use different techniques, fermion formulations, etc. In [90] the extra preliminary result z = 0.49(1)(1) is also quoted, which agrees with the results above. Some results are still preliminary and the study of systematics may not be complete. Indeed the spread from the central values is somewhat bigger than the quoted uncertainties. Averaging the results above we get z = 0.48 (1). Waiting for more complete results and a more systematic study of all uncertainties we used a more conservative error, z = 0.48 (3), which better captures the spread between the different computations.
Axion properties have a much weaker dependence on the strange quark mass which only enter at higher orders. For definiteness we used the value of the ratio r ≡ 2m s m u + m d = 27.4(1) , from [90].

ChPT low energy constants
For the value of the pion decay constant we used the PDG [43] value: f π = 92.21 (14) MeV , which is free from the leading EM corrections present in the leptonic decays used for the estimates. Following [27] the ratio f η /f π can be related to f K /f π , whose value is very well known, up to higher order corrections. Assuming the usual 30% uncertainty on the SU (3) chiral estimates we get f η /f π = 1.3 (1).
For the NLO low energy couplings we used the usual conventions of [26,27]. As described in the main text we used the matching of the 3 and 2 flavor Lagrangians to estimate the SU(2) couplings from the SU(3) ones. In particular we only need the values of L r 7,8 , which we took as computed at the scaleμ = 770 MeV. The first number has been extracted from the fit in [37] using the constraints for L r 4 in [38]. The second from [38]. A 30% intrinsic uncertainty from higher order 3-flavor corrections has been added. This intrinsic uncertainty is not present for the 2-flavor constants where higher order corrections are much smaller.
In the main text we used the values extracted from 3-flavor simulations in [38]. From the values above and using the matching in [27] between the 2 and the 3 flavor theories we can also extract: l 7 = 7(4) 10 −3 , h r 1 − h r 3 − l r 4 = −0.0048 (14) .