Axion-meson mixing in light of recent lattice $\eta$-$\eta'$ simulations and their two-photon couplings within $U(3)$ chiral theory

We study the mixing of the QCD/QCD-like axion and light-flavor mesons $\pi^0, \eta, \eta'$ within the framework of $U(3)$ chiral perturbation theory up to next-to-leading order in this work. The axion-meson mixing formulas are calculated order by order in the $U(3)$ $\delta$-expansion scheme, namely the joint expansions of the momentum, light-quark masses and $1/N_C$. We provide axion-meson mixing relations in terms of the $\pi^0$-$\eta$-$\eta'$ mixing parameters and their masses. The recent lattice simulations on the $\eta$-$\eta'$ systems turn out to be able to offer valuable inputs to constrain the unknown low-energy constants. The relation of the mass and decay constant of the axion is then further explored based on our updated calculations. The two-photon couplings of the light-flavor mesons, together with the axion, are also investigated in the $U(3)$ chiral theory up to next-to-leading order in the $\delta$-counting scheme.


Introduction
The hypothetical particle axion provides an elegant solution to the long standing strong CP problem [1][2][3]. Since its first proposal in late 1970s, tremendous efforts have been put in various fields of physics to search this intriguing hypothesized particle, ranging from particle physics, astronomy, cosmology to condensed matter and optical physics, etc [4].
In the Peccei-Quinn (PQ) picture, axion corresponds to the pseudo-Nambu-Goldstone boson (pNGB) resulting from the breaking of the global U PQ (1) symmetry. In the low-energy regime, the most characteristic interaction of the axion is the coupling with the topological gluon density, i.e., a/f a G µνG µν , where G µν is the gluon field strength tensor, its dual isG µν = ε µνρσ G ρσ /2 with ε µνρσ the Levi-Civita antisymmetric tensor and f a stands for the axion decay constant. The essential point of the PQ mechanism to solve the strong CP problem relies on the cancellation of the CP violating term θG µνG µν in QCD by the dynamical generation of a proper vacuum expectation value (VEV) from the axion field. Different ultraviolet (UV) models in the new physics sector can lead to rather different couplings of the axion with other standard-model (SM) particles, such as the leptons, quarks, photon and W/Z bosons [2,3,[5][6][7][8], although they can be constructed to give a universal GG coupling. Due to the nonperturbative nature of the gluons in the low-energy region, the common way to study the interactions of the axion is to first perform the axial transformation of the quark fields to eliminate the a/f a G µνG µν operator from the beginning and then match to the axion chiral perturbation theory (χPT), which encodes the axion field together with the pions in the SU (2) case, and the octet of π, K, η 8 in the SU (3) case. On the other hand, it is also possible to directly match the anomalous axion term a/f a G µνG µν with the χPT operators. In this regard, it is reminded that the large mass of the singlet pseudoscalar η 0 (the dominant component of the physical η ′ meson), compared to those of the pseudoscalar octet π, K and η 8 , can be mainly attributed to the anomalous breaking of the QCD U A (1) symmetry. Therefore, it implies that another way to introduce the axion field into χPT can be similar to the case of the singlet η 0 . The pioneer work to study the influence of the η ′ on the axion properties can trace back to Refs. [9,10]. Recently many works aiming at the extended descriptions of the axion-η ( ′ ) interactions are proposed from different points of view [11][12][13][14][15][16][17][18][19][20], including the axion production from the η/η ′ and kaon decays, CP violating axion signals, the axion-baryon couplings, the model-independent part of the axion-photon-photon coupling, etc. In this work, we will further pursue a systematical calculation of the axion-φ (φ = π, η, η ′ ) interactions order by order within the U (3) χPT by employing the δ expansion scheme [21][22][23], namely, a simultaneous expansion in powers of the momenta, light-quark masses and 1/N C , i.e. δ ∼ p 2 ∼ m q ∼ 1/N C . In addition, we also study the relations of the model-independent two-photon couplings of the axion, together with those of the π 0 , η and η ′ mesons, by performing the next-to-leading order (NLO) U (3) χPT calculations.
The layout of this paper is as follows. The relevant chiral Lagrangians up to NLO are elaborated in Sec. 2. We address the mixing formalism at leading order (LO) for the fourparticle system involving π 0 , η, η ′ and the axion in Sec. 3. The calculation of the NLO mixing, the fits to relevant lattice data and numerical analyses of the transitions matrix elements, and the masses of π 0 , η, η ′ and the axion are given in Sec. 4. The two-photon couplings of the aforementioned four particles are then studied up to NLO in the U (3) χPT in Sec. 5. A short summary and conclusions are given in Sec. 6.

Axion U (3) χPT up to next-to-leading order
The axion is characterized by the interaction with the gluons and its effective operators can be written as where G i µν andG i µν are the gluon field strength tensor and its dual, with i the color indices. The second term in Eq. (1) is considered to be model-independent, due to its relevance of solving the strong CP problem. For the bare mass m a,0 of the axion in the third term, its value is usually assumed to be vanishing in the minimal QCD axion setup [1][2][3], although it is also possible to have a nonvanishing m a,0 to solve the strong CP problem [4,[24][25][26][27]. The direct couplings of the axion with the photon and fermions heavily rely on the specific model-building considerations in the BSM sector. In this work, we focus on the axion interactions with the hadrons and photon that are purely induced by the effective Lagrangian in Eq. (1), i.e. the model-independent parts of the axion-hadron and axion-photon interactions.
The next step is to match the effective Lagrangian in Eq. (1) to the axion chiral perturbation theory [28,29], which provides a reliable framework to study the axion-hadron interactions order by order. In the low-energy QCD, apart from the chiral symmetry breaking, another distinct feature is the QCD U A (1) anomaly, i.e. the anomalous breaking of the U A (1) symmetry by topological charge density ω(x) = α s G µνG µν /(8π), which gives a natural explanation of the large mass of the singlet η 0 even in the chiral limit. In this work we will carry out the study sticking to the U (3) χPT by employing the large N C argument [30]. This implies that one can explicitly include the axion field into the χPT Lagrangian in a similar way as the situation of the η ′ . To be more specific, we use the δ-expansion scheme to arrange the various contributions in U (3) χPT by simultaneously considering the expansions in the momenta, light-quark masses and 1/N C [21][22][23].
To set up the notations, we briefly recapitulate the way to construct the LO operators of U (3) χPT. From the large N C point of view, the singlet η 0 would become the ninth pNGB in the large N C and chiral limits, since the instanton effect via the QCD U A (1) anomaly is 1/N C suppressed. Then, the dynamical fields in low energy QCD form the pNGB nonet, which is usually parameterized as U = exp i √ 2Φ/F , with the pNGB nonet matrix given by The QCD U A (1) anomaly effect can be incorporated in effective Lagrangian via the operator − τ 2 (−i log det U ) 2 , where τ corresponds to the topological susceptibility [23]. Taking into account that det U = exp(Tr log U ), the former operator can be cast as −3τ η 2 0 /F 2 , which is nothing but the mass term for singlet η 0 in the chiral limit. In practice it is convenient to introduce M 2 0 = 6τ F 2 as the LO mass squared for η 0 . In the large N C counting rule, the topological susceptibility τ and the pNGB decay constant F scale like O(1) and O( √ N C ), respectively [23,31], which indicate that M 2 0 behaves as O(1/N C ). Under the PQ assumption, the CP violating θGG term is canceled by the VEV part of the axion via the aGG term in Eq. (1). Therefore, this indicates that the anomalous aGG effect can be included in the chiral effective Lagrangian in a similar fashion as the θ term, the latter of which is discussed in great detail in Ref. [23]. According to the recipe in the previous reference, the axion field can be introduced to the LO U (3) χPT together with the anomalous mass for For phenomenological convenience, we will always use the parameter M 0 to replace the topological susceptibility τ in the following discussions. On general grounds, the axion chiral transformation on the quark fields needed to remove the aGG term in Eq. (1) is of the same type as the singlet axial chiral transformation U A (1) that is parameterized as exp(iη 0 2/3). This observation drives to the necessity of adding together the fields of the axion and η 0 in the large N C and chiral limit, where the latter is a pNGB of QCD parameterized as coordinates in the coset space [32]. Under these circumstances, the LO U (3) χPT including the axion field takes the form where the axion field is introduced via and other basic χPT building blocks are given by Here, s, p, v µ , and a µ are external scalar, pseudoscalar, vector and axial-vector external sources, respectively, introduced as spurion fields. The quark-mass corrections are introduced by fixing s = M q ≡ diag(m u , m d , m s ). The parameter M 0 in the last term of Eq. (3) corresponds to the leading-order mass of the η 0 , which squared is proportional to the topological susceptibility [23].
The leading scaling behavior of M 2 0 is O(1/N C ) in the δ-counting scheme [30]. Another approach frequently used in literature to introduce the axion field in χPT, e.g. as those in Refs. [28,29], is to first perform the axial transformation of the quark fields to eliminate the axion-gluon operator in Eq. (1), i.e., by taking where Q a is a 3×3 matrix spanned in the light three-flavor space. Such a transformation induces two additional terms in the effective Lagrangian (1): − aαs 8πfa G µνG µν Tr(Q a ) and − ∂µa 2faq γ µ γ 5 Q a q. The former term exactly cancels the axion-gluon operator in Eq. (1), after taking the constraint Tr(Q a ) = 1. While, the latter term describes the extra axion-quark interaction that is induced by the axial transformation (6). In addition, it is easy to demonstrate that this transformation will also introduce the axion field into the quark masses: Qa , which can lead to non-derivative axion interactions. In the SU (3) χPT, the LO mass mixing between the axion and the neutral unflavored π 0 and η 8 can be avoided by taking Q a ∝ M −1 q . However, in the U (3) case, even if one imposes the latter form for Q a , the mixing between the axion and the singlet η 0 still exists at LO. Therefore, in this work we will not perform the transformation of quark fields (6), and use instead the original U (3) Lagrangian in Eq. (3) to proceed the calculations. The physical results should be the same regardless of keeping the aGG term or replacing it via the axial transformation (6), although the intermediate steps in the χPT calculations can look different.
In the rest of the discussions, we will stick to the method by including the axion in the U (3) χPT through the X field of Eq. (4). In this case, the axion U (3) χPT Lagrangian coincides with the standard one with the obvious replacement of the X field. When restricting to the axion-meson mixing, the relevant NLO Lagrangian under the δ-counting rule consists of four operators where the first two terms accompanied by L 5 and L 8 share the same forms as those from the conventional SU (3) χPT [31]. Within the framework of χPT, the N C order of a given operator can be inferred by the number of traces in the flavor space [31]. Generally speaking, one extra trace will bring one additional 1/N C suppression order to the effective operator. By taking into account the identity of log det U = log U , it is straightforward to conclude that the first two terms in Eq. (7) are counted as O(p 4 , N C ), and the remaining two terms with Λ 1 and Λ 2 that only appear in the U (3) case [23] are counted as O(p 2 , N 0 C ). The two-photon decays of the light pseudoscalar mesons are governed by the Wess-Zumino-Witten Lagrangian and the relevant LO Lagrangian [33,34] is which is counted as O(p 4 , N C ) in the δ-counting scheme. The quantity Q stands for the matrix of the electric charges of the light quarks, i.e. Q = Diag( 2e 3 , − e 3 , − e 3 ), with e the magnitude of the electron charge, and A µ stands for the photon field. At NLO there are two additional terms [35,36] with f µν where denote the field strength tensors of left-hand and right-hand external sources, respectively. One can take F µ L =F µ R = −eQA µ to obtain the interaction vertexes with photons. The t 1 and k 3 terms in Eq. (9) belong to the O(p 6 , N C ) and O(p 4 , N 0 C ) types of operators in the δ counting, respectively. When restricting to the pseudoscalar-photon-photon case, the Lagrangian (9) reduces to According to the Lagrangians in Eqs. (3) and (8), the LO interactions of the axion with the pNGBs and photons are purely caused by the mixing of the axion and the neutral unflavored pNGBs. Therefore it is mandatory to first address the axion-pNGBs mixing problem.

Revisiting the mixing formalism at LO
In principle, there are both kinetic and mass mixings between different states [37,38]. At LO in the δ counting, the π 0 , η, η ′ and a will get mixing purely via the mass terms from the Lagrangian (3), while the kinetic mixing only starts to appear in the NLO Lagrangian (7). We calculate the mixing of the axion and π 0 , η, η ′ order by order below.
In Refs. [39][40][41][42], the η-η ′ mixing has been studied up to next-to-next-to leading order in the δ expansion in the isospin limit. As first pointed out in Ref. [39], it is convenient to use the η andη ′ that are diagonalized at LO, if one attempts to perform a systematical higher-order calculation in the U (3) chiral theory relying on the δ counting. The key reason is that the LO mixing vertex of the η 8 and η 0 is formally counted as the same order of their LO masses, implying that an infinity insertion of the LO mixing vertex of η 0 and η 8 in the loop calculations will not get suppressed in the δ-counting scheme, see Fig. 1 of Ref. [39] for illustrations. This cumbersome problem can be avoided by first performing the LO diagonalization of η 8 and η 0 and using the LO diagonalized fieldsη andη ′ , whose relations are given by with c θ = cos θ and s θ = sin θ. The reason behind is that the mixing strengths betweenη and η ′ are suppressed at least by one higher order in the δ-counting rule. From the Lagrangian (3), the LO mixing angle θ and masses ofη andη ′ can be calculated where ∆ 2 = m 2 K − m 2 π , and m π and m K are the LO masses of the pion and kaon in the isospin symmetric limit, respectively.
New subtle complexities will arise when simultaneously including the π 0 , the axion a, η and η ′ . On one hand, the isospin breaking (IB) contributions can not be ignored any more, since the mixings between the π 0 and the remaining fields η, η ′ and a are caused by the IB effects. It is noted that we only consider the IB effects arising from the strong interaction, namely the mass difference between the up and down quarks. On the other hand, the mixings between the axion a and other fields π 0 , η and η ′ are suppressed by the factor F/f a . Additional power counting rules will be introduced to simplify the discussions, apart from the conventional δ-counting that is widely employed in the η-η ′ case. Generally speaking, the magnitudes of the IB corrections are around 1 ∼ 2%. Compared to the correction in the δ counting, which can be naively estimated to be around 30% [1/(N C = 3)], it is justified to just keep the leading nonvanishing IB terms. The magnitude of the axion-related suppression factor F/f a is definitely much more smaller than that of the IB corrections. Therefore, at the current stage, it is safe to just take the leading F/f a terms. More precisely, we aim at the systematical calculations of the higherorder corrections in the δ-counting scheme by keeping the leading IB and F/f a terms in this work.
Similarly as the η-η ′ case, it is much more convenient to use the fields π 0 , η, η ′ and a that already diagonalize LO mass term to systematically calculate the higher-order corrections. Notice that we have usedη andη ′ to denote the LO diagonalized η and η ′ fields in the isospin symmetric case, while η and η ′ stand for the fields after including the IB effects. As stated before, the leading terms in the IB and F/f a expansions will be kept in our calculation and, in this case, one can obtain the diagonalized states π 0 , η, η ′ and a by performing the following field redefinitions where π 0 and a correspond to the bare states entering in the Lagrangian (3). The matrix elements v ij are determined by the mass mixing terms from the LO Lagrangian (3), and their explicit expressions are found to be where m π , m η and m η ′ correspond to the LO masses of the π, η and η ′ mesons, and ǫ = B(m u − m d ) corresponds to the leading strong IB factor. Regarding the matrix element v 23 in Eq. (16), it describes the IB contribution to the η-η ′ mixing, which is expected to play negligible roles comparing with the SU (3) symmetry breaking effects. Therefore in this work we will neglect the IB contributions to the η-η ′ mixing, i.e., v 23 will be set to zero throughout. At the LO electromagnetic (EM) correction, the Dashen's theorem tells us that the EM corrections to the kaon masses equal to the mass differences of the pions [36,43] . After the diagonalization at LO, the mass of the axion field a is found to be which by taking m a,0 = 0 reduces to the minimal setup of the QCD axion case It should be stressed that the neat analytical expressions in Eqs. (17)-(23) are obtained by keeping the leading IB effects in π 0 -a and π 0 −η ( ′ ) mixing and neglecting the IB contributions to the a-η ( ′ ) mixing and the masses of the axion and pNGBs. The concise formulas in Eqs. (22) and (23) give us direct access to estimate the influences from the pertinent η-η ′ mixing parameters to the axion mass in different scenarios. We take the minimal QCD axion case, i.e. Eq. (23), to carry out some exploratory phenomenological discussions. By taking the LO expressions of the η-η ′ mixing in Eqs. (13), (14) and (15), it is easy to demonstrate that the axion mass squared can be recast as which reduces to the well celebrated result m 2 a = F 2 m 2 π /(4f 2 a ) [2] by keeping the leading expansions of m 2 π /m 2 K and m 2 π /M 2 0 .
For the LO fit in Ref. [44], the value of M 0 is determined to be M 0 = 820.0 MeV, which leads to By substituting these values to Eqs. (19)- (22) with m a,0 = 0, we obtain m a = 6.1 µeV where the mass agrees well with the recent determinations [29,45].
with M 0 = 820 MeV and θ = −19.6 • . In turn, if we take M 0 = 820 MeV and θ = −10.0 • the result turns out to be m a = 14.5 µeV The obvious changes of the predictions in Eqs. (26)-(28) by taking different phenomenological inputs for the mη and mη′, give hints that there could be potentially noticeable higher order corrections. This urges us to further continue the NLO calculations and verify their impacts on the axion properties.

Mixing at NLO and Phenomenological discussions
Up to the NLO in the δ counting, the bilinear terms involving the π 0 , η, η ′ and the axion field a, which are diagonalized at LO, can be generally written as where the δ X k with subscript k denotes the corrections from higher orders to the kinetic terms and other δ j corresponds to the higher order contributions to the mass terms. It is noted that the higher derivative terms are contributed by the O(p 6 ) and higher-order operators [42] and hence are absent at NLO. Comparing with the results of only focusing the η-η ′ mixing in Ref. [42,44], we extend the calculations by simultaneously including the π 0 and the axion a in the above equation. All the δ i parameters in Eq. (29) can be calculated order by order in the U (3) chiral perturbation theory. Their explicit expressions calculated from the NLO Lagrangians of Eq. (7) are given in the Appendix.
Aside from the mass mixing terms, one also has to deal with the kinetic mixing at NLO. This can be done in a two-step procedure [42]. In the first step, not only one needs to eliminate the kinetic mixing terms, but also one should make sure that all the fields are in the canonical normalization, i.e., the bilinear derivative term of each field is multiplied by 1/2. In the second step, the mass mixing will be handled by an orthogonal matrix that does not spoil the already accomplished diagonalization of the kinetic-energy term [37,38]. We will stick to the δ expansion up to NLO in all of these procedures. The diagonalized canonical fields, which are labeled with hats on top of each state, can be obtained via the following field redefinitions where the x ij and y ij are introduced to deal with the kinetic and mass mixing terms, respectively. The matrix elements of x ij are found to be and the expressions of the y ij read By keeping the contributions up to NLO, the masses of the diagonalized canonical states take the form and the kaon mass squared up to NLO is similarly given by Up to NLO the relations between the physical states (denoted by the hatted fields) and the LO diagonalized ones reduce to To combine the LO relation (16), the mixing matrix between the physical states and the bare ones are given by where the NLO corrections are collected in the z ij and their explicit forms are given by As a result of combining Eqs. (39) and (12), one can obtain Alternatively one could also use the quark-flavor bases of η q and η s , and they relate with the octet-singlet η 8 , η 0 bases via where the constituent quark contents of η q and η s are (ūu +dd)/ √ 2 andss, respectively. The LO diagonalizedη,η ′ states can be decomposed in terms of the quark-flavor base with φ qs = θ id + θ and the ideal mixing angle θ id = arcsin( √ 2/ √ 3). Therefore, by combining Eqs. (39) and (43), one can in turn obtain the relations between the physical states and those in the quark-flavor bases. The results share the same forms as those in Eq. (41) with the explicit replacement of the angle θ by the angle φ qs defined in Eq. (43).
The two-angle mixing formula proposed in Ref. [22] gives an intuitive relation between the physical states denoted byη ( ′ ) here and the ones in octet-singlet basis which naturally reduces to the conventional mixing relation with one angle by taking F 8 = F 0 = F and θ 0 = θ 8 = θ. Similarly one could also introduce the two-angle mixing formalism to the quark-flavor base Remarkable progress has been made in recent years for the lattice calculation of the η-η ′ system. Both their masses and mixing angles at different pion masses are obtained by many lattice collaborations. In this work, we fix the unknown low-energy constants (LECs) of χPT by performing fits to relevant lattice data, including the η/η ′ masses and their mixing angle at different m π , the weak decay constants of pion and kaon, and the pion-mass dependences of m K . Apart from the lattice data of the η/η ′ masses at unphysically large m π from the ETMC [46], UKQCD [47], RBC/UKQCD [48], HSC [49] that have been analyzed in Refs. [42,44], we also include the new results from the RQCD collaboration [50] in the present study. In order to make direct comparisons of other data, only the lattice simulations with physical strange quark mass from Ref. [50] are taken into account here. The various lattice data of the masses for η and η ′ are shown in the left panel of Fig. 1, together with their mixing angles in the quark-flavor base in the right panel. To be in accordance with the lattice setup, we take the average of θ q and θ s to fit the latter data set. It is noted that the phenomenological determinations also favor rather similar values for θ q and θ s [42,51,52]. The decay constants F q and F s defined in the mixing matrix of the quark-flavor base as functions of pion masses are provided by the lattice study [46] and we will take these kinds of data in our fits too, as shown in Fig. 2. For the pion-mass dependences of the F π , F K and m K , we take into account the lattice data up to m 2 π = 0.1 GeV 2 , as explicitly illustrated in Fig. 3.   The lattice data of the η/η ′ masses are taken from ETMC [46], UKQCD [47], RBC/UKQCD [48], HSC [49] and RQCD [50]. For the η/η ′ masses from Ref. [50], only the results from the ensemble with approximately physical mass of strange quark are included. The lattice data of mixing angles in the quarkflavor base are taken from Ref. [46].  [53,54]. The data in the right panel are from Ref. [55].

Parameters NLO Fit
The values of the LECs from the fits are summarized in Table 1. The resulting curves for the masses of η and η ′ and their mixing angles from the fits are given in Fig. 1. It is noted that the LO fit to the masses of η and η ′ with just the single parameter M 0 can reasonably reproduce the lattice data with M 0 = 820 MeV [42,44]. Therefore we will fix M 0 at this value during the NLO fits. It is verified that by releasing the M 0 the fits do not improve and the values of the parameters in Table 1 barely change. The resulting parameters from the revised NLO fits are close to the previous ones given in Refs. [42,44]. With the fitted parameters in Table 1, we are ready to predict the important quantities related with the axion.
Since the bare mass of the axion, i.e. m a,0 in Eq. (1), is explicitly kept throughout our calculations, e.g. Eqs. (36), (17)-(21) and (32), it is straightforward for us to explore the socalled axion-like scenarios by assigning some specific nonvanishing values to m a,0 . Nevertheless, in this work we will mainly focus on the phenomenological predictions for the QCD-axion scenario by taking m a,0 = 0 in Eq. (1). The explicit values of the transition matrix elements in Eq. (41) between the physical states ofπ 0 ,η,η ′ ,â and the bare states of π 0 , η 8 , η 0 , a are determined to be where the first and second entries in each matrix element correspond to the LO and NLO contributions, respectively. For the numbers in the last row and last column accompanying 1/f a , they are given in units of MeV, while the number accompanying 1/f 2 a is given in units of MeV 2 . The main focus of our work is the relative corrections from the NLO part compared to the LO one.
For the mixing strengths between theπ 0 and η 0,8 , and also the mixing between theη ( ′ ) and the π 0 , which are all proportional to the leading IB factor ǫ, the NLO parts in the δ counting gives rather large relative corrections compared to the LO results, ranging from around 60% up to 100%. For the η-η ′ mixing, it is also found that the NLO contribution can be as large as around 60% relative to the LO part. However, it is interesting to note that all the mixing strengths between the axion and the light-flavor pseudoscalar π 0 , η, η ′ , i.e. the numbers in the last column and the last row, the relative NLO corrections in the δ counting to the LO parts are small, ranging from around 2% up to around 15%. The uncertainties of the mixing strengths given in Eq. (47), estimated by the using the error bands of the LECs in Table 1, turn out to be mild as well.
The contributions from the LO and NLO parts to the mass squared of the light-flavor pseudoscalar mesons and the axion in Eqs.
where the first and second entries inside the square brackets denote the results from the LO and NLO terms, respectively. Comparing with the LO case, it is clear that the NLO correction brings the η mass much closer to its physical value. When compared with the LO case, the NLO correction brings the η mass much closer to its physical value, although it somewhat worsens the description of the η ′ mass. However, we should notice that while the NLO contribution to the η ′ mass is quite small, less than 3% of the LO value, it is a 12% for the η mass, which allows us to match its experimental value. The NLO contribution to the axion mass turns out to be rather small.

Two-photon couplings
Relying on the previous results of the mixing relations, we are ready to study the twophoton decays of the π 0 , η, η ′ and a, based on the LO and NLO WZW operators in Eqs. (8) and (11), respectively. By inserting the mixing relations (39) into the LO WZW Lagrangian (8) and neglecting all the IB terms, the two-photon couplings of the physical states read where all of the NLO terms are introduced through the mixing. Since we stick to the NLO calculation, only the LO mixing relations (16) are needed to insert into the NLO WZW Lagrangian (11) and the resulting two-photon interactions are given by The two-photon decay amplitude of φ → γ(k 1 )γ(k 2 ) with φ = π 0 , η, η ′ and a can be written as where the two-photon coupling strengths F φγγ read We use the decay widths of π 0 → γγ, η → γγ and η ′ → γγ from the most recent PDG average [56] to estimate their two-photon couplings F Exp F Exp Those couplings can be then exploited to determine the NLO LECs from the WZW Lagrangian (9) entering in Eqs. (52)- (54). Their explicit values turn out to be leading to which perfectly reproduce the experimental inputs.
With the values of the NLO LECs of L NLO W ZW in Eq. (59) and the fitted parameters in Table 1, we can now give our prediction to the two-photon coupling of the axion F aγγ up to NLO in the δ counting where the first entry in the numerator on the right hand side corresponds to the LO contribution and the second one denotes the NLO contribution. The two-photon coupling F aγγ is related with the g aγγ used in Refs. [29,45] via where the numbers inside the bracket can be compared to the results of 1.92 ± 0.04 [29] and 2.05 ± 0.03 [45] from the SU (2) and SU (3) χPT analyses up to NLO, respectively. The determination of the magnitude of g aγγ from the NLO U (3) calculation looks a bit smaller than those from the NLO studies in the SU (2) and SU (3) cases. It is noted that the chiral loops start to appear at NLO in the SU (2) and SU (3) χPT in the conventional chiral power counting. While, in the δ-counting scheme, the chiral loops only enter at NNLO in the U (3) χPT. To consistently take into account all the pieces of the two-photon coupling of the axion at NNLO in the U (3) χPT, one also needs to extend the axion-meson mixing calculations at the same order and we leave this task to a future work.

Summary and Conclusions
In this work, U (3) chiral perturbation theory is demonstrated to be able to provide a useful framework to calculate the matrix elements of the mixing between the π 0 , η, η ′ and the axion order by order in the joint expansions of the momenta, light-quark masses and 1/N C , i.e. the δ-counting scheme. We have performed the complete calculation of the axion-meson mixing and the two-photon couplings of the axion and light-flavor pseudoscalar mesons up to NLO in the δ-counting rule within the framework of the U (3) chiral perturbation theory. The unknown chiral low-energy constants are fixed through the fits to a large amount of lattice QCD data, consisting of the pion-mass dependences of η-η ′ mixing data, kaon mass, and the decay constants of the pion and kaon. Reasonable reproductions of the various lattice data are achieved with the U (3) contributions up to NLO. The mixing matrix elements of the π 0 , η, η ′ and the axion are then further used to calculate their two-photon couplings, together with the NLO Wess-Zumino-Witten Lagrangian, where the parameters in the latter Lagrangian are fixed by the experimental two-photon couplings of the π 0 , η and η ′ .
The determined chiral low-energy constants are then used to predict the QCD-axion masses, the mixing strengths of the axion and the π 0 , η, η ′ , and the two-photon coupling of the QCD axion. The NLO contributions in the δ counting to the various axion quantities relative to the LO ones are found to be small. This work paves the way to systematically calculate the interactions between the axion and the light-flavor pseudoscalar mesons π, K, η and η ′ order by order in the joint expansions of momenta, light-quark masses and 1/N C .
Appendix: Next-to-leading order coefficients for bilinear terms Substituting the LO mixing relations of Eq. (16) into the NLO Lagrangian (7), one can calculate all the NLO pieces of the bilinear terms in Eq. (29) and their explicit expressions read Although the NLO formulas in Eqs. (77)-(86) can be found in Ref. [42], we show their explicit expressions here for the sake of completeness.