Non-derivative Axionic Couplings to Nucleons at large and small N

Among the possible CP-odd couplings of the axion to ordinary matter, the most relevant ones for phenomenology are the Yukawa couplings to nucleons. We analyze such non-derivative couplings within three different approaches: standard effective field theory, the Skyrme model and holographic QCD. In all the cases, the couplings can be related to the CP-odd non-derivative couplings to nucleons of the low-lying mesons and the $\eta'$. Using the effective field theory approach we discuss how to derive the expressions for the CP-odd interaction terms as functions of the parameters of the effective Lagrangian at generic number of colors $N_c$ and flavors $N_f$. Then, we compute the CP-odd couplings to nucleons of the axion, the $\eta'$ and the pseudo-Goldstone mesons in both the Skyrme and the holographic QCD model with $N_f=2,3$. We present model-independent expressions for the coefficients of the non-derivative axion-nucleon couplings. This allows us to provide quantitative estimates of these couplings.

D. On-shell nucleon vertices in the limit of large mass 53 [25]. The same D-branes, wrapped around some appropriate Euclidean cycle, provide instanton effects in string theory [26]- [29]. The nature of these effects depends on the amount of supersymmetry. In the case of maximal supersymmetry they do not generate a potential, but affect higher derivative terms like the R 4 corrections, [30,31]. In all cases, the end result is that the original shift symmetry is broken to a discrete subgroup. 2 In the above descriptions (and in other QFT constructions), ALPs are fundamental fields in the theory. However this is not the only possibility. ALPs can also be composite or emergent. In its most common realizations [34]- [40], [15], the composite axion is identified with a neutral pseudo-Goldstone boson of an extra, strongly coupled, confining gauge theory, whose dynamical scale Λ a determines the axion coupling f a . In these models, the axion can be rendered "invisible" only if Λ a Λ QCD , [15]. As a common feature, the Peccei-Quinn symmetry in these models acts on some extra fermion fields. Moreover, as a difference with fundamental axions, in simple composite axion models there is no direct coupling between the axion and the lepton sector.
In a more general setup, recently considered in [41] and motivated in [42], ALPs emerge as composites of a hidden sector with a large number of colors (N c 1), which couples to the Standard Model (SM) via interactions that are irrelevant in the IR. This fact ensures that one can have an invisible emergent axion even if the strong coupling scale Λ h of the hidden gauge theory is much smaller than Λ QCD . These emergent/composite axions do not rely necessarily on fermionic Peccei-Quinn like symmetries, but rather on the approximate Peccei-Quinn symmetry associated to instanton densities in gauge theories.
The holographic correspondence [43] provides a paradigm for such a type of composite/emergent axions. In the holographic context, in fact, composite operators in a QFT are mapped into fundamental fields in the dual string/gravity theory. Hence, standard string theory axions can be related to composite operators in QFT. This is the case for the QCD topological θ term [44]- [51].
The holographic QCD axion model recently proposed in [52] provides a further different way of realizing a composite axion. In the model there is no hidden gauge theory sector, hence no extra dynamical scale. The Peccei-Quinn symmetry acts on an extra massless quark flavor with a (non-local) Nambu-Jona-Lasinio (NJL) quartic interaction. The latter has a suitable higher dimensional UV completion in the holographic model and it induces the condensation of the extra quark at a scale M a , of the order of f a , set by the NJL coupling. If M a is taken to be much higher than the QCD chiral symmetry breaking scale, the extra quarks and the related massive hadrons essentially decouple from the Standard Model. The only exception is the pseudo-Goldstone boson of the extra chiral symmetry breaking. This is identified with the composite QCD axion. Its interactions with the low energy QCD modes are precisely described by the axion-dressed chiral Lagrangian [14].
The couplings of the QCD axion and other ALPs to ordinary matter are crucial in axion searches. They typically receive contributions both from the UV realization of the ALP theory and from the QCD chiral Lagrangian and the structure of the axial anomaly [10,9,11,13]. These are combined with some input from QCD matrix elements that can be extracted from experiments or, more recently, from lattice data [53]. State of the art derivative axion couplings to nucleons are obtained in this way.
There exist also non-derivative couplings to nucleons that so far have not received much attention (see however [54]). There are two reasons for this. The first is that such couplings, when they arise from the CP-preserving effective Lagrangian for the nucleons, are effectively derivative couplings. This is the case when the nucleons are near the mass-shell, as it happens at low enough energies. If the couplings arise in the CP-violating sector of the effective Lagrangian, then they are suppressed by the small CP violation due to a residual θ angle and/or to SM effects. 3 The second reason is that such non-derivative couplings to nucleons depend on non-trivial and non-perturbative QCD physics, and chiral symmetry considerations (and the chiral Lagrangian) are not helpful enough to pin them down.
In this paper we will undertake an analysis of such couplings. There are two main motivations behind our analysis. The first is that ALPs are generically speaking light particles and the range of allowed masses towards the lower end has expanded substantially in recent years, especially as candidates for dark matter, [17,18,20,59]. In such cases, non-derivative couplings, even if small, can become important in coherent collections of ALPs. The second is that we now have a new theoretical tool to address non-perturbative nucleon physics and this is the holographic correspondence.
In the holographic approach, QCD is modeled by a large N c strongly coupled gauge theory, which can be studied by means of a dual weakly coupled string model on a classical gravity background. The top-down holographic model which better reproduces the main features of low energy QCD is due to Witten, Sakai and Sugimoto (WSS) [60,61]. Nucleons in this theory have been first studied in [62,63,64]. They are realized in a way which resembles the traditional Skyrme picture [65], where nucleons arise, at large N c , as solitons of the chiral Lagrangian [66]. In the WSS model they are instanton solutions of a five dimensional U (N f ) gauge theory (N f being the number of flavors) on a certain curved background. After reduction to four dimensions this theory encodes not only the chiral Lagrangian (including the Skyrme term) for the QCD-like pseudo-Goldstone bosons and the η , but also the effective Lagrangian for the massive axial (vector) mesons. Moreover, couplings and masses for the mesons are determined by the few parameters of the model. The holographic QCD axion in [52] is embedded in the WSS model.
There are moreover bottom up models for holographic meson and baryon physics that have been analyzed in the past. The first and crudest of them named AdS/QCD, [67,68], reproduced the physics of the chiral Lagrangian to leading order and has given rather good predictions of parameters for the chiral Lagrangian beyond the two-derivative level. There are more sophisticated bottom-up models for QCD, like V-QCD, [69,47] that improves importantly on AdS/QCD.
In this paper we embark in an effort to determine the non-derivative couplings of the axion to nucleons under three different perspectives: standard effective field theory, Skyrme model and holographic WSS model. 4 In our analysis of the Skyrme and holographic model, the large N c limit is taken by construction. In all the three approaches, the axion couplings can be related to those of the η and the pions to the nucleons, after considering the mixing between the axion and these mesons. One loop corrections in the effective theory of nucleons and mesons can contribute too. The mixing arises diagonalizing the mass matrix in the low energy effective action. It is worth emphasizing that meson-nucleon couplings have been largely studied in chiral perturbation theory. For instance, the CP-odd pion-nucleon coupling has been computed in e.g. [70,71,72,73] for N f = 2, 3 at various orders in the perturbative expansion. We will show how our results in the Skyrme and holographic model compare with those findings.
A word of caution is in order. The Skyrme and holographic models, while reproducing many features of low-energy QCD, are certainly not suitable for precision physics at percent level as there is no formal parameter that controls their accuracy. Nevertheless, their use alongside the standard Chiral Lagrangian approach is justified by the following considerations. To begin with, the CP-odd sector of the axion physics is difficult to be analyzed directly within any first-principle model (e.g. one has a sign problem on the lattice). Moreover, while in the original Skyrme model the nuclear binding energies are notoriously overestimated, its extension including vector mesons improves significantly this issue [74]. Finally, bottom-up holographic models for meson physics have been surprisingly successful in reproducing experimental data, and organizing (subleading) parameters that are theoretically unconstrained in the chiral Lagrangian. The holographic WSS model is a top-down model and has a proven discrepancy with experimental results of less than 30% (in average) when dealing with static properties of nucleons, including couplings similar to the ones computed in this paper [75]. It is also worth mentioning that if, in principle, Chiral Lagrangian calculations (combined with large N c arguments) can be systematically improved, in practice they are limited by unknown low-energy constants.
Large N c models like Skyrme and WSS can provide model-dependent estimates of such unknowns. 5 Our main results are those for the CP-odd non-derivative couplings of the pseudoscalar mesons and the axion to the nucleons. They can be summarized as follows.
• Using chiral perturbation theory, we derive the expressions for the CP-odd interaction terms, at any number of colors N c and flavors N f , between η or the axion and nucleons (transforming under generic spin-half representations of the gauge group) in terms of the parameters of the effective Lagrangian. The main results are the Lagrangian terms in equations (3.22)-(3.23) (for contributions arising through a "residual" θ-angle, possible even in the presence of the axion 6 ).
• We estimate, in turn, the one-loop contributions to the non-derivative couplings in effective field theory using previous results, [76] and extending them to general N c , N f .
• In chiral perturbation theory, we write the CP-odd non-derivative couplings of the axion to the nucleons both in terms of the CP-odd couplings of the pions and η , equation (4.14), and, alternatively, in terms of the pion-nucleon sigma term and the strong contribution to the neutron-proton mass splitting, equation (4.21), extending the formula in [54].
• We compute the CP-odd couplings (due to the effective θ angle and other sources of CP-violation in the QCD sector) of the axion, the η meson and pions to nucleons in the Skyrme model. The main results for N f = 2 are given in equations (5.19) and (5.20), and for N f = 3 in equations (E.8), (E.11) and (E.12). To the best of our knowledge, these couplings have not been computed in the literature so far.
• We compute in turn novel results for the same couplings in the WSS model.
The results for N f = 2 are given in equations (6.23), (6.24), and (6.25). The generalization to N f = 3 is analogous to the computation in the Skyrme model.
In section 7 we will compare the numerical estimates of the axion couplings to nucleons obtained in model-independent expressions from chiral perturbation theory plus lattice and phenomenological inputs, with the ones obtained in the Skyrme and holographic WSS models. The first ones are affected by some uncertainties on the value of the sigma term, while the latter depend on the specific models. Although there is a certain discrepancy in the results, a reasonable range of the value of the couplings appears to bec having in mind an effective Lagrangian term of the form δL =c N aN N , where N is the nucleon doublet and a is the axion. Here we do not distinguish the proton and neutron couplings, whose difference is smaller than the proposed window of values. Also, we have taken into account the one loop corrections in the effective theory, discussed in section 3.3. They are somewhat smaller compared to the numbers in The paper is organized as follows. In section 2 we review the construction of the QCD chiral Lagrangian, containing the η meson and the pions, coupled with the axion. In particular we write down the field redefinitions suitable for diagonalizing the related mass matrix. In section 3 we discuss the effective theory describing the coupling of mesons and axions to nucleons and deduce the general structure of the non-derivative axion-nucleon couplings. One loop corrections are discussed here too. These CP-odd non-derivative axion couplings are expressed in terms of the pion and η couplings to nucleons, and in turn of the pion-nucleon sigma term and protonneutron mass difference in section 4. In section 5 we compute the couplings in the Skyrme model with N f = 2 QCD flavors. In section 6, after a short review of the description of nucleons in the holographic setup, we perform analogous computations in the WSS model with N f = 2 QCD flavors. A summary and an analysis of the results obtained in sections 4, 5 and 6 is given in section 7. Review material, further technical details and the results for N f = 3 QCD flavors in the Skyrme and WSS models are provided in the appendices.

The effective Lagrangian couplings of the axion to the η
The general structure of axion models and their UV completions are summarized in appendix A. Here, we examine the effective Lagrangian (containing the pions and the η meson) for QCD coupled to the axion. Our aim is to write the effective theory for generic axion couplings rather than choosing a specific "frame", which allows us to check that the final results, in term of physical states, are frame independent. This will be useful in particular when considering the couplings of the η and the axion to the nucleons as we will do in Sec. 3.
The QCD Lagrangian with generic renormalizable axion couplings 7 is given by Here q L/R = (1±γ 5 )q/2 are the quark fields, we suppressed flavor indices, used mostly plus convention for the Minkowski metric, and also included the axion dependence in the mass terms M (a) as well as the axion kinetic term. The θ-angle can be absorbed in the constant term of a so we set it to zero -effects due to finite θ will be considered later on. We also included the current term ∝ ∂ µ aJ µ q5 here, but as we argue in Appendix B, it does not affect those mixing terms between the axion and mesons which we are interested in. Therefore we will set the current to zero in the following.
We shall then derive the low energy effective action for pions and the η mesons which contains the axion couplings determined by (2.1). Axion couplings in effective theories for QCD have been considered in earlier works [77,11], see for example the fairly recent detailed analysis [13]. In our analysis, it is essential to also include the effect of the η meson for two main reasons. First, we are interested in the couplings of the axion both at small and large N c , and the axial anomaly of QCD is suppressed at large N c so that the η meson is as light as the pions. Second, as we shall see, the mixing between the η and the axion is particularly important for the non-derivative couplings of the axion to the nucleons which we will analyze below. Earlier work considering the effect of the η in the axion-meson mixing include [12,78]. Since we are mostly interested in the mixing of the axion with the mesons, we will concentrate on the terms giving rise to quadratic terms in the action.
The first interaction term in (2.1) couples, in the effective theory, the axion to the pseudoscalar glueball field [77,79]. The glueball can be integrated out, giving rise to the "standard" effective Lagrangian with the θ-angle replaced by the axion field. The meson fields are included in where Π = π a T a . The matrix U transforms as where L, R are U (N f ) matrices. We normalize the group generators as Notice that η is not a mass eigenstate but will mix with the other fields, as we will discuss below.
With these conventions, the effective Lagrangian that is the low-energy theory of (2.1) (with the axion derivative/current term omitted) can be written as, 8 where M π is the pion mass matrix, and χ YM is the Yang-Mills topological susceptibility.
We allowed the η to have a different decay constant than the pion [80] by including the last term in (2.5). Requiring the standard kinetic term normalization fixes The meson mass matrix is proportional to the quark mass matrix, M π (a) = B 0 M (a) up to higher order chiral corrections. We assume that it is possible to write the mass matrix in the form M (a) = e iaQ/fa M 0 e iaQ/fa , (2.7) (possibly after vectorial axion-valued transformations) where 9 M 0 is Hermitean and independent of a. It is then straightforward to solve explicitly for the mixing between the mesons and the axion. Effects due to an anti-Hermitean component in the mass matrix are discussed in Appendix B. We define the topological susceptibility of full QCD, χ, and the "gauge invariant" coupling of the axion,c s , as In order to diagonalize the Lagrangian, we first carry out the following field redefi-nitions (see also [12]): 10 which remove the mixing of the pions and η with the axion up to terms suppressed by higher powers of 1/f a . More precisely, the mixing terms in Eqs. (2.9) and (2.10) are determined by the non-diagonal elements of the mass matrix in (2.5). Notice that the terms involvingQ are (up to nonlinear terms in the pions) equivalent to a chiral rotation which removes the a dependence of M (a). Then the other terms are sourced by the gauge invariant couplingc s . After applying (2.9) and (2.10), the transformation (2.11) removes the non-diagonal kinetic terms. Inserting these relations in the Lagrangian (2.5), the quadratic terms (with terms suppressed by 1/f a dropped) are included in 13) and the axion mass squared is given by (2.14) The fieldsη andπ b are however not mass eigenstates (for a generic quark mass matrix) because they are subject to the "usual" mixing of neutral states imposed by the quark masses. The final physical states η p and π p are found by diagonalizing the mass matrix arising from (2.12) which in general needs to be done numerically.
CP-odd couplings and a finite θ-angle (which may be induced by linear couplings of the axion to other sectors not considered here) introduce small VEVs for the mesons. Such VEVs will be important (apart from the mixing between the axion and the mesons) when considering their couplings to the nucleons in Sec. 3. We discuss these VEVs in Appendix B.1. We also include results for the CP-violating vertices of three pseudoscalar mesons in Appendix B.2.
Finally, the scaling of the QCD parameters at large N c and/or N f for our conventions is the following. The decay constants scale as while the term ∝ c η , which induces the splitting of the pion and η decay constants, is suppressed by a power of 1/N c N f since it is the coefficient of a double trace term. We normalized this coefficient such that it scales as In this case, the difference of the decay constants also obeys We also note that

The effective couplings of η and the axion to nucleons
In this section we proceed further and discuss the effective field theory that couples the mesons and the axion to the nucleons. In particular we analyze the CP-violating non-derivative couplings of the η and the axion to the nucleons. 11

Generic N f and N c
We shall discuss here the effective theory for the nucleons and their couplings to the pions, η and the axion for generic values of N c and N f . Specific results for N f = 2 will be given in Sec. 3.2. We start by writing down the couplings in the absence of the axion and the θ-angle, and then discuss how they are introduced in a generic frame. Working in a generic frame will help us to check that the results behave correctly under chiral transformations. We shall work at any odd value of N c , and for simplicity only consider spin 1/2 baryons but in arbitrary representations of the flavor group. In order to describe the couplings of the baryons to the mesons and the axion, the standard method, [81], is to first define a new matrix through u 2 = U . Because the square root of a matrix is ambiguous, this is not enough to make the definition precise; we can take This matrix transforms as under the chiral transformations (see, e.g., [82]). Here K is a unitary matrix depending on the pion fields and defined 12 through the equation in (3.2). The baryons belong to the N c -index representations 13 of (the vectorial) SU (N f ) and transform according to [87] 3) The baryon bilinears may therefore couple to meson operators X transforming as X → KXK † . Such operators include where M 0 is the quark mass matrix and u was defined in (3.1). The derivative couplings of the pions and the η arise from the termŝ where "k ∼ contractions" stands for all possible (invariant) contractions of the flavor indices of the two nucleons and the meson matrix. We have also scaled out the leading expected behavior of the coupling constants on N c and N f so that the expectation iŝ Notice however that we may still have subleading dependence involving terms ∼ 1/N c and ∼ 1/N f . We will follow similar conventions for other couplings below, unless stated otherwise. In the second term of (3.5) we include a sum over various contractions (with couplings g a ,. . . ) of the flavor indices of the baryons and u µ before taking 12 For a precise definition and discussion, see, e.g., [81,82]. 13 One should keep in mind that in the limit of large N c there would be more relations than what we impose here: the nucleon states are arranged according to the spin-flavor symmetry SU (2N f ), with the ground state being the fully symmetric N c -index representation in this limit. These representations break into towers of states classified according to their spin and flavor representations separately [83,84,85,86]. the trace, which are consistent with chiral symmetry. Not all of the possible terms are independent, but their relations are nontrivial for baryons transforming under generic representations of the chiral symmetry. The nontrivial piece amounts to combining the baryon, the part of the meson term transforming under the adjoint representation, and the anti-baryon in the conjugate representation into a singlet. We will explain the structure in more detail in Appendix C.3, considering N f = 2 and N f = 3 as examples.
The precise knowledge of the structure of these terms is not needed for most of our results. In the limit of large N c it has been solved in [86]. Moreover, for the coupling of the flavor singlet η meson, these terms reduce to the standard trace, i.e., they have the same form as the first term, but are leading by a factor of N c . The suppression of the double trace term is consistent with the expected enhancement of the chiral symmetry 14 to U (N f ) L × U (N f ) R in the limit of large N c . In other words, (3.5) may be replaced by (3.7) where we projected to the adjoint representation in the second term by rearranging the singlet contribution. Therefore  In this article, however, we are interested in the non-derivative couplings to the nucleons, which we will discuss next.

Non-derivative couplings
We go on to discuss the non-derivative couplings of the nucleons to the η , the pions, and the axion. Such couplings are absent when CP is conserved, but they can arise through a nonzero θ-angle which in turn can result from additional UV terms in the action, which drive the expectation value of the axion away from zero. Other CP-violating contributions than those appearing effectively through the θ-angle are also possible, we analyze them in Appendix C. We first write down the relevant terms in the effective Lagrangian in the absence of the axion and the θ-angle, and then discuss how the axion and θ are included. A class of non-derivative couplings arises from the mass terms involving χ ± (that were defined in (3.4)): (3.11) Several comments are in order. The terms in (3.10) give rise to corrections to the nucleon masses, and they are ∼ O(N c ), in agreement with the scaling of the nucleon mass. Notice that we separated again a piece from the flavor singlet contributions in order to write the more complicated non-singlet terms in terms of the combinations χ ± − Tr χ ± /N f which vanishes (up to higher order interaction terms in the pions) for flavor independent quark masses.
The terms containing an extra γ 5 are suppressed by a power of 1/N c . Considering the singlet term as an example, this can be seen from the identity where the first equality is exact when the nucleons are on-shell and we used the fact that the baryon mass m B ∼ N c (see also Appendix D). From (3.12) we also observe that the terms involving γ 5 are, moreover, derivative terms in disguise and they are suppressed in the non-relativistic limit. This is not entirely true for the single trace "contraction" terms in (3.11), because of the flavor structure: the single trace terms of (3.10) give rise to flavor dependence of the nucleon mass which leads to a non-derivative contribution from (3.11). These contributions are of the second order in the differences χ ± − Tr χ ± /N f . The single trace term in (3.11) will not, however, contribute to the linear axion couplings because the axion dependence also cancels at linear order due to the subtraction of the singlet piece after substituting (2.9) and (2.10) as we shall show below. This holds true also for the mixing effects between the pions (arising from expanding the χ − ) and the axion. Therefore we shall not discuss the terms in (3.11) further.
Another interesting class of terms are the flavor singlet terms which couple the η ∝ log det U to the nucleon singlets. Possible terms are (3.14) The scaling with N c was determined as follows: every factor of i det log U brings a power of 1/N c with respect to the nucleon mass term and the terms in (3.10) and in (3.11) [88]. Recall that the term d 1 is in addition suppressed by 1/N c due to the γ 5 factor. In order to extract the non-derivative couplings from these Lagrangian terms, we consider the axion couplings induced by the axion terms in the QCD Lagrangian (2.1) and turn on a finite θ-angle. We replace 15 and where M 0 is assumed to be Hermitean. In terms of the axion, the θ-angle is therefore introduced effectively by replacing Notice that also the meson mixing in (2.9) and in (2.10) will be affected: a finite θ-angle leads to a small displacement of vacuum expectation values (VEVs) of the pions and the η from zero. 16 Other sources of CP-violation than the θ-angle can also lead to such VEVs. We analyze them in Appendix B.1 and as it turns out, they do not lead to terms with essentially different structures than those terms which we derive in this section. In Appendix B we also analyze the effects of non-hermiticity of the quark mass matrix and show that they can be absorbed in the θ-angle.
In order to extract the contributions involving θ, we therefore need to study the nucleon couplings which are quadratic in the pion fields, the axion, and η (at vanishing θ). First we notice that where χ ± (a) = uM (a) † u ± u † M (a)u † , the anticommutator is normalized as {A, B} = AB + BA, and we already inserted the relations (2.9)-(2.11) but so far kept θ = 0. We remark that in (3.20) the axion dependence is proportional to the identity matrix in flavor space. 17 This verifies the fact that it is indeed absent in the second term of (3.11) (which we have already omitted). Notice that we discuss here only the trace term ∼ Tr [χ + (a)], see Appendix C.2 for the analysis of the second term in (3.10).
The terms which give rise to quadratic couplings of η and the axion to the nucleons are those in (3.10), the second term in (3.13), and the term in (3.14). Inserting (3.18), (3.20) and (3.21) in these expressions, and adding θ through the replacement (3.17), we find the following couplings: where we neglected second-order terms in the quark masses. Since the axion coupling is proportional to c m1 , which also controls the nucleon mass corrections, there is a relation between the two (in the absence of other sources of CP-violation than the θangle). We will discuss this in more detail in Sec. 4. Notice also that theη couplings on the second line (3.23) are suppressed by N f /N c with respect to the couplings on the first line. Therefore the couplings on the first line should be compared to the Skyrme and WSS model analyses carried out in sections 5 and 6, which employ the 't Hooft limit. Contributions from the single trace term in (3.10) are analyzed in Appendix C.2: these include the CP-odd nucleon-pion couplings (C.12), and a correction to the nucleon-axion coupling (C.13), which is nonzero if quark masses depend on flavor. The CP-odd nucleon-η coupling is independent of these terms. Recall also that the mixing of η and the pions due to flavor dependent quark masses is not included in these results.

Results for N f = 2
We now discuss the structure of the couplings at finite, fixed N f , taking N f = 2 (and any odd N c ) as an example. This will be useful in order to compare to the results of sections 4, 5 and 6. In addition, results for N f = 3 are given in Appendix C. First, we analyze the derivative couplings in (3.5). The three pions and the η 18 form a four-dimensional representation of U (N f = 2) which breaks into the iso-singlet η and the iso-triplet of pions at finite N c . The ground state, spin-1/2 fermions, are in the fundamental isospin-1/2 representation: In this case the derivative couplings of (3.5) and (3.7) read a +ĝ a /N c , and τ a are the Pauli matrices. 19 There is only one independent term in the sums over contractions in (3.5) and (3.7). The flavor structures of the mass terms (3.10)-(3.11) are analogous to the derivative terms. We note that in the limit of large N c , the derivative couplings of the η and the pions become identical since f η /f π → 1 in this limit, signaling the enhancement of the chiral symmetry from SU (2) to U (2).
When N f = 2, it is natural to take the light quark masses to be equal, m u = m d ≡ m ud , so that isospin is unbroken and η does not mix with the pions. In this case the couplings of the previous subsection apply directly for the physical mass eigenstates a p , η p =η , and π b p =π b . The couplings of the axion and the mesons to the nucleons therefore read The terms on the first line arise from (3.22) and (3.23). We also included the pion coupling on the second line, which depends on the single trace coupling c (1) m1 in (3.10); see eq. (C.12) in Appendix C.2. Recall that in this case the pion mass is given by m 2 π ≈ B 0 m ud in the chiral limit. We note that the same coupling also controls the QCD contributions to the mass splitting between the proton and the neutron, so that there is a relation between the the physical one, but it can be related to the physical η P and η P fields. In [89] for instance, taking into account the pseudoscalar mixing angle, it is taken as η ≈ 4 5 η P + 3 5 η P , so its mass is around 549 MeV. 19 In this action we neglected interactions with other nucleon representations like the ∆.
mass splitting and the CP-odd pion-nucleon-nucleon coupling (see, e.g., [72]). We will discuss this relation in more detail below. Similarly to the derivative couplings above, c m1 equals c (1) m1 at large N c up to corrections suppressed by 1/N c . Since also f η ≈ f π in this limit, the couplings of the pions and that of the η in (3.28) to the nucleons become equal, signaling the enhanced chiral symmetry.

The one-loop corrections
The meson-nucleon couplings of the effective Lagrangian described above give rise to further contributions to the couplings of interest via quantum effects. These quantum effects are determined by the effective meson and nucleon fields. In the 't Hooft limit, they are suppressed by powers of N c , while in the Veneziano limit, they may be unsuppressed.
We will discuss such effects below. In the case of N c = 3, N f = 2, they were calculated in [76] but we will extend here this calculation to arbitrary N c and N f .
The starting Lagrangian for the pion-nucleon interactions was given already in the previous subsections, but we will abstract here the couplings that are of interest. The η and η couplings to pions from the chiral Lagrangian are where the couplings behave as in the chiral and large N c , N f limits (see Appendix B.2). Similarly the coupling to the axion can be written as where f aππ scales as ∼ θm ud /N 2 f . The couplings to the nucleons are 20 We also need the kinetic terms of the spin 1/2 nucleons and the spin-3/2 ∆, which in the conventional N f = 2 case are  N above stands for the spin-1/2 nucleon iso-doublet while ∆ µ for the ∆ iso-quartet. The covariant derivatives D µ and u µ were defined in (3.4) and the matrix u in (3.1). T a are the Clebsch-Gordan matrices that couple the doublet, the quartet and the triplet of SU (2). From (3.26) we have where the last equality above with the coupling of the nucleons to the axial current being valid in the chiral limit. At arbitrary N c and N f , the field N represents the lowest-lying flavor representation of spin 1/2, that is the Nc+1 2 , Nc−1 2 , 0, · · · , 0 in Dynkin notation. The spin 3/2 representation is the Nc+3 2 , Nc−3 2 , 0, · · · , 0 flavor representation and so on. In the general case, there may be also higher-spin baryons with spin larger than 3/2 and up to N c /2. The parameters m N , The flavor structure of the pion-nucleon interaction involves T a ij where a is an adjoint index of U (N f ) and i, j are indices in the flavor representations R, R of the nucleons. At general N f , these can be any representation in the tensor product (⊗ ) Nc . We will be assuming here exact flavor symmetry so that all quark masses are the same.
We discuss first the one-loop couplings of the η and η mesons to the nucleons and consider the axion coupling below. The diagrams that contribute at one loop are shown in figure 1. Diagram (a) has ground-state nucleons as intermediate states while diagram (b) has higher spin nucleons. Diagram (c) gives a vanishing contribution due to isospin symmetry as charge conjugate pions give opposite sign contributions. Notice that the CP-violation must arise from the three-meson vertex because this vertex vanishes if CP is conserved.
For the rest, we will call R the flavor representation of the ground-state nucleons and R I all the others that give rise to the analogues of ∆. The flavor couplings that enter in L ∆N π in (3.34), F a i,α , are the Clebsch-Gordan coefficients linking the adjoint, R andR I of the flavor group. The index i is that of the ground-state nucleon while α is that of the ∆. Of course for higher spin nucleons, ∆ µ 1 µ 2 ··· is a tensorial spinor, and the space-time factors in the couplings in (3.34) will change appropriately. These are O(1) changes though and we will not worry about them at the moment.
For the computation of the loop diagrams we follow Ref. [76]. The covariant extended on-mass shell scheme is used for renormalization [90,91]. The diagram (a) gives at N f = 2, [76] where the renormalization scale was set to the nucleon mass and the external nucleons were put on-shell. In the general N f , N c case, where the coupling is between the nucleon states N i and N j we must multiply (3.36) by (T a T a ) ij = I(R)δ ij of the flavor group where I(R) is the standard quadratic Casimir of the relevant representations. The axial coupling in (3.36) scales as The ∆-mediated diagram has the contribution where the dimensionless function Z can be found in [76]. In the large N c limit it scales as At general N c and N f , the diagrams of type (b) corresponding to the coupling of external nucleon states N i and N j are multiplied by where (F a I ) i,α are the Clebsch-Gordan coefficients linking the adjoint (a index), the ground-state flavor representation (i index) and the higher spin intermediate flavor representation I (index α). There are also O(1) differences due to the different spin structure of the intermediate ∆ states. The conclusion is that all such individual contributions scale at large N c as in (3.37). 22 For N f = 2, N c = 3 from [76] we obtain The results in the extended on-mass shell scheme are suppressed by about 30% with respect to the heavy-baryon limit [76]. The one-loop couplings to the axion can also be extracted from the above results: this amounts to replacing f Hππ m H by f aππ /f a in (3.36) and in (3.38). We find that Finally, we can estimate these loop contributions as a function of the θ angle as follows. By using the expressions for the three meson couplings in chiral perturbation theory, (B.31) and (B.32), we find that where we took the limit m u + m d m s , and used m π 135 MeV, m η 958 MeV, f π 92 MeV, and θ ηη −π/9 for the mixing angle defined in (B.21). We obtain for the CP-violating one-loop couplings Similarly, from (B.34) (settingc s = 1 and taking the chiral limit with two light quarks) we find that so that the one-loop axion couplings can be estimated as

General IR relations for the couplings
Before considering large N c models in detail in the following sections, let us consider model-independent low energy relations between the couplings valid also at small N c (partially recollecting what we have learned in the previous sections). We consider the chiral Lagrangian (2.5) choosing the scheme where the axion is entirely contained in the topological term [14] where U is defined in (2.2). We concentrate on N f = 2 flavors, where the quark mass matrix reads M = diag(m u , m d ) and T a = τ a are the Pauli matrices. The GMOR relation sets As we have seen in section 2, (4.1) generates a mixing of the η with the axion, and in turn the latter with the pions. Neglecting the mixing terms between η and pions, the mass eigenstates at leading order in 1/f a are now (see equations (2.9)-(2.11)) is the full topological susceptibility of the theory. The previous expressions show that the mixing of the axion with pions only occurs if the flavors are non-degenerate in mass. The couplings of the η and the pions to nucleons are defined via the following effective Lagrangian terms 5) where N is the two component vector for the spin 1/2 nucleons N = (p, n), the first two couplings are CP-even, while the last two are CP-odd and arise in the presence of any effective non zero θ-angle or any other CP-violating source.
The axion couplings to nucleons, to leading order in 1/f a , can be obtained from the ones above, simply replacing η and π a with the respective physical fields defined in (2.9), (2.10). Notice that in the non-relativistic limit, which is justified at energies much smaller than the nucleon mass m N , the non-derivative CP-even couplings can be traded for derivative ones since, for any scalar field φ i∂ µ φN γ µ γ 5 N ≈ 2m N φN γ 5 N . (4.6) These meson-nucleon couplings were already analyzed using effective field theory in Sec. 3; the results for derivative couplings are given in equation (3.27) and for the non-derivative CP-violating couplings in equations (3.28) and (3.29).
If we define the derivative axion-nucleon couplings as we find, using (4.6), (4.5) and (4.3) To leading order in the chiral limit, where the mass of the pion is much smaller than that of the η (i.e. when f 2 , the expressions above can be rewritten as where the isoscalar and isovector axial couplings, obtained from the matrix elements of the axial current between baryon states (the axial form factors), are related to those of the η and the pion by means of the generalized Goldberger-Treiman relations The couplings (4.10), (4.11) are, for the case N f = 2, precisely the "universal" ones obtained in any axion model in the KSVZ class [6,7]. The non-derivative axion-nucleon couplings δL aN N,non−der =c N aN N , (4.13) are present whenever the Peccei-Quinn mechanism is not perfect (or if weak interaction CP-breaking effects are taken into account). As already recalled, explicit breaking of the Peccei-Quinn symmetry by higher-dimensional operators shifts the VEV of the axion. These operators are generically present in most UV completions of axion models. Even if suppressed by powers of the Planck mass, up to dimension ten their coefficient must be very small, in order to cope with the experimental bounds (see e.g. [92,93]) setting θ 10 −10 , [94,95]. The extent of fine-tuning needed for these coefficients determines the so-called "quality" of the axion. Therefore, on general grounds one can expect that the axion does not cancel completely the θ-angle. In such a case, we can study linear CP-breaking couplings of the axion and η to nucleons, with a (non-vanishing) θ-dependent coefficient.
Taking the axion-meson mixing (4.3) into account, we can writē (4.14) Therefore, to leading order in the above mentioned chiral limit Let us introduce the shorthand notation for the first and second term in each line of (4.15)c p ≡c 1 −c 2 ,c n ≡c 1 +c 2 . (4.16) In [54] a relation between the common termc 1 in (4.15) and the pion-nucleon sigma term δM N in the chiral limit has been derived. We will re-derive it in the degenerate case in the Skyrme and WSS models in sections 5, 6. We denote by M N the nucleon mass in the isospin-conserving limit, where the two light quarks have a common mass equal to mu+m d 2 . By virtue of the Feynman-Hellmann theorem, the sigma term is then The relation of the sigma term with thec 1 part of the axion coupling 23 reads [54] 24 23 Note that this implies a relation between the couplingḡ η N N and the sigma term. 24 In [54] the formula is obtained by considering the nucleon matrix element of the interaction operator θ fa mum d mu+m d (ūu +dd +ss). The latter is extracted, in the m u,d m s limit, from the expansion of the quark mass term in the Hamiltonian, once the θ-term has been rotated in the mass sector and the energy has been minimized.
Moreover, working again in the chiral limit, the pion coupling (and so the termc 2 defined above) is known in chiral perturbation theory to be related at tree level to the strong contribution to the neutron-proton mass splitting, which we denote as (M n − M p ) (L M ) (see e. g. [71] and sections 5, 6), as 25 Thus, using (4.15), (4.18) and (4.20), we can write universal low energy relations, valid for N f = 2 flavors in the chiral limit (m π m η ) at tree level in chiral perturbation theory and up to third order in the quark mass difference , generalizing the expression in [54], In summary, in these relations the CP-odd couplings of the axion to neutrons and protons are expressed in terms of the pion-nucleon sigma term and the strong contribution to the neutron-proton mass splitting. The latter can be estimated with phenomenological extractions or lattice calculations. We are going to comment on the numerical evaluation of these couplings in section 7.
Notice that what we have previously denoted byḡ aN N is the value of thec N couplings in the quark mass degenerate case, which corresponds to = 0 and (M n − M p ) L M = 0. We will keep this definition in the rest of the paper.

The Skyrme model picture
In this section we consider the couplings discussed above in the context of the Skyrme model, where the nucleons can be seen explicitly as solutions in the low-energy Lagrangian. In the Skyrme model, one treats baryons as solitons of the chiral Lagrangian, stabilized with the aid of a quartic term in U , the so-called Skyrme term [65]. The axion-dressed Lagrangian (4.1) thus contains now a further term The coupling in the Skyrme term scales with N c as e ∼ 1/f π ∼ 1/ √ N c . Moreover, we work at large N c and finite N f , so that f η = f π and double trace terms are discarded 25 In the N f = 3 case an analogous relation holds involving the mass difference of Ξ and Σ baryons [70]. in (4.1). Hence, for instance, in the N f = 2 mass-degenerate case m u = m d = m ud , the relations (4.3) reduce to Up to corrections of order O(1/f 2 a ), the masses of the physical states are where in the first line we have introduced the Witten-Veneziano mass m W V . The Skyrme model without the axion allows to compute the meson-nucleon couplings defined in eq. (4.5). From these, using the general formulae collected in the previous section, the axion-nucleon couplings can be deduced.
Notice that in the Skyrme model, both g πN N [66] and g η N N [89] have been computed (the latter one in a suitable extension of the Skyrme model including some vector mesons). To the best of our knowledge, the CP-odd couplingsḡ η N N andḡ πN N have not been computed in the Skyrme model so far. The couplingḡ πN N in the Skyrme model with a finite θ-angle has been argued to be suppressed in the large N c limit [96]. In the chiral limit (which in general does not commute with the large N c one) it has been computed in e.g. [70,71,72,73].
In the following subsection, considering the N f = 2 non-degenerate case, we will directly extract the non-derivative axion-nucleon couplings and, in turn,ḡ η N N and g πN N , showing in particular that the relations (4.18), (4.20) and (4.21) are indeed reproduced.

Nucleon mass terms and CP-odd couplings for
The starting point of our analysis is the computation of the quark mass contribution to the nucleon mass in the N f = 2 non-degenerate case. It can be extracted from the mass term of (4.1), set on-shell on the Skyrme solution, once we subtract the vacuum configuration U = 1 2 and using the GMOR relation B 0 (m u + m d ) = 2m 2 π . This gives describe the ansatz for the extended Skyrmion solution of [97,98] which we truncate to the η -pion sector. The generalized pseudoscalar-vector meson Lagrangian considered in [98] contains in fact various mixing and isospin-breaking terms involving the pion, the η and other vector mesons. These terms contribute to the Skyrmion equation of motion and to the nucleon mass terms. However they are expected to provide subleading corrections (see e.g. comments in [99]) and we will neglect their effect here. In eq. (5.7) J is the angular momentum and Θ the moment of inertia of the nucleon. Moreover 1, 2, 3) are the Pauli matrices and 1 2 is the two-dimensional identity matrix. To leading order in 1/N c and in the small quark mass limit, the function F is the solution of the equation [66] where F = F (y) with y = 2ef π r, and the conditions for having baryonic number equal to one are F (0) = π, F (∞) = 0. The equation can be solved numerically. The solution asymptotes as y −2 at large y.
The static Skyrmion solution is described by the matrix U s , and the quantization is achieved in a non-relativistic limit passing to a slowly rotating time-dependent solution a(t)U s a(t) −1 where the matrix a accounts for the time-dependent SU (2) (for N f = 2) moduli. Some detail on the instanton moduli quantization, valid mutatis mutandis for both the Skyrme and the holographic WSS model, can be found in appendix F.
In the extended Skyrme model of [97,98], where the Skyrmion emerges as a soliton solution of an effective Lagrangian including ρ and ω vector mesons, the quantization includes a time-dependent solution for the η soliton component [97], whose form is given in eq. (5.7) above. The function η(r) is extracted from the η equations of motion. The latter include (at least) two parameters, related to the vector-vector-pseudoscalar coupling g V V φ and the vector-pseudoscalar 3 coupling h. These have been extracted in the model at hand in [89].
The η soliton component is proportional to J and thus to the soliton angular velocity χ ∼Tr[ȧ τ a −1 ] which clearly vanishes in the static limit. Time derivatives of the moduli are usually 1/N c suppressed in the model and therefore the contribution of this term will be generically suppressed too. However, this term provides the leading order contribution to the CP-even g η N N coupling. The explicit result can be found in section E of [89] for the N f = 2 case. By standard parameter fitting the authors of [89] obtain g η N N = 1.61.
Since ϕ ∼ 1/N c , we can expand the Lagrangian term (5.6) around ϕ ≈ 0. To first order in ϕ we get sin F x a r Tr τ 3 aτ a a −1 ϕ .
The corresponding Hamiltonian term is given by where I 3 is the third component of the isospin operator and we have used the relation From (5.10), recalling that I 3 = 1/2 (resp. −1/2) on the proton (resp. neutron) state, we get the pion-nucleon sigma term and the quark mass contribution to the neutron proton mass splitting [100,98] δM (5.12) Notice that the sigma term can be given explicitly as where the constant α is determined by the numerical solution for F (y) through Moreover, the neutron-proton mass difference can be obtained by numerically solving the equation for η(r) [98]. Notice that in the more general case of N f degenerate quark masses, one would get the same expression as in (5.13), as a result of the standard embedding of the SU (2) Skyrmion solution into SU (N f ).
Taking the previous results into account, there is a quick way to get the CP-odd axion-nucleon couplings.
Following the procedure outlined in section 2, we can remove the axion-meson mixing terms in the vacuum solution in such a way that the net effect is just a redefinition of the quark mass matrix. In the "chiral limit" m π m W V (with m W V defined in (5.4)), treating the axion as an external field, a convenient choice is 26 where 2m = m u + m d . Hence, in the Skyrme plus axion case, after shifting the axion by f a θ to account for a residual θ angle, the nucleon mass term arises from the Lagrangian term where ϕ and U s are given in (5.7). Working to first order in ϕ, expanding for small a, θ and keeping only the terms in θa, we get Treating the axion as an external field, the corresponding Hamiltonian term reads which can be expressed in terms of the Skyrme model parameters using (5.12) and (5.13). This expression gives the Lagrangian CP-odd axion couplingsc N precisely as in eq. (4.21). In terms of the model parameters For the more general case of N f degenerate quark flavors, the common factor 4 in the denominator is replaced by N 2 f . Taking into account the axion-meson mixing, and the general IR relations (4.15), the previous expressions allow us to immediately findḡ η N N andḡ πN N . In the chiral limit they read where we have used eqns. (5.12) and (5.13). As a countercheck, we have controlled that the expressions above can be also obtained directly starting from the θ-dependent Skyrme Lagrangian 27 and considering the equations of motion for the η and pion fluctuations around the Skyrmion solution.
If the Peccei-Quinn mechanism works perfectly, so that no residual θ term is left over, the results in (5.19) can be interpreted as providing the quartic axion-axionnucleon couplingsḡ a 2 N N and those in (5.20) the meson-axion-nucleon couplings: it just suffices to trade back θ for a/f a . In practice in the chiral limit this gives In appendix E we will present related results in the N f = 3 non-degenerate case.

Large N c estimates in the WSS holographic model
In this section we consider a second model which allows to estimate the non-derivative couplings to nucleons of axion, η and pions. It is the most successful top-down holographic cousin of planar QCD, know as the Witten-Sakai-Sugimoto model (WSS) [60,61]. Although it has a higher-dimensional UV completion, it gives reasonable quantitative results for many QCD observables. Similarly to the Skyrmions, baryons are solitons of this theory. Moreover, the model allows to include systematically the contribution of the (axial) vector mesons.
In the next subsection, we will briefly review the model for the reader who is unfamiliar with it. In the following subsections we will compute the above mentioned CP-odd couplings in analogy to what we have done for the Skyrme model.

The WSS model
The WSS model is based on a D 4 − D 8 brane setup in type IIA string theory. There are N c D 4 -branes wrapped on a circle of radius R 4 = 1/M KK where antiperiodic boundary conditions for the fermions are imposed. Then N f D 8 -anti-D 8 -branes are placed at antipodal points on that circle. 28 At energies much smaller than M KK , the dynamics on such branes is given by a 3 + 1 dimensional large N c SU (N c ) gauge theory with N f massless quarks. The massive sector contains scalars and fermions (with mass scaling as M KK ) in the adjoint representation of the gauge group. In the limit where a simple holographic gravitational dual description of the model can be given, this massive sector cannot be decoupled from the QCD-like one. The classical gravity regime, in fact, amounts on taking N c 1 and λ 1, where the 't Hooft coupling λ sets the ratio between the confining SU (N c ) string tension and M 2 KK . The flavor sector in the model is described by the low energy modes of the D 8 -branes. When f ∼ λ 2 (N f /N c ) 1 the latter can be treated as probes of the background sourced by the D 4 's: this corresponds to the quenched approximation for the quarks. 29 In this limit, the D 8 -anti-D 8 -branes actually join in the background, realizing geometrically chiral symmetry breaking. Their effective action reduces to a U (N f ) Yang-Mills theory with Chern-Simons terms on a curved space-time in five dimensions 28 This is not the only possibility. More general non-antipodal configurations can be also considered. However, it is only in the antipodal case that the glueball and meson mass scales coincide. 29 See [103] for an account of the flavor backreaction to first order in f .

where (in units
Here z is the radial holographic direction and we have omitted the wedge product symbol "∧" in the CS terms. We shall mostly focus on such an action in the N f = 2 case. The matrix for the Goldstone modes is given by the path ordered holonomy matrix The effective action for U precisely reduces to the chiral Lagrangian with the Skyrme term [61], with the pion decay constant f π and the coupling e given by The pseudoscalar meson η acquires a mass according to the Witten-Veneziano relation where χ Y M is the topological susceptibility of the pure SU (N c ) theory, given by This implies that in the WSS model the Witten-Veneziano mass, is much smaller than M KK . Baryons in the model are instanton solutions [63] with baryon number n B given by where B is the space spanned by (x 1,2,3 , z) and F is the U (N f ) field strength. Nucleons correspond to n B = 1 solutions.
In the N f = 2 case with massless flavors, using an SU (2) × U (1) notation for the gauge field a simple static instanton solution can be given around z = 0 where the curvature of the background can be neglected. The related solution corresponds to a charged BPST instanton where M = 1, 2, 3, z and (6.11) This solution depends on eight parameters: the instanton center of mass position X M = ( X, Z), the instanton size ρ and implicitly, three SU (2) "angles" related to the fact that the solution can be rotated by means of a global gauge transformation. Precisely as in the Skyrme case of section 5.1, they are encoded in the matrix Substituting the solution (6.10) into the action (6.1) one finds where, up to O(λ −2 ) corrections where M B0 gives the baryon mass in the λ → ∞, N c → ∞ limit. This implies that ρ and Z are not genuine moduli; in fact they are classically fixed by minimizing M B as These relations imply that the instanton size ρ ∼ 1/ √ λ is very small (but not zero) in the λ 1 regime, and that the center of the instanton is classically localized at Z = 0. The full expression for the static instanton solution valid for any value of z in the λ 1 limit can be found in [64]. In [48,49], the solution above has been modified by including a (small) mass term m q for the flavor fields (taken to be degenerate in mass) as well as a (small) non-zero topological θ term. The deformation is driven (in units M KK = 1) by the (small) parameter cm q κ θ ≡ m 2 π π θ . (6.16) In the mass-and-θ-deformed instanton solution, two new gauge field components are turned on, to first order in the above parameter. One is the Abelian component A mass z , which gives a non trivial vacuum expectation value to the η meson, while the other is the non-Abelian component A mass 0 which can be interpreted as a dipole electric potential in the bulk. It is the latter component which is responsible for the nucleon electric dipole moment to be different from zero [48,49], but it is the former which will be useful for our purposes below.
In appendix F we review the instanton moduli quantization, which will be necessary for the computation of the CP-odd couplings in section 6.3.

The derivative axion-nucleon couplings
As we have pointed out in section 5 the derivative axion coupling to nucleons, for the N f = 2 case, can be simply deduced once the CP-even couplings g η N N and g πN N are known. In the WSS model in the chiral limit these couplings have been already computed in [64]. They read clearly satisfying the generalized Goldberger-Treiman relations. Using (4.10), (4.11) the derivative axion-nucleon couplings follow. It is worth noticing that the leading contribution to the isoscalar coupling g η N N comes from the time-dependent term, proportional to the angular velocity χ, in the solution for A z (F.5).

The non-derivative axion-nucleon couplings for N f = 2
In the holographic WSS model, the nucleons correspond to leading order to the static instanton solutions A = A cl (6.10) for the five-dimensional gauge field on the flavor branes. Moreover, the quantization of the solution proceeds in analogy with the Skyrme case and, as we review in appendix F, it generates (among others) a nontrivial time-dependent profile for the Abelian componentÂ z which is related to the η field. Furthermore, the quark mass term in the model, in the small quark mass limit, turns out to have precisely the same form as in the Skyrme case [104,105,48,49]. Finally, and crucially, the WSS model allows for the explicit introduction of an axion [52], giving the standard modification of the low energy chiral Lagrangian as written in (4.1). The holographic set-up provides a UV origin of the axion model, but f a is an essentially free parameter.
All in all, focusing on the pseudoscalar sector, one can extract the quark mass contribution to the baryon Hamiltonian and the CP-odd couplings precisely as in the Skyrme model in section 5.
Let us consider the "chiral limit" m π m W V with N f = 2 non-degenerate flavors. The starting point is thus the same Lagrangian mass term as in (5.6) where now and, from the time-dependent solution (F.5), where P Z is the momentum conjugate to Z and χ is the instanton angular velocity defined in eq. (F.7). We shall not need the explicit expression for F (r) since, when computing the VEV of the above expression on the nucleon ground state, this will not contribute. We thus set P Z = 0 in the following. Notice that ϕ scales like 1/N c , just as in the Skyrme case. The above ingredients are sufficient to extract the pion-nucleon sigma term and the neutron-proton mass splitting following precisely the same steps as in the Skyrme model case. To leading order in the holographic limit they are given by [106,107] where γ ≡ Since at low energy the effective axion-pseudoscalar Lagrangian in the WSS model [52] has precisely the same form as in the Skyrme case, the CP-odd axion-nucleon couplings are expressed in terms of the nucleon mass properties (6.20) precisely as in (4.21). In terms of the model parameters we thus get (6.23) In the "chiral limit", the CP-odd meson-nucleon couplings are given in turn by precisely the same relations as in (5.20). In terms of the WSS model parameters, they are given bȳ and byḡ which is precisely the minimal scaling argued in [96] to occur in the Skyrme model. It is worth noticing that the above expressions can be equivalently obtained by considering the instanton solution for the gauge field components corresponding to the pseudoscalar mesons in the WSS model at finite θ angle [48,49].
Let us notice that, as in the Skyrme model case, the above results have been obtained by neglecting the effects of vector mesons which are expected to be subleading. 31 In appendix G we will collect analogous results for the N f = 3 case.

Discussion and numerical estimates of the couplings
Let us summarize what we have learned in the previous sections. As a first result, we have observed that the non-derivative CP-odd axion couplingsc p,n with protons and neutrons can be expressed, purely from low energy considerations valid for any number of colors N c , in two equivalent ways. The first one is in terms of the CP-odd couplingsḡ η N N andḡ πN N of the η and pions to the nucleons, as in formula (4.14) or its chiral limit (4.15). Then, by exploiting known results, the axion couplings can also be expressed in terms of the pion-nucleon sigma term and the strong contribution to the neutron-proton mass splitting, as in formula (4.21). These have to be thought of as tree level relations in the language of chiral perturbation theory. Their relevance is that the axion-nucleon couplings are expressed in a model-independent way in terms of nucleon properties which, if not directly measurable from experiments, can be computed on the lattice or can be deduced from phenomenology. There is a vast literature on lattice QCD results for the pion-nucleon sigma term. The most recent papers, from e.g. the Wuppertal-Marseille [108] and the ETMC collaborations [109] report values of δM N of 38(3)(3) and 41.6(3.8) MeV, respectively. These lattice results exhibit a certain tension with phenomenological derivations from pionic atom data (see e.g. [110]) which essentially point towards a value of the sigma term which is around 60 MeV. 32 Lattice QCD results on the strong force contribution to the neutron-proton mass splitting can be 31 This expectation appears to be confirmed by explicit computations in the complete model. We thank Lorenzo Bartolini for discussions on this point and for communication of his preliminary results. 32 Further discussions on the sigma term issue can be found in e.g. [112,113].
found in e.g. [111], giving (M n − M p ) L M ≈ 2.52 (17)(24) MeV. Using this result, the value ≡ (m d − m u )/(m d + m u ) ≈ 0.35(3) (see e.g. [13]) and the most recent ETMC computation of the sigma term we get, from our formulae (4.21), the estimates in the first column of table 1. The values are expressed in MeV times θ/f a . Instead, by employing the phenomenological extraction of the sigma term in the second paper in [110], i.e. δM N ∼ 58 (5)   Moreover, we estimated the contributions for N f = 2 mass degenerate flavors toc N (ḡ aN N ) arising through CP-violating a − π − π coupling at one-loop order in chiral perturbation theory in section 3.3, equation (3.47). 33 These contributions are smaller than the estimates for the "direct" contributions which we discussed so far, but they affect the results in the first two columns of table 1 by increasing their values by 5.5 MeV times θ/f a .
In view of the tension between lattice and phenomenological data for the sigma term, we find it interesting to write down other possible expressions for the CP-odd axion couplings. This is one of the reasons why it is interesting to have found explicit expressions for those couplings, in terms of the model parameters, in Skyrme and WSS models. Moreover, trying to improve the leading order relations we have found in chiral perturbation theory is in practice problematic, due to unknown low energy constants. In this situation, one could view the Skyrme and WSS results as providing model-dependent estimates of the latter.
Note that our computations of the CP-odd meson-nucleon couplings in both the Skyrme and WSS model have shown agreement, to leading order in the large N c limit, with formulae (4.18) and (4.20) already known from chiral effective field theory. But below we want to avoid the explicit use of the value of the sigma term and the neutron-proton mass difference. Instead, we use the Skyrme and holographic models for direct estimates of the couplings. Although, as mentioned in the Introduction, we cannot expect these models to have an accuracy at percent level (e.g. in WSS it is expected to be around 30% on average in our type of observables), at the moment we cannot claim to have a better precision with other means. A-posteriori, the numerical values obtained in these models appear to be comparable with the Chiral Lagrangian ones, once we include the one-loop corrections in the latter.
Let us begin with the Skyrme model. In the N f = 2 case the coefficient of the Skyrme term is usually fixed to e ≈ 4.84 by fitting the masses of baryons and the pion; the fit also sets f π ≈ 54 MeV [100]. With these values of parameters we obtain the estimate for the magnitude of the axion couplings (5.19) in the third column of table 1. Up to the first decimal digit, isospin breaking corrections are not visible.
One can also proceed in an alternative way as follows. Let us first notice that the numerical value ofḡ πN N in (H.5) is so small, as compared toḡ η N N (see equation (H.2)), that we can safely discard its contribution. 34 Therefore, since [66], we can rewrite the CP-odd η -nucleon coupling in (5.20), with δM N given in (5.13), to leading order in the 1 limit, as 3) The advantage of expressions (7.3) with respect to the form presented in the previous section is that they do not explicitly depend on the parameters of the Skyrme model anymore, so they are suitable to provide a straightforward estimate for the couplings. We conform the values of the physical quantities to the ones in [13], where c p , c n are calculated 35 Remembering that α ≈ 18.25 (see eq. (5.14)), using relations (7.4) and (7.3), we get the numerical values in the fourth column of table 1. Working in the "chiral limit", with the leading order mass term and at first order in increases the error in table 1, from 5 MeV to 6 MeV. 34 We do not expect an order of magnitude correction to these couplings by using slightly different values of the parameters. 35 The error on m π is irrelevant.
Coming to the Holographic WSS model, fitting the model parameters to the values for f π and m π reported in (7.4) and to the ρ-meson mass m ρ ≈ 776 MeV implies setting [61] λ ≈ 16.63 , M KK ≈ 949 MeV .
The axion non-derivative coupling (6.23), including the known quantum corrections (see appendix H), is reported in the fifth column of table 1. Moreover we can obtain an alternative estimate of the couplings in analogy with what we have done for the Skyrme model above. As we have seen, the WSS model allows for the computation of the four couplings g η N N , g πN N ,ḡ η N N ,ḡ πN Nsee formulae (6.17), (6.24), (6.25) for the N f = 2 case. At leading order in N c , such that the modulus ρ is treated as classical, and to leading order in the 1 limit, these couplings satisfy the relations where γ ≈ 1.10 was defined in (6.21) andγ ≈ 1.05 in (6.22). Therefore, from relations (4.10), (4.11), (4.15) we can express the non-derivative axion-nucleon couplingsc p,n in (6.23) in terms of the derivative ones c p,n Taking into account the values reported in (7.4) we obtain the results in the sixth column of table 1.
On top of the errors reported in table 1, we can give an estimate of the errors implied by our approximations. The estimate is going to be rough, since lacking an explicit calculation of the discarded contributions, we do not control their coefficients. Nevertheless we can get an idea of the magnitude of these errors for order one coefficients. Working in the "chiral limit" implied discarding corrections of order (m π /m η ) 2 ∼ 0.02. Moreover, the mass term in the Lagrangian is a valid approximation at leading order in m ud /Λ QCD 0.04. Finally, we have discarded corrections in 2 ≈ 0.12. 36 The three effects mentioned above increase the error ofc p from 4 MeV to 5 MeV. 37

A. The structure of axion models
In this appendix we review the general structure of axion models. An axion field may have several meanings in the literature, but for our purposes it has two properties: it has a perturbative shift symmetry and couples linearly to an instanton density (here the QCD density) of the Standard Model (SM).
Such a field was first introduced by Peccei and Quinn to solve the strong CPproblem and render the θ-angle of the strong interactions a dynamical variable.
There are several realizations of axions in QFT, with two extreme classes: the fundamental axions, [6]- [11], and the emergent or composite axions, [41]. We will produce here a typical simple model of a fundamental axions. A general prescription of theories of composite axions and further references can be found in [41].

A simple class of UV complete examples of axions can be constructed as follows.
A new sector that realizes two global U (1) symmetries is used and has two scalar field transforming under them. An anomaly free combination of the two U (1)s is gauged. The scalars obtain a vev so that they break the U (1) gauge symmetry but leave the orthogonal global U (1) symmetry unbroken. This will play the role of the PQ symmetry when coupled to the SM.
We therefore consider a set of new chiral fermions ψ i (we take them to be all left-handed) that are charged under two U (1) symmetries, U (1) 1 and U (1) 2 with charges (Q 1 i , Q 2 i ) and two Higgs-like scalars, H 1 , H 2 with charges (q 1 1 , q 2 1 ) and (q 1 2 , q 2 2 ). We gauge U (1) 1 with a corresponding gauge boson A µ and we assume that the full spectrum is anomaly free.
This implies that that guarantees the absence of the gauge anomaly and mixed anomaly in four dimensions. The bosonic part of the action is 38 We can obtain symmetry breaking where both H 1 and H 2 condense with the same vacuum expectation value for simplicity We will assume that the remaining scalar degrees of freedom, beyond the phases, have high masses and we will neglect them in the low energy dynamics. The bosonic part of the action now becomes Changing basis to we may rewrite the action as where The fieldσ will become the axion in the low energy theory and as we defined it, it has mass-dimension 1. We may normalize it as Consider now the Yukawa couplings of the fermions with the selection rules g i,j I = 0 unless q 1 We would like to have v M EW where M EW is a typical EW scale. For finite and small couplings g, g ij I , it is clear that we can take m A M EW , and that we can arrange the charges so that all the fermions ψ i acquire masses that are also well above the EW scale.
The field φ is the longitudinal component of the very massive U (1) gauge boson, and can be set to zero by a gauge transformation. The only particle from this sector that remains in the IR theory (the SM) is the massless PQ scalar σ.
The scalar σ will have direct couplings to the SM fermions (quarks and leptons) if they are charged under the two UV U (1) symmetries and if the high energy Higgs H 1,2 couple to them. The remaining global symmetry shifts σ by constants.
We have also assumed that the gauged U (1) symmetry is anomaly free. 39 However the remaining global U (1) symmetry associated with σ may have mixed anomalies with the standard model. We denote the relevant charge traces as where Q 2 is the charge generator of U (1) 2 , T a are the gauge generators of the SM and the trace is on the SM fermions. If such traces are non-zero, then in the effective field theory we have the following axionic couplings where G is the SU (3) Yang-Mills field strength and F is the U (1) electromagnetic one. We have ignored the coupling to the W-instanton density, as its effects are heavily suppressed. Therefore we can parametrize the effective couplings of σ to the SM matter as in [13] Here σ plays the role of the PQ axion for the SM. V mass−terms is a non-derivative interaction coupling the axions to mass terms of fermions and bosons. The ellipsis in (A.14) stands for higher order and higher derivative interactions that have been omitted.
We have defined the axion to be dimensionless, and its effective interaction strength is set by the axion scale f a that is assumed to be much larger than the SM scales. Therefore all axion interactions are weak. This is the reason why higher non-linear interactions are not so important in the effective action.
S axion does not contain quantum corrections associated with the SM. Due to the protected nature of the axion such corrections are very special. At weak coupling they arise because of YM instantons, that will generate a potential for the axion, [13]. This is the analogue of the U (1) A anomaly for the PQ symmetry.
The couplings of the axion to leptons and or photons are straight-forward as they can be treated in perturbation theory. However, the couplings to the strong interaction operators are trickier.
With a chiral rotation of quarks, the coupling to the instanton density can be transferred to the mass matrix of the quarks. Therefore, an invariant collection of couplings involve the quark masses and the derivative coupling to J µ . The quark component of J µ generates out of the QCD vacuum various vector mesons. Their content is axion-model dependent. The axion dependence of the quark couplings generates axion couplings to scalar mesons. Their content is also axion-model dependent.
There is however a subset of model independent couplings, if we want the axion to solve the strong CP problem. Namely c s = 0 and then there is a universal coupling of the axion to all the quark states, that amounts to a coupling to the so-called singlet η meson. It is this class of couplings that we address in this work.

A.2 Other axion models
Generically, the origin of the axion coupled to the SM may vary. In the fundamental axion case, like the simple model described in the previous subsection, the axion involves a high energy degree of freedom protected by the PQ symmetry.
The effective Lagrangian of the axion coupling it at low energies to the SM degrees of freedom is the same as in (A.14) as it is dictated by (anomalous) symmetry alone. Therefore different fundamental theories yielding axions will have a similar effective action, albeit with different charge assignments and coefficients, [10,9].
The situation is somewhat different when the axion is composite/emergent. In such cases the axion is a composite generated by the instanton density. 40 The general quadratic effective action of the emergent axion was derived in [41]. In that case there is also a mass contribution to the axion effective action originating in the UV theory. In a QCD-like UV theory this is controlled by the topological susceptibility of the hidden sector. The axion in that case has a mass with two contributions that are a priori independent: one from the hidden sector and the standard one from QCD.

B. Current terms and η -axion mixing in chiral Lagrangians
In this Appendix we address the mixing between the axion and mesons in a somewhat more general case than that described in section 2.
In general, the derivative of the axion couples to a current J µ involving contributions from leptons, vectorial quark current and the axial quark current as We write the axial current already appearing in (2.1) as J µ q5 =qQγ µ γ 5 q , (B.1) withQ a Hermitean matrix in flavor space. The vectorial current is conserved and therefore couplings to such current play no role in our analysis: we can neglect them without loss of generality. 40 Due to anomalies an analogue of the η will be as good. In general, at finite N c and N f these two mix strongly, [47]. At weak mixing it is usually the (mostly) η that is lighter in QCD-like theories but in more general theories it could be otherwise. Which of the two is lighter depends of a hidden QFT coupled at high energies to the SM.
It is however simpler to start with the action where the derivative coupling to the current J µ q5 is absent (i.e., the case already considered in Sec. 2). Therefore we first perform a (flavor dependent) axion valued chiral rotation which removes the derivative coupling. It will be reinstated using a chiral rotation in the effective theory below. After this transformation, the Lagrangian obtained from (2.1) is without the current term After this transformation, the final form of the chiral Lagrangian, corresponding to the QCD Lagrangian of (2.1) with the current term included, is where the last term is suppressed by both the flavor masses m q and 1/N c [88,115]. We show below that its effect on the axion mixing can be absorbed into redefinitions of other parameters up to quadratic terms in c aη . We also observe that the chiral rotation gives rise to kinetic mixing between the axion and the pions. Comparing the derivatives of the QCD Lagrangian and the effective Lagrangian with respect to M ij at a = 0 we obtain We then discuss how the mass eigenstates are computed to leading nontrivial order in 1/f a . First we make an observation about the kinetic mixing. We can apply a shift so that the change of the axion kinetic term exactly cancels the kinetic mixing and the remaining terms in (B.5) are unchanged up to corrections suppressed by 1/f a . The outcome is that at leading nontrivial order in 1/f a the only effect due to the O(f −1 a ) current terms, is the mixing of the physical axion with the pions implied by the above shift. Because the non-derivative couplings involving the axion in (B.5) are already suppressed by 1/f a , the non-derivative mixing terms introduced by this shift are O (1/f 2 a ), i.e., negligibly small. Therefore the effect of the kinetic mixing will be irrelevant for the analysis of the axion-nucleon couplings in Sec. 3. Higher order terms in the effective action also include a derivative coupling between the axion and η which is present even in the absence of the current (B.1) but suppressed by 1/N c [115]. The same applies to such a coupling as the kinetic mixing in (B.5): its effect is irrelevant for our studies.
We then discuss the diagonalization of the quadratic terms in (B.5). We assume that M (a) = e iaQ/fa M 0 e iaQ/fa where M 0 is independent from the axion and Hermitean. We consider the effects due to a non-Hermitean matrix below. In order to eliminate the O (1/f a ) non-derivative mixing terms with the axion (up to higher order corrections in 1/f a ), we first carry out a field redefinition √ 2η Inserting this in the action (B.5) and requiring that the non-diagonal mass terms involving the axion vanish, fixes The nondiagonal kinetic terms are then removed by a subsequent redefinition If c aη is treated as a small correction and quadratic terms in this parameter are neglected, we may write the result for the coefficients k i in a compact way in terms of the topological susceptibility (which now includes a 1/N c correction): Within this approximation, the rules (2.9)-(2.10) therefore generalize tô Notice that the last (fourth) term in (B.14) is suppressed with respect to the third term by the quark mass. The mass eigenstates are then found by diagonalizing the pion mass matrix as explained in Sec. 2. The axion mass squared is given by (B.16) We then discuss the effect of c aη for the axion-nucleon couplings studied in Sec. 3. It is sufficient to include it in the expressions (3.18) and (3.20) using the approximation in (B.12). For the former expression we find that (B.17) Since f 2 π ∼ N c the correction due to c aη is subleading in 1/N c . For the latter expression we find In this case, the correction due to c aη is suppressed by the quark masses. Since χ also contains a factor proportional to the quark masses, this correction is actually of quadratic order in the quark masses and can be safely ignored.
Notice that the last term in (B.5) also affects the pion mixing and masses even in the absence of the axion. For example, for a flavor independent quark mass m q , the masses of the pions and the η are given by Finally, we recall that for N f = 3, the diagonalization of the pion mass matrix boils down to the "standard" analysis in the literature (see, e.g. [116,77]). The η meson, the η meson, the pions and the kaons form a nonet of the vectorial U (3), which is broken to various subgroups both by 1/N c effects and by finite quark masses. Taking the light quark masses to be the same, m u = m d ≡ m ud , only the η and η mesons mix. The relevant mass matrix is given as where we dropped the terms ∝ c aη for simplicity. We adopt a convention where the physical states are given bŷ η = η p cos θ ηη − η p sin θ ηη , π 8 = η p cos θ ηη + η p sin θ ηη .
Notice that we chose the normalization ofc η similarly as that of the θ-angle:c η = O (N c ). This term could arise from the CP-violation in the CKM matrix, [55,56,57] or BSM sources. In fact, one immediately sees thatc η can be absorbed by a shift in the θ-angle (up to highly suppressed contributions) so we can set it to zero without loss of generality. We then find for the VEVs of the pions after a straightforward computation where . In order to analyze the couplings to nucleons, we also need the VEVs in Tr [χ ± ]. For consistency with chiral symmetry, in these couplings we also need to include the contributions from M 0 . We find that We also consider the terms arising from Tr [χ + ] which are linear in the mesons and the axion and are present due to the VEVs. First notice that The linear contributions are given by the second term on the second line and the terms arising from the third line after inserting the VEVs. After some cancellations, we obtain Therefore, as expected, all contributions due to M 0 can be absorbed in theθ-angle, so it is enough to consider the couplings of the mesons and the axion arising (only) from the VEVs due to the θ-angle as we do in Sec. 3.

B.2 Three-meson couplings
Having extracted the VEVs of the mesons due to a finite θ-angle, we can also extract the (CP-violating) couplings of (non-derivative) vertices with three mesons. Up to corrections suppressed at large N c , these arise from the mass term When computing these couplings it is important to note that the mixing (2.9) and in (2.10) is (for generic quark mass matrices) only consistent with chiral symmetry if cubic and higher interactions are included. That is, we write the field redefinitions of the pions and the η in general aŝ , (B.30) whereÛ was defined in (2.13). This definition differs from (2.9) and (2.10) if M 0 , Q and Π do not commute. Notice that the transformation (B.30) corresponds to an axion-valued axial transformation of the pion fields, which ensures that the field redefinition does not generate undesired higher order interactions that would be inconsistent with chiral symmetry.
After inserting the mixing with the axion from (B.30) and introducing the θ-angle through (3.17) we obtain the three pion couplings [70]: whereΠ =π a T a . The coupling of theη to the pions is In the limit of large N c and in the chiral limit with flavor independent quark masses the η -pion-pion coupling simplifies to 3/2 f f πη π aπa ≡ḡ η ππη π aπa , (B.33) where we inserted χ ≈ m 2 π f 2 π /2N f . Notice that the expressions of these couplings are identical in the WSS model and in chiral perturbation theory at large N c , so (B.33) also holds in the WSS model. The coupling between the axion and the pions is found to be C. Details on the nucleon couplings in chiral perturbation theory C.1 Non-derivative couplings from external CP-violation In this Appendix we discuss the non-derivative couplings of nucleons to the mesons and the axion which arise due to such sources of CP-violation that are not included in the θ-angle, such as the CKM matrix in the quark or lepton sector. For example, the mass terms (3.10) and (3.11) can be generalized to read The coefficients with a bar are CP violating and in the absence of the θ-angle should arise from non-QCD CP-violating sources. We only consider here CP-odd couplings which are singlets under the chiral transformations. 41 Similarly, an extra term should be included in the Lagrangian (3.13) containing flavor singlet couplings: For the CP-violating term we assumed a scaling which is analogous to the CP conserving terms, as above. We then consider the axion couplings induced by the axion terms in the QCD Lagrangian (2.1). The important terms are the first term in (C.4) and the second term in (C.5).We use notation defined in Sec. 2 and include the axion through the replacements (3.15) and (3.16).
After this, the terms of interest become where χ − (a) = uM (a) † u−u † M (a)u † . Inserting here the definition of U , and using the relations (2.9)-(2.11), we may extract the non-derivative couplings (i.e., the terms involvingN N ) of the fieldsη and a p to the nucleons.
The terms linear in a and η are 42 then Notice that since f 2 π ∼ N c , both axion terms are of the same order in 1/N c . The η term in (C.8) is suppressed by a factor of N f /N c with respect to the term in (C.7) so it is negligible in the 't Hooft limit but should be included in the Veneziano limit.

C.2 CP-odd nucleon-meson couplings
Here we extend the discussion of the nucleon couplings to mesons in Sec. 3 to include in particular the contributions from the single trace terms in (C.1). As it turns out, for this it is convenient to also redefine the nucleon fields. First we define the transition matrix K B through (see (3.2 Then the redefined baryon fieldsB are obtained by acting on B with K B as in (3.3).
Only the couplings involving χ + in (C.1) are affected by the baryon redefinition at the order we are interested in. The combination coupling to the redefined baryons is where we also noted that a = a p up to subleading corrections in 1/f a which are irrelevant here. After introducing the θ-angle through (3.17), we obtain the CP-odd couplings C.3 Group theoretic structure of the nucleon couplings to the mesons We discuss here the structure of the single trace terms including sum over contractions, which appear in Eqs. (3.7), (3.10), and (3.11). We start by considering the structure at N f = 2 as an example. The number of terms in the sums can be understood in the following way, which can be generalized to higher N f . First, we classify the pions and the η according to the underlying quark representations 43 1⊕3 = 2⊗2. All possible flavor structures of the η −N − N and π −N − N couplings are then obtained by adding the nucleons, and consequently by identifying the singlets in the product state. We write When expanding the last expression, only products of a representation with its conjugate contain singlets (and exactly one each), 1 ⊗1 = 1 and 3 ⊗3 = 1 ⊕ 3 ⊕ 5 .
(C. 15) These two singlets map to the two terms in (3.26) and in (3.27). We proceed by commenting on the couplings at N f = 3. We also first fix N c = 3 for simplicity. Then, in addition to the nonet of light mesons, the ground state of spin-1/2 baryons transforms under the adjoint representation of SU (3). The usual way [82] to write the couplings corresponding to (3.5) and to (3.7) iŝ where g a = −D/2 +ĝ a and with the baryon octet included in the traceless 3 × 3 matrix transforming as B → KBK † . We renamed the couplings g (1) a and g (2) a as D and F following standard conventions.
In order to analyze the number of terms, we may again pack the mesons into 3 ⊗3 and consider the product of the nucleon and quark representations: 8 ⊗ 3 = 3 ⊕ 6 ⊕ 15. When multiplied with the similar antiquark representations each of the representations 3, 6, and 15 again gives rise to a singlet. Therefore we obtain three singlets which map to the three terms of (C. 16) and (C.17). The same argument actually works for the ground state spin-1/2 baryons for any N f and any odd N c . In this case the baryons belong to a representation defined through a Young diagram having exactly two rows, one with N c /2 + 1/2 boxes and the other with N c /2 − 1/2 boxes [85,86]. For excited states of baryons in generic representations one obtains O (N c ) independent terms.

D. On-shell nucleon vertices in the limit of large mass
We first consider the bilinears of Dirac spinors in the limit of large nucleon mass m N . For leading contributions in this limit we expect that the exchanged three-momenta are ∼ m 0 N and the variations in the nucleon energies are ∼ 1/m N . In general the Dirac propagator can be decomposed as The second term corresponds to a nucleon propagating backward in time and is suppressed in the limit of large m N (assuming that we are discussing a nucleon state rather than an anti-nucleon).
This suggests that up to corrections suppressed by 1/m N the nucleon-meson interactions are determined by the (on-shell) amplitudesū(p 2 , s 2 )Γu(p 1 , s 1 ) where Γ = I, γ 5 , γ µ , and γ µ γ 5 . We take the incoming nucleon to be in the rest frame with momentum p 1 = (m N , 0) and write the outgoing nucleon momentum as p 2 = (m N , k) where we dropped terms ∼ 1/m N in the energies. We use a representation where the spinors are given explicitly by (D.2) where a possible basis for the spin wave functions χ(s) is given by the spinors (1,0) and (0, 1), but we leave the basis unspecified. By direct computation we find that where χ i = χ(s i ). We notice the following identities which hold for any value of m N . Notice also that because v(p, s) = γ 5 u(p, s) the nucleon creation and annihilation amplitudes immediately follow from the above results. These identities suggest that the operators ∂ µ η N γ µ γ 5 N ∝ ∂ µ η N S µ N , and 2m N η N γ 5 N , where S µ is the nucleon spin operator, are identical at least to leading order in the limit of large m N . The spin S µ is fixed at least to leading order as the amplitudes do not flip spin in a basis aligned with k: the transverse components of the amplitude (D.6) contain also potentially spin-flipping elements but these vanish when contracted with the only available momentum k.
We also notice that at O (1/m N ) the couplings involving γ 5 above are ∝ k. This seems to be due to parity: whenθ = 0, the only available parity-odd scalar is k · S. 44 Notice however that e.g. the nucleon pair creation amplitude ∝ū(p 2 , s 2 )γ 5 v(p 1 , s 1 ) is not proportional to k · S. E. CP-odd couplings in the Skyrme model with N f = 3 In this section we will provide a derivation of the CP-odd couplings of the axion and of the η, η mesons to the nucleons, in the Skyrme model with N f = 3 non-degenerate 44 There are two independent momenta, but invariance of the amplitudes under O (1/m N ) boosts requires that at linear order only momentum differences (i.e. k) appear. Notice also that S is only affected by O (1/m N ) corrections under such boosts. flavors. We use the notation in [101], where the Lagrangian, once the substitution θ → a/f a is performed, is already in a diagonal basis for the axion. Hence, for the rest of the discussion we drop the subscript "p" of the fields for notational convenience. Moreover, in this section we work in the "chiral limit" m π m W V , m K (m W V is defined in (5.4) and m K is the kaon mass).
Actually, the Lagrangian in [101] is written at linear order in θ. Generalizing that expression to the non-linear case, the mass term reads .

(E.2)
In this expression M 1 = m 2 π , M 3 = 2m 2 K − m 2 π andγ (called α in [101]) is a certain combination of the masses and the topological susceptibility which will be irrelevant in the "chiral limit". Finally λ θ , proportional to θ, will be given in a few lines. The leading effect of SU (2) isospin breaking will be introduced later on.
The Skyrmion is just the embedding of the SU (2) one into SU (3), so the configuration we need is 45 U = U s = e i λ·xF (r) , where with λ we denote the vector composed by the first three Gell-Mann matrices. The function F is the same as for the N f = 2 model reported in section 5.1. One can calculate the quark contribution to the nucleon mass by setting λ θ = 0, which corresponds to setting to zero the axion field, and putting (E.1) on-shell on the Skyrmion solution. The result is where α is given in (5.14). Note that in this formula we have not taken the chiral limit. In the degenerate mass case, eq. (E.4) exactly reduces to the N f = 2 expression (5.13).
In order to extract the axion couplings we expand (E.1), on-shell on the Skyrmion, up to second order in λ θ . In the "chiral limit", the result turns out to be very simple In the same limit, we have [101] λ θ = 1 2 βm 2 π θ , (E.6) 45 In this expression we have dropped a factor Exp[− i  From the results in [101] one can as well extract the CP-breaking η and η couplings to the nucleons. In fact, the relevant Lagrangian term in the "chiral limit" reads 46 At the order we are working, it is sufficient to consider the η − π 8 mixing in (B.21), where the mixing angle is θ ηη ≈ −π/9 (see e.g. [101]). Using these expressions in (E.10) we obtainḡ Therefore, apart from the last result, which takes into account the η − η mixing, the N f = 3 case in the "chiral limit" gives the same formal outcome of the N f = 2 model for the axion couplings. This suggests that, not unexpectedly, the corrections due to the inclusion of the strange quark are small.

F. Instanton moduli quantization
We remind the reader that the quantum description of the instanton solution in the WSS model, is obtained in analogy with what is done for the Skyrmion [66]. Since M B0 = 8π 2 κ ∝ λN c 1, a non relativistic limit can be taken in which the system reduces to a quantum mechanical model for the time-dependent instanton (pseudo) moduli. The slowly rotating instanton solution around z ≈ 0 is given by U cl → a(t)U cl a(t) † , (F.1) 46 Note the different convention on the sign of θ with respect to [101].
where, setting X = 0 we have Here, α = α(r) = π 1 + ρ 2 /r 2 , (F.3) with r = | x|. Time-dependent rotations act non-trivially on the non-Abelian components of the gauge field. Consistency with the equations of motion, in turn, induces a non trivial time dependence on the Abelian components too, except that on A 0 [64].
In the z 1 region, the time dependent non-Abelian field components read The Abelian field components are given in turn by and The quantum Hamiltonian for the baryon states is obtained by substituting the above time-dependent solution into the WSS action.
Notice that the structure of the instanton solution above recalls that of the Skyrme model with vector mesons examined in [89]. There, the time-dependent soliton solution includes non-trivial expressions for the time component of the isotriplet (ρ 0 ), the space components of the iso-singlet vector (ω i ) and the η-type meson. The contribution of the latter, which turns out to be proportional to the angular velocity of the spinning soliton, turns out to be crucial for many instances. In [107] it has been used to compute the mass splitting between neutron and proton in the WSS model with non-degenerate quark masses. It is crucial for providing non zero values for both the CP even coupling g η N N and the CP-oddḡ πN N one. In turn, these couplings are related, using the mixing terms, to further axion-nucleon couplings.
In the present setup, the angular velocity is given by where J j is the spin operator J k = i 2 −y 4 ∂ ∂y k + y k ∂ ∂y 4 − klm y l ∂ ∂y m , (F. 8) for the axion content of the η, η . Here and in the following, we have made use of the GMOR relations [117] m 2 π f 2 π = 2c(m u + m d ) , m 2 K ± f 2 π = 2c(m u + m s ) , m 2 K 0 f 2 π = 2c(m d + m where the dots stand for a term in x 3 sin α(r) which is vanishing once we integrate the result in space. Up to the latter term, the result (G.7) is exactly the same one of the N f = 2 computation, so the couplings are still the ones in section 6.3. This is consistent with what we found in the Skyrme model. As it happened there, the "chiral limit" m π m W V , m K renders subleading the corrections due to the third quark flavor.

H. Alternative numerical estimates of the couplings
In this Appendix we collect some alternative numerical estimates of the couplings calculated in the main body of the paper in the Skyrme and WSS models, using a different parameterization with respect to the one employed in section 7. Here we plug in the formulae for the couplings the standard values of the parameters usually employed in the two theories we have considered.
Let us begin from the Skyrme model. For the N f = 2 case the coefficient of the Skyrme term is usually fixed to e ∼ 4.84 by fitting the masses of baryons and the pion; the fit also sets f π ∼ 54 MeV [100]. Therefore, with these values of parameters the estimate for the magnitude of the axion couplings (5.19) in the mass degenerate case of the N f = 2 Skyrme model is For what concernsḡ πN N , the standard choice of parameters in the extended model of [89,98], to which we refer for the notations, is It is worth noticing that a first estimate of the same coefficient, in chiral perturbation theory with N f = 3, gave |g πN N | ≈ 0.027|θ| [70]. The parameters extracted from hadron mass fittings in the N f = 2 + 1 case are instead [101] e = 3.87 , f π = 44.5 MeV , (H.6) so, remembering that α ≈ 18.25 and using x = 0.48 [13], the numerical estimate of the axion coupling (E. 8) in the N f = 2 + 1 Skyrme model is (H.7) The sizable numerical difference between (H.1) and (H.7) is due to the difference in the two cases of the values of the parameters f π , e usually employed to fit nuclear data. This is a known drawback of the Skyrme model. Therefore, in the formulae we replace ρ n cl with ρ n = ρ 3+n R(ρ) 2 dρ ρ 3 R(ρ) 2 dρ . (H.13) Setting N c = 3 and fixing the parameters as above one finds for the quantum corrected couplingsĝ A ≈ 0.527 , g A ≈ 0.734 . (H.14) Concerning the axion non-derivative coupling (6.23) in the mass-degenerate case, by keeping ρ to its classical value we obtain the estimatē taking into account the quantization of the instanton moduli.