QCD θ-vacuum energy and axion properties

At low energies, the strong interaction is governed by the Goldstone bosons associated with the spontaneous chiral symmetry breaking, which can be systematically described by chiral perturbation theory. In this paper, we apply this theory to study the θ-vacuum energy density and hence the QCD axion potential up to next-to-leading order with N non-degenerate quark masses. By setting N = 3, we then derive the axion mass, self-coupling, topological susceptibility and the normalized fourth cumulant both analytically and numerically, taking the strong isospin breaking effects into account. In addition, the model-independent part of the axion-photon coupling, which is important for axion search experiments, is also extracted from the chiral Lagrangian supplemented with the anomalous terms up to Op6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({p}^6\right) $$\end{document}.


Introduction
A CP-violating topological term, i.e., the θ-term, is allowed in the Quantum Chromodynamics (QCD) Lagrangian. It can be written as where α s is the QCD coupling constant, G µν,c is the gluon field strength tensor, with c a color index, andG c µν = ε µνρσ G ρσ,c /2 its dual. Because none of the quarks is massless, physical observables only depend on a combination of the θ 0 parameter and the phases present in the quark mass matrix M q , i.e., θ = θ 0 + arg detM q . Being a dimensionless parameter, the natural value of θ is expected to be O(1), which would significantly affect physical systems such as atomic nuclei, and lead to measurable effects, as nucleons, for instance, would possess a nonvanishing electric dipole moment [1]. However, the so-far negative results of experimental searches for the nucleon electric dipole moment lead to a tiny upper limit: |θ| 10 −10 [2][3][4][5][6][7]. To understand why the value of θ is so small is the socalled strong CP problem. One elegant possible solution of this problem is the Peccei-Quinn (PQ) mechanism [8,9], which introduces a global U(1) symmetry, called PQ symmetry. This symmetry is spontaneously broken at energies much higher than the typical QCD scale of order O(1 GeV) and is also broken by an anomalous coupling to gluon fields. The axion appears as the corresponding Goldstone boson [10,11] which has an anomalous coupling to GG. The parameter θ is then dynamically driven to zero at the minimum of the axion potential, giving rise to a possible solution to the strong CP problem.

JHEP05(2020)001
Some important quantities in axion physics, such as the axion mass and self-coupling, are dictated by the axion potential. The visible axion models [10,11] with the axion decay constant at the electroweak scale or even smaller are believed to have been ruled out by experiments. For the invisible axion [25][26][27][28], its mass window is usually assumed in the range from about 10 −6 eV to 10 −2 eV. According to constraints from astrophysical observations, the present bounds on the axion decay constant is 10 9 GeV f a 10 12 GeV [29,30] (we refer to refs. [24,[31][32][33][34] for several recent reviews). 1 Within the above available parameter space, the axion may be the main source of cold dark matter in the universe [36][37][38][39][40][41]. In addition, it may form a Bose-Einstein condensate [42] or even compact boson stars [43][44][45][46][47][48][49]. The axion can couple to the Standard Model (SM) particles like electrons, nucleons, photons and so on. However, all these couplings are suppressed by the axion decay constant f a , which is remarkably large, resulting in the invisible axion which has very weak couplings to the SM particles [36]. Since the axion-photon coupling vertex, see eq. (4.1) below, allows for the production of an axion from the interaction of a photon with a background magnetic field, the axion-photon coupling g aγγ plays a central role in axion searches in both laboratory experiments and stellar objects [24]. In this case it is very useful to study the axion properties, especially the axion-photon coupling, at a high precision from the theoretical point of view.
At low energies in QCD, all hadronic degrees of freedom are frozen and thus can be neglected except for the pseudo-Goldstone bosons of the spontaneous chiral symmetry breaking. Chiral perturbation theory (CHPT) [50][51][52], as the low-energy effective theory of QCD, can be used to describe the vacuum properties as well as the dynamics of QCD in the non-perturbative regime reliably. In this paper, we will calculate the θ-vacuum energy density, or equivalently the QCD axion potential, up-to-and-including next-to-leading order (NLO) in SU(N ) CHPT. Setting N = 3, the mass and self-couplings of the axion can then be extracted from a Taylor expansion of the axion potential. In addition, we also compute the NLO corrections to the axion-photon coupling.
Before continuing, we would like to stress that a similar study was performed in ref. [53], where the QCD axion potential derived in two-flavor CHPT up to NLO (the QCD θvacuum energy density up to NLO was first derived in ref. [54]) is used, and a matching between two-flavor and three-flavor CHPT is performed to determine the axion-photon coupling. Here, the calculations are explicitly done in SU(N ) CHPT for the θ-vacuum energy density and with N = 3 for the other quantities. In SU(3) CHPT, the topological susceptibility as well as the fourth cumulant of the topological distribution up to NLO have been calculated before using the Goldstone boson masses at θ = 0 [55][56][57]. Very recently, the topological susceptibility and axion mass are calculated up to next-to-next-to-leading order and including electromagnetic corrections up to O (α em ) in SU(2) CHPT [58]. The axion-nucleon coupling is also calculated up to the leading one-loop order in ref. [59]. Here, we derive the one-loop contribution to the SU(N ) θ-vacuum energy density by a direct calculation of the logarithm of the functional determinant for the Goldstone bosons in a 1 It was recently argued that there is still a possibility for a viable QCD axion model with a mass in the MeV range [35].

JHEP05(2020)001
θ-vacuum, extending the two-flavor treatment in ref. [54] to the case of N non-degenerate flavors. 2 This study is useful when the up and down quark masses take values close to the strange quark one. This could happen in lattice QCD calculations where the quark masses are parameters that can be chosen freely.
The outline of the paper is as follows. In section 2, we generalize the calculation of the θ-vacuum energy density in the framework of SU(2) CHPT [54] to the SU(N ) case with N non-degenerate quark masses. In section 3, we derive the axion properties, including the mass and self-coupling, in detail for N = 3. In section 4, the model-independent part of the axion-photon coupling is determined from the chiral Lagrangian supplemented with the odd-intrinsic-parity sector of the chiral effective Lagrangian. Section 5 contains a brief summary and discussion. The appendix provides a relatively detailed derivation of the recursion relation giving rise to the general solution of the vacuum angles φ f .

θ-vacuum energy density up to NLO
The QCD axion potential as a function of a/f a has the same form as the QCD θ-vacuum energy density as a function of θ. In this section, we compute the θ-vacuum energy density in SU(N ) CHPT with N non-degenerate quark masses, which is an extension of ref. [54], where the θ-vacuum energy is computed up to NLO in the SU(2) and SU(N )-symmetric cases.

Leading order
The discovery of instantons not only solved the U(1) A problem, but also implied that there is a θ-term in the QCD Lagrangian. In order to study the physics with a θ parameter, it is common to rotate away the θ-term by performing a chiral rotation on the quark fields. At low energies, we can then match the resulting Lagrangian to the chiral Lagrangian since now the relevant degrees of freedom are the pseudo-Goldstone bosons [51,52]. The Lagrangian density of the SU(N ) CHPT at leading order (LO) in a θ-vacuum is where χ θ = 2B 0 M q exp[iX a θ] contains the θ angle and the diagonal and real quark mass matrix is M q = diag{m 1 , m 2 , . . . , m N }, and . . . denotes the trace in the flavor space.
The matrix X a takes the following general form: which arises from a U(1) A chiral rotation on the quark fields eliminating the θ-term in the QCD Lagrangian. In this case, the θ-dependence is completely captured by the quark mass term. The U(1) A chiral rotation can be distributed to different quark flavors, leading to different choices of X a . F 0 is the pion decay constant in the three-flavor chiral limit, and B 0 = − qq /F 2 0 is related to the scalar quark condensate. U(x) is the field configuration for the vacuum and the Goldstone bosons of the spontaneous breaking of chiral symmetry.

JHEP05(2020)001
It can be written as U(x) = U 0Ũ (x), whereŨ (x) collects the Goldstone bosons, and U 0 describes the vacuum, parameterized as [55,61]. For the SU(3) case,Ũ = e iΦ/F 0 , with Φ given by Note that the neutral flavor eigenstates in the octet of the pseudoscalar mesons as shown above, i.e. π 3 and η 8 , are not mass eigenstates. Diagonalizing the mass matrix of the meson fields, one gets the physical mass eigenstates π 0 and η, which are mixtures of π 3 and η 8 . By expanding the LO Lagrangian in terms of the meson fields to quadratic order, the LO θ-dependent meson masses including isospin breaking effects are obtained as where for convenience we have defined (m 1 , m 2 , m 3 ) ≡ (m u , m d , m s ) and The parameter ξ is given by with θ the π 0 -η mixing angle in the θ-vacuum, which arises due to strong isospin breaking. Diagonalization of the mass matrix requires Obviously, the above θ-dependent Goldstone boson masses reduce to the standard SU(3) relations [52] by taking the limit θ = 0 and setting φ f = 0. The dependence of φ f on the θ angle needs to be determined by minimizing the vacuum energy to be discussed below.
To determine the ground state, i.e. the vacuum, we setŨ = 1. Performing the trace in eq. (2.1), one obtains the LO potential energy density

JHEP05(2020)001
Moreover, minimizing eq. (2.9) with respect to the parameters φ f with the constraint f φ f = θ gives the following equations 3 for SU (3), and similar equations for SU(N ), i.e., m f sin φ f is the same for all flavors. The above equations depend only on the linear combination φ f given in eq. (2.6), instead of on X f and ϕ f separately. This implies that φ f is physical while X f and ϕ f are not. One can use this freedom to choose the "gauge" most convenient for the question of interest. One possible choice is to choose X a = 1/N, which is commonly used in the literature (see, e.g., refs. [54,55,[61][62][63]). Noticing that the only constraint on X a is X a = 1, one may also choose the U(1) A rotations to be to simultaneously shift the θ angle to the quark mass matrix phase and align the vacuum properly. This is a convenient choice for the aγγ coupling (with θ changed to the dynamical axion field a/f a ) to be discussed in section 4 since this removes the leading order a-π 0 and a-η mixing. The equations (2.10) do not admit an analytical solution in terms of elementary functions in a compact form. 4 In the isospin symmetric case, the up and down quark masses are degenerate m 1 = m 2 ≡ m but m = m 3 , we have φ 1 = φ 2 ≡ φ, and then eq. (2.10) becomes [64] m sin φ = m 3 sin(θ − 2φ), (2.12) which allows for analytic solutions, though complicated ones. If one focuses on the cumulants of the QCD topological distribution, which are derivatives of the vacuum energy density, e vac (θ), with respect to θ, , n ∈ N, (2.13) one may solve eqs. (2.10) by expanding in powers of θ. Specifically, c 2 corresponds to the topological susceptibility. Up to O θ 3 , one gets [55] φ f =m where we have introduced For the vacuum alignment in SU(2) CHPT up to NLO, we refer to the appendix of ref. [62], which also shows that it is sufficient to consider the LO vacuum alignment for the computation of the cumulants up to O p 4 . 4 In the SU(2) case, there is an analytic solution [61], which then allows to derive a closed form of the vacuum energy density up to NLO in the chiral expansion [54].

JHEP05(2020)001
with i running over all the flavor indices considered in the theory. The solutions in eq. (2.14) are not restricted to the three-flavor case but also valid for N > 3. Consequently, the θ-dependence of the vacuum energy density at LO can be obtained by substituting the solution in eq. (2.14) into eq. (2.9), which gives [55] In appendix A, we work out a recursion relation for φ f up to an arbitrary power of θ, where k j are non-negative integers, K t ≡ t j=1 k j , (k 1 ,...,kt) means that the sum runs over all possibilities of k j satisfying In the next subsection, we will compute the one-loop contribution of the Goldstone bosons to the energy density.

Next-to-leading order
To study the θ-vacuum energy up to the NLO, O(p 4 ), one has to include both the treelevel diagrams from L (4) and the one-loop diagrams with innsertions from L (2) . The SU(N ) chiral Lagrangian at NLO is given by where we only display the terms relevant for the vacuum energy. The L i and H 2 are the so-called low-energy constants (LECs) and the high-energy constant (HEC), respectively. The latter is only required for renormalization and does not appear in observables. After setting U = U 0 and evaluating the traces, one gets the tree-level contribution to the NLO vacuum energy density (2.20) The LECs and HEC contain both ultraviolet (UV) finite and divergent parts. They are related to the renormalized ones, denoted by an upper index r, by [52,65] the UV divergence at the space-time dimension d = 4, where µ is the scale of dimensional regularization. The UV divergence in the NLO tree-level contribution exactly cancels the one arising in the one-loop contribution, as will be seen below. Now let us calculate the one-loop contribution to the θ-vacuum energy density. In the classical CHPT papers [51,52], the one-loop effective generating functional is expanded around the free-field configuration at θ = 0. This treatment is then applied to derive the topological susceptibility and the fourth cumulant in SU(N ) CHPT in refs. [55,56,63]. The expression for the vacuum energy density at NLO in SU(2) with non-degenerate quark masses, as well as that in SU(N ) with degenerate quark masses, is derived in ref. [54], where the generating functional is expanded around the free-field configuration in the θvacuum. The result allows for an evaluation of any cumulant of the QCD topological charge distribution, and is the QCD axion potential at NLO [53]. Here, we generalize the result in ref. [54] to SU(N ), with N non-degenerate quark masses. The effective action for the free-field configuration in the θ-vacuum is where "Tr" denotes traces over both the flavor (in the adjoint representation) and the coordinate spaces, and the differential operator D 0 (θ) takes the following form where P, Y = 1, . . . , N 2 − 1 are the flavor indices of the Goldstone bosons, andM P (θ) are θ-dependent meson masses at LO given in eq. (2.5). Within dimensional regularization, one gets the one-loop contribution to the vacuum energy density as [54] e (4,loop) where V is the space-time volume, the P runs over the Goldstone boson mass eigenstates (for the SU(3) case, they are given in eq. (2.5)), and the term proportional to λ collects all the UV divergences in the one-loop contribution.
Noticing that the matrix elements of the diagonalized mass-squared matrix of the Goldstone bosons are given by where we have used the SU(N ) version of eq. (2.10) to replace all m i sin φ i by m 1 sin φ 1 , and have neglected the θ-independent terms. From the above θ-vacuum energy density, the lowest two cumulants of the topological charge distribution up to NLO can then be easily extracted. It can be checked from eq. (2.29) that we can reproduce the expression of topological susceptibility at NLO keeping all orders in strong isospin breaking exactly given in ref. [66]. We are more interested in the axion mass and its self-coupling, and thus we will extract them from the axion potential based on the relation between the θ-vacuum energy and axion potential in the following section. Numerical values of the topological susceptibility and the normalized fourth cumulant will also be given for reference.

Axion mass and self-coupling
Both the axion mass and self-coupling are important quantities, since they directly affect experimental searches for the axion. For example, one tries to detect the axion in microwave cavities by stimulating their conversion to photons via the Primakoff effect within an external magnetic field [24]. The axion self-coupling plays an important role in the formation of an axion Bose-Einstein condensation [42] as well as possible boson stars [43-45, 47, 67]. This motivates the study of these two quantities in this section to high precision. Before we proceed to derive the axion mass and self-coupling up to NLO, let us discuss a little bit about the axion solution to the strong CP problem, and start with the effective Lagrangian,

JHEP05(2020)001
where in addition to the θ-term, a pseudoscalar axion field is introduced which couples to gluons. As shown by Peccei and Quinn [8,9], the periodicity of the vacuum expectation value (VEV) GG in θ+a/f a forces the minimum of the axion VEV to be at θ+ a /f a = 0, and thus the θ-dependence is eliminated. Expanding the axion field around its VEV, one sees that the θ-vacuum energy density derived in the previous section, with θ being replaced by a phys /f a , gives the axion potential, where a phys = a − a is the physical axion field. In the following we will denote a phys as a for simplicity, and then the axion potential is given by V (a) = e vac (a/f a ).
Expanding V (a) in powers of the axion field around the vacuum, we obtain Comparing the above equation with the definition of cumulants of the QCD topological distribution in eq. (2.13), one finds the following relations for the axion mass and axion self couplings: where c 2n are the cumulants defined in eq. (2.13) with n ≥ 2. Thus, the axion mass and four-axion self-coupling at LO are given by respectively, wherem = (m u + m d )/2, and we have replaced B 0 and F 0 by M 2 π + /(2m) and F π , the physical pion mass squared and decay constant, respectively, which is legitimate at LO. One sees that at LO, the difference between the SU(3) and SU(2) expressions resides merely in the definitions ofm andm [3] in eq. (2.15).
In the same way we have calculated the axion mass and self-couplings at LO. Their expressions at NLO, including the higher order corrections, can be extracted from eq. (2.29). The former reads where we have used the NLO expressions for the pion mass and decay constant [52]:

JHEP05(2020)001
Similarly the self-coupling up to NLO can be easily obtained as The numerical evaluation requires the values of the quark mass ratios and of the LECs, which have been determined by the lattice QCD calculations and experimental data. A review of the present knowledge of the LECs appearing in the chiral Lagrangian for the meson sector can be found in ref. [68]. Using the input values listed in table 1, we find the axion mass and the quartic axion self-coupling at NLO to be m a = 5.89(10)µeV · 10 12 GeV f a , (3.8) respectively. Here we have used the charged pion mass in eq. (3.6) for eliminating the overall B 0 (m u +m d ) factor in m 2 a and λ 4 . Although the difference between the charged and neutral pions from QCD is of O δ 2 , the charged pion receives an electromagnetic contribution at LO. Such an effect to the quantities of interest here can be eliminated if using the neutral pion mass instead, which amounts to replacing M 2 π + by M 2 π 0 in eqs. (3.5) and (3.7) and adding the following terms inside the curly brackets of these two expressions [52]:  Table 1. Numerical inputs used in this paper. The pion decay constant F π , and experimental meson masses M P are in units of MeV, and are taken from ref. [33]. The renormalized LECs L r i are in units of 10 −3 ; they correspond to values at scale µ = 770 MeV and are taken from ref. [68]. The NNLO anomalous LECs C W 7 and C W 8 are given in units of 10 −3 GeV −2 ; for their determinations, see the text. For the quark mass ratios defined as z = m u /m d and r = m s /m, we take the FLAG average of the N f = 2 + 1 lattice results [70].
where the electromagnetic effects have been taken into account. As a result, the values in eqs. (3.8) and (3.9) become m a = 5.71(9)µeV · 10 12 GeV f a , (3.12) which are regarded as our results for these quantities and will be used in the following. As we mentioned earlier, both the axion mass and its self-coupling are tightly related to the cumulants of the QCD topological charge distribution through the θ-vacuum energy density, see eq. (3.3). Thus, from eq. (2.13) or (3.3) we can further extract the numerical values of the topological susceptibility χ t and the normalized fourth cumulant b 2 = c 4 /(12χ t ) [56] with the inclusion of isospin breaking effects at zero temperature, i.e.,  [69]. This indicates that the explicit inclusion of the strange quark degree of freedom does not induce large differences on the axion properties. There are at least two compelling reasons accounting for this feature. First, the effects from the heavier quark flavors have been largely included in the corresponding SU(2) LECs. Second, in ref. [53] the authors performed their numerical calculations with a matching between two-flavor and three-flavor CHPT LECs. Thus, the inclusion of the strange-quark degree of freedom does not change the results

JHEP05(2020)001
sizeably. Yet, the expressions given here should be useful for chiral extrapolation of lattice results performed at unphysical quark masses, in particular when the up and down quark masses are close to the strange quark one.

Axion-photon coupling
The axion-photon coupling is defined by the following Lagrangian (see, e.g., refs. [53,71,72]), whereF µν = 1 2 ε µνρσ F ρσ , with F µν the electromagnetic field tensor with the sign convention 0123 = +1. Specifically, the axion-photon coupling is given by where E/C is the ratio of the electromagnetic and color anomaly coefficients, which is given by n (Q PQ Q 2 )/ n (Q PQ T 2 ), with the sums running over all fermions with PQ charges Q PQ , and T a the QCD color generators satisfying T a T b = T 2 δ ab /2. The value of E/C depends on the specific axion models. The first term in g QCD aγγ is the contribution from the axial rotation of the quark fields, q → exp i a 2fa X a γ 5 q with X a = 1 (here we use the convention γ 5 = iγ 0 γ 1 γ 2 γ 3 ), which was introduced to eliminate the term a fa αs 8π G c µνG c,µν from the axion Lagrangian. The second term in g QCD aγγ , g mix aγγ , is the contribution from the a-π 0 and a-η mixings, with the π 0 and η coupled to two photons.
As discussed below eq. (2.10), there is a freedom of choosing the diagonal matrix X a satisfying X a = 1. If it is chosen as X a = diag {m/m u ,m/m d ,m/m s } =mM −1 q as in refs. [29,73], then U =Ũ = e iΦ/F 0 , see eq. (2.11), and there is no a-π 0 or a-η mixing term in the LO chiral Lagrangian. One obtains the O p 4 contribution to the model-independent aγγ coupling to be g QCD,(4) This result recovers the one derived in SU(2) CHPT [53] at O(p 4 ) in the limit of m s → ∞.
The same result can also be obtained by using other choices of X a . In that case, one needs to consider a-meson mixing. The Wess-Zumino-Witten (WZW) Lagrangian [74,75] with an external photon field can be used to get the mixing contribution. The Lagrangian is given by [76][77][78] where e > 0 is the electric charge unit, Q and N c denote the usual diagonal quark charge matrix, Q = diag{2/3, −1/3, −1/3} for the three-flavor case, and the number of quark JHEP05(2020)001 colors, respectively. Here the convention is such that U transforms under SU(3) L ×SU(3) R as U → g R U g † L with g L and g R elements in SU(3) L and SU(3) R , respectively. According to Weinberg's power counting scheme, the above WZW Lagrangian starts to contribute from O(p 4 ). The axion-meson mixing contribution can be obtained by substituting U in the above Lagrangian by exp −iY a a fa with Y a = X a −mM −1 q . One finds Using eq. (4.2), one again gets the expression given in eq. (4.3).
Our goal in this section is to compute the axion-photon coupling to O(p 6 ). The chiral Lagrangian with a minimal set of terms in the anomalous-parity strong sector at O(p 6 ) has been given in ref. [79], not only for SU(2) but also for SU(N ) with N ≥ 3. Based on the anomalous Lagrangians, several works have been done in the anomalous-parity sector [80,81]. In this work, only the terms proportional to C W 7 and C W 8 are relevant to the axion-photon coupling, which read where C W 7 and C W 8 are two LECs. We have taken the same notation as in ref. [79]. In the following, we choose X f =m/m f and U =Ũ for the computation of the aγγ coupling. With this convention the diagrams relevant for the computation of the O(p 6 ) corrections to g aγγ are depicted in figure 1: (a) the axion-pion and axion-eta mass mixing from the NLO tree-level Lagrangian; (b) the tree-level diagram from L (6) ano ; (c) one-loop diagrams with one vertex taken from L WZW and the other one taken from the LO chiral Lagrangian; the contributions from diagrams (d) and (e) exactly cancel with each other with the upper photon line in diagram (d) being on shell. It is interesting to note that for the anomalous processes such as π 0 , η and η decaying into two photons, the one-loop contributions vanish when the up-down quark mass difference is neglected [82,83]. Likewise, in the SU(2) case the sum of all one loop corrections vanishes when both photons in the final state are on-shell [53]. However, in the SU(3) case, diagram (c) does contribute to the axion-photon coupling at O(p 6 ) when taking isospin breaking effects into account. Note that the pion-eta mixing needs to be considered in order to keep g aγγ scale-independent and UV finite.
Putting together all the pieces, we obtain the axion-photon coupling keeping all orders in strong isospin breaking up to O(p 6 ) as JHEP05(2020)001 where The functions f ± (sin, cos) are equivalent to f ± (cos, sin) with the sine and the cosine interchanged, i.e., f ± (sin, cos) = f ± (cos → sin, sin → cos), and is the LO pion-eta mixing angle in the vacuum as can be obtained by setting θ = 0 in the expression of θ in eq. (2.8).
For the parameters C W 7 and C W 8 , it was argued in ref. [80] that C W 7 is largely suppressed compared to C W 8 as the latter receives a strong contribution from the η while the former does not. The authors also suggested |C W |T η | 2 with the η → γγ amplitude given by [80] to extract the value of C W 8 from the measured value of the η → γγ width: (0.516 ± 0.020) keV [33]. Following ref. [80], we take F η = (118.4 ± 8.0) MeV and assign a 30% uncertainty for the O m 2 s contribution compared to that of O (m s ), we get C W 8 = (0.60 ± 0.20) × 10 −3 GeV −2 , as listed in table 1. We have set C W 7 to 0 as its effect can be absorbed into the uncertainty of C W 8 . With the input parameters presented in table 1, one gets The error for the axion-photon coupling is also dominated by the uncertainties of C W 8 , r and L r 7 , which are of similar size. From eq. (4.10), we obtain g aγγ 1.2 × 10 −16 . . . 1.2 × 10 −13 GeV −1 for the axion mass in the range 1 . . . 1000 µeV with E/C = 8/3. Especially for m a = 6.7 µeV, this equation predicts g aγγ 8.1 × 10 −16 GeV −1 for models with E/C = 8/3 like the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) model [28], which is still in the allowed region by the recent axion dark matter search, with m a around 6.7 µeV [84].
The Primakoff effect plays a key role in axion searches. For example, the working principle for an axion helioscope [85,86] is that axions produced in the core of the Sun are converted back into photons in a strong magnetic field. Clearly, if the ratio E/C = 2, which is   (2) CHPT [53]. For the axion-photon coupling g aγγ , only the model-independent part, denoted by g QCD aγγ , is shown.
quite a possibility as shown by Kaplan in ref. [71], then the g aγγ would be highly suppressed. The axion detection using the Primakoff effect, such as microwave cavity experiments, or light shining through wall experiments (for a recent review, see ref. [87]) would thus be extremely difficult. Here, we present the reference values of g aγγ for E/C = 2 and 8/3: +0.71(4) × 10 −3 /f a , E/C = 8/3.

(4.11)
With the expressions of the axion mass and the axion-photon coupling, it is straightforward to estimate the axion lifetime, namely, (4.12) As the axion lifetime is inversely proportional to m 5 a , the axion is more stable when its mass is smaller. The axion lifetime is estimated as τ a→γγ 10 33 s if the lower limit f a 0.5 × 10 9 GeV is employed. Such a cosmologically stable particle is a well-motivated cold dark matter candidate [32,88].

Summary
In this paper, we have calculated the QCD θ-vacuum energy and in turn the axion potential up-to-and-including NLO corrections in SU(N ) CHPT. Unlike the SU(2) case, no analytic solutions exist for SU(N ) with N ≥ 3. We work out for the first time a recursion relation for φ f , up to an arbitrary order in θ. Then, as an extension of ref. [54], by expanding the oneloop effective generating functional around the free-field configuration in a θ-vacuum, we have calculated the θ-vacuum energy density up NLO, including the one-loop contribution, in SU(N ) CHPT with N non-degenerate quark flavors. With the recursion relation for the φ f angles, one can compute any-order cumulants of the QCD topological charge distribution as well as the axion mass and self-couplings.
Since the QCD axion potential as a function of a/f a takes the same form as the QCD θ-vacuum energy as a function of θ, we have also calculated the axion mass and selfcoupling to NLO from the SU(3) θ-vacuum energy density taking into account the strong isospin breaking effects. With the determination of the LECs from experimental data and JHEP05(2020)001 lattice simulations, we have further evaluated the numerical values for axion mass and self-coupling up to NLO, which are similar to those obtained in the SU(2) case in ref. [53].
We also computed the axion-photon coupling up to O(p 6 ). Numerically, it is given by g aγγ = αem 2πfa [E/C − 2.05 (3)], which implies that if E/C = 2, the axion-photon coupling would be extremely small. In this case the axion searches using g aγγ , such as light shining through a wall or microwave cavity experiments, would be very difficult. This might also have an important impact on the axion electrodynamics as well as the possible existence of boson stars, in which the axion-photon coupling plays a crucial role. We use the following expansions, Once we solve all the coefficients C f,2m+1 , we then get the general solution of φ f . Let Notice that for the expansion of sin φ f in powers of θ in eq. (A.3), the terms at O θ 2n+1 are closely related to the partition of 2n + 1 into odd parts (e.g., the partitions of 5 into odd parts include 5, 3 + 1 + 1 and 1 + 1 + 1 + 1 + 1, see the left side of eq. (A.8)) studied in number theory. One can go to higher orders and solve for C f,2n+1 in the same way. Finally, one gets the recursion relation for all the coefficients as C f,2n+1 = n t=1 (k 1 ,...,kt) where k j are nonnegative integers, K t ≡ t j=1 k j , (k 1 ,...,kt) means that the sum runs over all possibilities of k j satisfying k 1 + · · · + (2t − 1)k t = 2n + 1, and Kt k 1 ,...,kt = K t !/(k 1 ! · · · k t !) are the multinomial coefficients.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.