Precision calculation of the axion-nucleon coupling in chiral perturbation theory

We derive the axion-nucleon interaction Lagrangian in heavy baryon chiral perturbation theory up to next-to-next-to-leading order. The effective axion-nucleon coupling is calculated to a few percent accuracy.


I. INTRODUCTION
More than forty years after the proposal to add another symmetry in QCD, viz. the Peccei-Quinn (PQ) symmetry U(1) PQ [1, 2], the pseudo-Nambu-Goldstone boson resulting from the spontaneous breakdown of this symmetry, the QCD axion, remains one of the most favored candidates for a Beyond the Standard Model (BSM) particle. The reasons are manifold: Originally introduced as a resolution of the so-called strong-CP problem, i.e. the question why the observable QCD vacuum angleθ = θ QCD + arg det M q (with M q the quark mass matrix) is such a small quantity (current measurements of the neutron electric dipole moment imply |θ| 10 −11 [3][4][5]), its experimental detection would not only unequivocally solve the strong-CP problem, but potentially also provide an answer (or complement the answer) on the question of the nature of the cosmological dark matter, another pressing issue in contemporary physics research. At the same time a model with PQ symmetry breaking can lead to massive Majorana or Dirac neutrinos depending on the choice of assigning PQ charges to the SM particles and Higgses [6][7][8][9][10][11][12][13][14][15]. Moreover, the fact that axions and PQ symmetries arise quite naturally in superstring theory [16] increases their popularity further.
If the QCD axion indeed exists, its couplings to Standard Model particles, i.e. matter particles and gauge bosons, and hence to composite particles as nucleons, must be very weak, because these are controlled by the very large axion decay constant f a . Currently only lower and upper bounds on these couplings can be given. If these bounds are determined from nuclear processes, the exactness of the determination of f a then strongly depends on the accuracy of our knowledge on the effective axion-nucleon coupling strength.
Our work is organized as follows: In Sec. II, we recapitulate the interaction Lagrangian between quarks and the axion. Then, in Sec. III, we derive the axion-nucleon interaction to the third order, that is including all terms up to next-to-next-to-leading order and pion loop contributions. We also give the numerical values of the axion coupling to neutrons and protons. We end with a short summary in Sec. IV.

II. AXION-QUARK INTERACTION LAGRANGIAN
Consider the QCD Lagrangian including the axion field a(x) at energies below the PQ where L QCD,0 contains all terms that are not of interest in what follows, including the axionphoton interaction term [28,37,38]. Furthermore, the quark mass matrix, f a is the axion decay constant, g the strong interaction coupling constant, G µν = G a µν λ a /2 is the conventional gluon field strength tensor with λ a the Gell-Mann matrices, andG µν = 1 2 ǫ µναβ G αβ its dual, where the trace hence acts in the color space. The PQ current is given by from which the first term gives rise to the kinetic term of the axion, whereas the second term describes the axion-quark interactions proportional to the model-dependent coupling constants combined in the matrix X q = diag(X q ) acting in the flavor space. These are given by for the KSVZ axion and the DFSZ axion, respectively, and x = cot β is the ratio of the vacuum expectation values (VEVs) of the two Higgs doublets. We exclude the PQWW axion from the analysis since it has been ruled out experimentally [17,38].
Note that we do not integrate out the heavy quarks from the beginning. As, depending on the model, the axion couplings to heavy quarks are quite a possibility, they might contribute to the axion-nucleon interactions due to sea quark effects. Additionally, as the couplings X q are scale-dependent quantities [28], running effects enter the couplings to nucleons, which can only be recovered if the axion interactions with heavy quarks are taken along within the calculations, to wit in the form of isoscalar currents.
It is advisable to perform an axial rotation on the quark fields in order to remove the term ∝ a Tr G µνG µν in Eq. (1) by transforming with where z = m u /m d and w = m u /m s . The particular form of Q a has been chosen in order to avoid the π 0 -a mass mixing [22].
With that transformation, the Lagrangian can be written as where now the non-derivative axion-quark interactions are entirely shifted into the phase of the mass matrix, whereas the derivative axion-quark interactions are present in the coupling to the axion which is now anomaly-free [17,39]. The terms in J a µ must be split into isoscalar and isovector pieces in order to translate it later into an effective field theory (EFT) language. Consider the two-dimensional subspace of Eq. (8) with q = (u, d) T : where τ 3 is the conventional third Pauli matrix and we have introduced the abbreviations Setting furthermore and inserting Eq. (9) into Eq. (6), one finds the general axion-quark interaction Lagrangian which is now expressed in a suitable basis so that the isovector and isoscalar parts of the axion-nucleon interaction can easily be extracted.

III.1. The Lagrangian
We construct the HBCHPT Lagrangian with the additional axion field and its interactions by adapting the one developed in Ref. [40] (including the notation) and adding additional terms allowed from symmetries containing the isoscalar axial currents. Usually, the axial currents entering the HBCHPT Lagrangian as external sources are taken to be traceless in order to avoid subtleties arising from the U(1) A anomaly. However, here our model is anomaly-free by construction and we now have to add the isoscalar axial currents appearing separately in the Lagrangian in Eq. (12) that are not traceless. This is done in complete analogy to the traceless axial currents.
We introduce which contains the three pseudo-Nambu-Goldstone bosons of the spontaneously broken chiral symmetry, with the index a = 1, 2, 3 and summation implied. Furthermore, F π is the pion decay constant in the chiral limit, for which we will take the physical value 92.1 MeV for the difference to the chiral limit value only amounts to effects of higher orders than those considered here. The isovector axial current a µ enters the theory by means of the chiral connection Γ µ of the covariant derivative, the so-called vielbein, and the field-strength tensor, where we have set v µ = 0 for the external vector field. Note that the field-strength tensor vanishes in the present model because a µ ∝ ∂ µ a τ 3 as can be read off from the Lagrangian (12).
Introducing thus the isoscalar axial current a s µ,i ∝ c i ∂ µ a ½, we can construct similar objects: which vanishes due to a s µ,i u = ua s µ,i , and a vielbein equivalent The corresponding field strength tensors, of course, vanish as in the case of the isovector axial current. The index i = (u + d, s, c, b, t) runs over all isoscalar quark combinations, cf.
Eq. (12). Furthermore, we need as the last building block, where χ = 2B(s − ip) includes the external scalar and pseudoscalar fields s(x) and p(x), and B is a constant related to quark condensate Σ = − ūu via B = lim mu,m d →0 (Σ/F 2 ).
Collecting the proton and neutron fields in the isodoublet N(x) = (p, n) T , the most general HBCHPT Lagrangian up to order O (p 3 ) in the low-energy expansion, reads where we only show terms that finally lead to interaction vertices with only one single axion, because interactions containing n axions are suppressed by factors 1/f n a and can hence be neglected. Note that there is no isoscalar counterpart to the d 17 term proportional to a low-energy constant (LEC) d i 17 , since such a term would have the same structure as the d i 16 term, becauseũ i ∝ ½, and thus not independent. The terms in the first line are the leading order and next-to-leading order terms, while all the other terms are the next-to-next-toleading order terms. In the Lagrangian (21), v µ is the nucleon four-velocity and N = N v are velocity-dependent nucleon fields with mass m. Strictly speaking, m is the nucleon mass in the two-flavor chiral limit, often denoted asm N . We will suppress the index v in what follows. The axial couplings g A and g i 0 should also be taken in the chiral limit, but we will later match them with the nucleon matrix elements ∆q, which refer to the physical values of the quark masses, see below. S µ is the covariant spin-operator, which has the following properties in d dimensions needed later, employing dimensional regularization to deal with the appearing divergences: Besides the parameters already mentioned, a number of new LECs d where d In these equations, L(λ) contains the divergence at space-time dimension d = 4, and the β-functions are set to cancel the divergences of the one-loop functional, as discussed below.
In order to derive the full axion-nucleon coupling at O(p 3 ), the Lagrangian (21)  All interaction terms of the Lagrangian (21) contribute. The expressions for the external sources can be read off from the axion-quark interaction Lagrangian (12): Setting hence u = ½ and expanding the exponential in M a , cf. eq. (7), up to O(a 1 ), we find where we have inserted the leading order pion mass M 2 π = B(m u + m d ). Introducing the abbreviation the single axion-nucleon interaction Lagrangian reads From that, we can derive the corresponding tree-level NNa-vertex Feynman rule, where q is the momentum of the outgoing axion, p the momentum of the incoming nucleon, and p ′ = p − q the momentum of the outgoing nucleon. Furthermore, we have set ω (′) = v · p (′) . Note that this expression (31) contains divergences due to the terms ∝ d (i) 16 (λ) and ∝d (i) 25,29 (λ).

III.3. Pion-loop contributions
According to the usual power counting scheme, see, e. g., Ref.
we find and thus This yields the following Feynman rules with k the outgoing pion momentum, and q the outgoing axion momentum: In what follows, we will make use of the loop functions given in the appendix of Ref. [41]: where L(λ) is given in eq. (26), and where in particular the expression of J 0 (ω) is valid for ω < M π , which is the region we are interested in.
Using the Feynman rules given above, the first loop diagram is calculated as where we have set (summation over i, j implied) Using the identity the equation (41) can be written as The terms ∝ J 3 (ω (′) ) vanish because of Eq. (23). The final result is found using the anticommutator (23) and inserting the loop function (39), where we have separated the finite, scale-independent terms (the first and second lines) from the divergent or scale-dependent ones (the third line).

III.4. Axion-nucleon coupling at O(p 3 )
In order to remove the divergences appearing in the coupling (31) and diagrams (a) and (c), we utilize the following set of β-functions: where β 16 , β 25 , and β 29 have been calculated already before within the theory without axions [42]. The remaining β-functions β i 16 , β i 25 , and β i 29 are new in the theory with axions and have been worked out here for the first time. We thus have a finite NNa-vertex with the Feynman rule: For the following estimation of the coupling strengths of the axion-proton and the axionneutron vertex, we assume the rest frame of the incoming nucleon, i.e. v = (1, 0, 0, 0) T , so where s µ ∆q = p|qγ µ γ 5 q|p , with s µ the spin of the proton. The proton and neutron matrix elements are related by isospin symmetry, i.e. p|ūγ µ γ 5 u|p = n|dγ µ γ 5 d|n , p|dγ µ γ 5 d|p = n|ūγ µ γ 5 u|n , and p|qγ µ γ 5 q|p = n|qγ µ γ 5 q|n for q = s, c, b, t. In particular we find g p a = − (1 + 2z)∆u + (2 + z)∆d + 3w∆s 1 + z + w for the case of axion-proton interaction, and g n a = − (2 + z)∆u + (1 + 2z)∆d + 3w∆s 1 + z + w for the case of axion-neutron coupling.
We extract the respective quantities from the recent FLAG review [43], which are here given in the MS scheme at the scale µ = 2 GeV (we utilize the nucleon matrix elements calculated on the lattice with N f = 2 + 1 excluding isospin breaking effects). The LEC d 18 is taken from Ref. [44] (fixed by the Goldberger-Treiman discrepancy), whereasd 16 has been adopted from Ref. [45]: where the nucleon mass in the chiral limit m =m N has been estimated via the third order relation (which is the accuracy to which we are working) with c 1 = −1.07(2) GeV −1 at O(p 3 ) from Ref. [44]. Note that the pion mass is taken to be the leading order pion mass M π = B(m u + m d ). Thed i 16 's are of course hitherto undetermined since they are new LECs in the theory with axions.
Inserting these values, one finds for the interaction of axions with neutrons. Note that g p a and g n a are nothing but the pure leading order coupling strengths, which were reported already in Ref. [ Since the values of the remaining LECs appearing in this expression are unknown, we estimate the strength of the contribution of g N 2 LO a by assuming that the undetermined LECs are of O(GeV −2 ) following the conventional naturalness arguments, see e.g. [45,46]. We expect the values of these LECs to be comparable to the ones known fromd 16 Note that we are still working at the matching scale µ = 2 GeV (in contrast to [28]).
Collecting all contributions, one gets for the KSVZ model with X q = 0, The coupling of axions to nucleons hence is always non-zero in both models, even though g KSVZ ann = 0 is possibile within the error range. In the DFSZ model, the strength of the coupling to protons g app can range from −0.581(58) at sin 2 β = 0 to −0.143(69) at sin 2 β = 1.
The coupling to neutrons may take on values from +0.283(55) at sin 2 β = 0 to −0.132 (67) at sin 2 β = 1, which means that in the DSFZ model g ann might still vanish depending on the value of β.

IV. SUMMARY
In this work, we have calculated the axion-nucleon couplings at the next-to-next-toleading order in two-flavor non-relativistic baryon chiral perturbation theory. Including all phenomenological knowledge from a variety of sources, we find the N 2 LO corrections of a few percent only. These couplings are therefore pinned down to a high precision.
Although we have reached a higher accuracy in the framework of chiral perturbation theory, the errors are still relatively large for two reasons. Firstly, there are still sizeable uncertainties stemming from the LO nucleon matrix elements calculated on the lattice, and secondly there are considerable uncertainties from the undetermined LECs. In fact, a possible future detection of the axion could be used to determine these LECs by applying the method of Bayesian inference. This is, however, not expected to be the case in the near future, because one would need more precise determinations of all other involved quantities such as the nucleon matrix elements ∆q or the quark masses (which might be achieved in a few years in lattice QCD), and at the same time very precise measurements of the axion-nucleon coupling would be necessary. The formula in Eq. (49) is hence primarily of relevance particularly for any future study on the axion-nucleon interaction, since the numerical values for the respective couplings can always be brought up to date by using the presented formulas and inserting the most recent estimations for the involved quantities.
In the future, it would be interesting to work out explicitly the chiral corrections to the axion-photon interaction and the influence of the strange quark.