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

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_{\rm 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^{\mu\nu}\tilde{G}_{\mu\nu}$, where $G^{\mu\nu}$ is the gluon field strength tensor, its dual is $\tilde{G}_{\mu\nu}=\varepsilon_{\mu\nu\rho\sigma}G^{\rho\sigma}/2$ with $\varepsilon_{\mu\nu\rho\sigma}$ 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 $\theta G^{\mu\nu}\tilde{G}_{\mu\nu}$ 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 $G\tilde{G}$ 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^{\mu\nu}\tilde{G}_{\mu\nu}$ operator from the beginning and then match to the axion chiral perturbation theory ($\chi$PT), which encodes the axion field together with the pions in the $SU(2)$ case, and the octet of $\pi,K,\eta_{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^{\mu\nu}\tilde{G}_{\mu\nu}$ with the operators. In this regard, it is reminded that the large mass of the singlet pseudoscalar $\eta_{0}$ (the dominant component of the physical $\eta^{\prime}$ meson), compared to those of the pseudoscalar octet $\pi,K$ and $\eta_{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 $\chi$PT can be similar to the case of the singlet $\eta_{0}$. The pioneer work to study the influence of the $\eta^{\prime}$ on the axion properties can trace back to Refs. [9, 10]. Recently many works aiming at the extended descriptions of the axion-$\eta^{(^{\prime})}$ interactions are proposed from different points of view [11, 12, 13, 14, 15, 16, 17, 18, 19, 20], including the axion production from the $\eta/\eta^{\prime}$ 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-$\phi$ ($\phi=\pi,\eta,\eta^{\prime}$) interactions order by order within the $U(3)$ by employing the $\delta$ expansion scheme [21, 22, 23], namely, a simultaneous expansion in powers of the momenta, light-quark masses and $1/N_{C}$, i.e. $\delta\sim p^{2}\sim m_{q}\sim 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 $\pi^{0},\eta$ and $\eta^{\prime}$ mesons, by performing the next-to-leading order (NLO) $U(3)$ 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 four-particle system involving $\pi^{0},\eta,\eta^{\prime}$ 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 $\pi^{0},\eta,\eta^{\prime}$ 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)$ $\chi$PT in Sec. 5. A short summary and conclusions are given in Sec. 6.

The axion is characterized by the interaction with the gluons and its effective operators can be written as

$\displaystyle\mathcal{L}_{a}^{G}=\frac{1}{2}\partial_{\mu}a\partial^{\mu}a+\frac{a}{f_{a}}\frac{\alpha_{s}}{8\pi}\sum_{i=1}^{8}G^{i}_{\mu\nu}\tilde{G}^{i,\mu\nu}-\frac{1}{2}m_{a,0}^{2}\,a^{2}\,,$

(1)

where $G^{i}_{\mu\nu}$ and $\tilde{G}^{i}_{\mu\nu}$ 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 [24, 25, 26, 27, 4]. 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 $\omega(x)=\alpha_{s}G_{\mu\nu}\tilde{G}^{\mu\nu}/(8\pi)$, which gives a natural explanation of the large mass of the singlet $\eta_{0}$ even in the chiral limit. In this work we will carry out the study sticking to the $U(3)$ by employing the large $N_{C}$ argument [30]. This implies that one can explicitly include the axion field into the Lagrangian in a similar way as the situation of the $\eta^{\prime}$. To be more specific, we use the $\delta$-expansion scheme to arrange the various contributions in $U(3)$ by simultaneously considering the expansions in the momenta, light-quark masses and $1/N_{C}$ [22, 23, 21].

To set up the notations, we briefly recapitulate the way to construct the LO operators of $U(3)$ . From the large $N_{C}$ point of view, the singlet $\eta_{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}\big{(}i\sqrt{2}\Phi/F\big{)}$, with the pNGB nonet matrix given by

$$\Phi\,=\,\left(\begin{array}[]{ccc}\frac{1}{\sqrt{2}}\pi^{0}+\frac{1}{\sqrt{6}}\eta_{8}+\frac{1}{\sqrt{3}}\eta_{0}&\pi^{+}&K^{+}\\ \pi^{-}&\frac{-1}{\sqrt{2}}\pi^{0}+\frac{1}{\sqrt{6}}\eta_{8}+\frac{1}{\sqrt{3}}\eta_{0}&K^{0}\\ K^{-}&\overline{K}^{0}&\frac{-2}{\sqrt{6}}\eta_{8}+\frac{1}{\sqrt{3}}\eta_{0}\end{array}\right)\,.$$

(2)

The QCD $U_{A}(1)$ anomaly effect can be incorporated in effective Lagrangian via the operator $-\frac{\tau}{2}(-i\log\det U)^{2}$, where $\tau$ corresponds to the topological susceptibility [23]. Taking into account that $\det U=\exp({\rm Tr}\log U)$, the former operator can be cast as $-3\tau\eta_{0}^{2}/{F^{2}}$, which is nothing but the mass term for singlet $\eta_{0}$ in the chiral limit. In practice it is convenient to introduce $M_{0}^{2}=\frac{6\tau}{F^{2}}$ as the LO mass squared for $\eta_{0}$. In the large $N_{C}$ counting rule, the topological susceptibility $\tau$ and the pNGB decay constant $F$ scale like $O(1)$ and $O(\sqrt{N_{C}})$, respectively [23, 31], which indicate that $M_{0}^{2}$ behaves as $O(1/N_{C})$. Under the PQ assumption, the CP violating $\theta G\tilde{G}$ term is canceled by the VEV part of the axion via the $aG\tilde{G}$ term in Eq. (1). Therefore, this indicates that the anomalous $aG\tilde{G}$ effect can be included in the chiral effective Lagrangian in a similar fashion as the $\theta$ 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)$ together with the anomalous mass for $\eta_{0}$ via $-\frac{\tau}{2}(-i\log\det U+a/f_{a})^{2}=-\frac{\tau}{2}(\sqrt{6}\eta_{0}/F+a/f_{a})^{2}$. For phenomenological convenience, we will always use the parameter $M_{0}$ to replace the topological susceptibility $\tau$ in the following discussions. On general grounds, the axion chiral transformation on the quark fields needed to remove the $aG\tilde{G}$ term in Eq. (1) is of the same type as the singlet axial chiral transformation $U_{A}(1)$ that is parameterized as $\exp(i\eta_{0}\sqrt{2/3})$. This observation drives to the necessity of adding together the fields of the axion and $\eta_{0}$ in the large $N_{C}$ and chiral limit, where the latter is a pNGB of QCD parameterized as coordinates in the coset space $U_{L}(3)\otimes U_{R}(3)/U_{V}(3)$ [32]. Under these circumstances, the LO $U(3)$ including the axion field takes the form

$\displaystyle\mathcal{L}^{\rm LO}=\frac{F^{2}}{4}\langle u_{\mu}u^{\mu}\rangle+\frac{F^{2}}{4}\langle\chi_{+}\rangle+\frac{F^{2}}{12}M_{0}^{2}X^{2}\,,$

(3)

where the axion field is introduced via

$\displaystyle X=\log{(\det U)}+i\frac{a}{f_{a}}\,,$

(4)

and other basic building blocks are given by

	$\displaystyle U=u^{2}=e^{i\frac{\sqrt{2}\Phi}{F}}\,,\qquad\chi=2B(s+ip)\,,\qquad\chi_{\pm}=u^{\dagger}\chi u^{\dagger}\pm u\chi^{\dagger}u\,,$
	$\displaystyle u_{\mu}=iu^{\dagger}D_{\mu}Uu^{\dagger}\,,\qquad D_{\mu}U\,=\,\partial_{\mu}U-i(v_{\mu}+a_{\mu})U\,+iU(v_{\mu}-a_{\mu})\,.$		(5)

Here, $s$, $p$, $v_{\mu}$, and $a_{\mu}$ 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}\equiv{\rm 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 $\eta_{0}$, which squared is proportional to the topological susceptibility [23]. The leading scaling behavior of $M_{0}^{2}$ is $O(1/N_{C})$ in the $\delta$-counting scheme [30].

Another approach frequently used in literature to introduce the axion field in , 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

$\displaystyle q\to e^{i\frac{a}{2f_{a}}\gamma_{5}Q_{a}}q\,,\qquad[q=(u,d,s)^{\rm T}]\,,$

(6)

where $Q_{a}$ is a $3\times 3$ matrix spanned in the light three-flavor space. Such a transformation induces two additional terms in the effective Lagrangian (1): $-\frac{a\alpha_{s}}{8\pi f_{a}}G_{\mu\nu}\tilde{G}^{\mu\nu}\,{\rm Tr}(Q_{a})$ and $-\frac{\partial_{\mu}a}{2f_{a}}\bar{q}\gamma^{\mu}\gamma_{5}Q_{a}q$. The former term exactly cancels the axion-gluon operator in Eq. (1), after taking the constraint ${\rm 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: $M_{q}\to M_{q}(a)=e^{i\frac{a}{2f_{a}}Q_{a}}M_{q}e^{i\frac{a}{2f_{a}}Q_{a}}$, which can lead to non-derivative axion interactions. In the $SU(3)$ $\chi$PT, the LO mass mixing between the axion and the neutral unflavored $\pi^{0}$ and $\eta_{8}$ can be avoided by taking $Q_{a}\propto M_{q}^{-1}$. However, in the $U(3)$ case, even if one imposes the latter form for $Q_{a}$, the mixing between the axion and the singlet $\eta_{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 $aG\tilde{G}$ term or replacing it via the axial transformation (6), although the intermediate steps in the calculations can look different.

In the rest of the discussions, we will stick to the method by including the axion in the $U(3)$ through the $X$ field of Eq. (4). In this case, the axion $U(3)$ 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 $\delta$-counting rule consists of four operators

$\displaystyle\mathcal{L}^{\rm NLO}=L_{5}\langle u^{\mu}u_{\mu}\chi_{+}\rangle+\frac{L_{8}}{2}\langle\chi_{+}\chi_{+}+\chi_{-}\chi_{-}\rangle-\frac{F^{2}\,\Lambda_{1}}{12}D^{\mu}XD_{\mu}X-\frac{F^{2}\,\Lambda_{2}}{12}X\langle\chi_{-}\rangle\,,$

(7)

where the first two terms accompanied by $L_{5}$ and $L_{8}$ share the same forms as those from the conventional $SU(3)$  [31]. Within the framework of , 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=\langle\log U\rangle$, 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 $\Lambda_{1}$ and $\Lambda_{2}$ that only appear in the $U(3)$ case [23] are counted as $O(p^{2},N_{C}^{0})$.

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

$\displaystyle\mathcal{L}_{WZW}^{\rm LO}=-\frac{3\sqrt{2}e^{2}}{8\pi^{2}F}\varepsilon_{\mu\nu\rho\sigma}\partial^{\mu}A^{\nu}\partial^{\rho}A^{\sigma}\langle Q^{2}\Phi\rangle\,,$

(8)

which is counted as $O(p^{4},N_{C})$ in the $\delta$-counting scheme. The quantity $Q$ stands for the matrix of the electric charges of the light quarks, i.e. $Q={\rm Diag}(\frac{2e}{3},-\frac{e}{3},-\frac{e}{3})$, with $e$ the magnitude of the electron charge, and $A_{\mu}$ stands for the photon field. At NLO there are two additional terms [35, 36]

$\displaystyle\mathcal{L}_{WZW}^{\rm NLO}=it_{1}\varepsilon_{\mu\nu\rho\sigma}\langle f_{+}^{\mu\nu}f_{+}^{\rho\sigma}\chi_{-}\rangle+k_{3}\varepsilon_{\mu\nu\rho\sigma}\langle f_{+}^{\mu\nu}f_{+}^{\rho\sigma}\rangle X\,,$

(9)

with

$\displaystyle f_{+}^{\mu\nu}=uF_{L}^{\mu\nu}u^{\dagger}+u^{\dagger}F_{R}^{\mu\nu}u\,,$

(10)

where $F_{L}^{\mu\nu}=\partial^{\mu}F_{L}^{\nu}-\partial^{\mu}F_{L}^{\nu}-i[F_{L}^{\mu},F_{L}^{\nu}]$ and $F_{R}^{\mu\nu}=\partial^{\mu}F_{R}^{\nu}-\partial^{\mu}F_{R}^{\nu}-i[F_{R}^{\mu},F_{R}^{\nu}]$, with $F_{L}^{\mu}=v^{\mu}-a^{\mu}$ and $F_{R}^{\mu}=v^{\mu}+a^{\mu}$, denote the field strength tensors of left-hand and right-hand external sources, respectively. One can take $F_{L}^{\mu}$=$F_{R}^{\mu}$= $-eQA^{\mu}$ 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_{C}^{0})$ types of operators in the $\delta$ counting, respectively. When restricting to the pseudoscalar-photon-photon case, the Lagrangian (9) reduces to

	$\displaystyle\mathcal{L}_{WZW}^{\rm NLO}=$	$\displaystyle t_{1}\frac{32\sqrt{2}B}{F}\varepsilon_{\mu\nu\rho\sigma}\partial^{\mu}A^{\nu}\partial^{\rho}A^{\sigma}\langle\big{(}M_{q}\Phi+\Phi M_{q}\big{)}Q^{2}\rangle$
		$\displaystyle+16k_{3}\varepsilon_{\mu\nu\rho\sigma}\partial^{\mu}A^{\nu}\partial^{\rho}A^{\sigma}\langle Q^{2}\rangle\bigg{(}\frac{\sqrt{2}}{F}\langle\Phi\rangle+\frac{a}{f_{a}}\bigg{)}\,.$		(11)

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.

In principle, there are both kinetic and mass mixings between different states [37, 38]. At LO in the $\delta$ counting, the $\pi^{0},\eta,\eta^{\prime}$ 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 $\pi^{0},\eta,\eta^{\prime}$ order by order below.

In Refs. [39, 40, 41, 42], the $\eta$-$\eta^{\prime}$ mixing has been studied up to next-to-next-to leading order in the $\delta$ expansion in the isospin limit. As first pointed out in Ref. [39], it is convenient to use the $\mathring{\overline{\eta}}$ and $\mathring{\overline{\eta}}^{\prime}$ that are diagonalized at LO, if one attempts to perform a systematical higher-order calculation in the $U(3)$ chiral theory relying on the $\delta$ counting. The key reason is that the LO mixing vertex of the $\eta_{8}$ and $\eta_{0}$ is formally counted as the same order of their LO masses, implying that an infinity insertion of the LO mixing vertex of $\eta_{0}$ and $\eta_{8}$ in the loop calculations will not get suppressed in the $\delta$-counting scheme, see Fig. 1 of Ref. [39] for illustrations. This cumbersome problem can be avoided by first performing the LO diagonalization of $\eta_{8}$ and $\eta_{0}$ and using the LO diagonalized fields $\mathring{\overline{\eta}}$ and $\mathring{\overline{\eta}}^{\prime}$, whose relations are given by

$\displaystyle\left(\begin{array}[]{c}\mathring{\overline{\eta}}\\ \mathring{\overline{\eta}}^{\prime}\end{array}\right)=\left(\begin{array}[]{cc}c_{\theta}&-s_{\theta}\\ s_{\theta}&c_{\theta}\end{array}\right)\left(\begin{array}[]{c}\eta_{8}\\ \eta_{0}\end{array}\right)\,,$

(18)

with $c_{\theta}=\cos{\theta}$ and $s_{\theta}=\sin{\theta}$. The reason behind is that the mixing strengths between $\mathring{\overline{\eta}}$ and $\mathring{\overline{\eta}}^{\prime}$ are suppressed at least by one higher order in the $\delta$-counting rule. From the Lagrangian (3), the LO mixing angle $\theta$ and masses of $\mathring{\overline{\eta}}$ and $\mathring{\overline{\eta}}^{\prime}$ can be calculated

	$\displaystyle m_{\overline{\eta}}^{2}$	$\displaystyle=$	$\displaystyle\frac{M_{0}^{2}}{2}+m_{\overline{K}}^{2}-\frac{\sqrt{M_{0}^{4}-\frac{4M_{0}^{2}\Delta^{2}}{3}+4\Delta^{4}}}{2}\,,$		(19)
	$\displaystyle m_{\overline{\eta}^{\prime}}^{2}$	$\displaystyle=$	$\displaystyle\frac{M_{0}^{2}}{2}+m_{\overline{K}}^{2}+\frac{\sqrt{M_{0}^{4}-\frac{4M_{0}^{2}\Delta^{2}}{3}+4\Delta^{4}}}{2}\,,$		(20)
	$\displaystyle\sin{\theta}$	$\displaystyle=$	$\displaystyle-\left(\sqrt{1+\frac{\big{(}3M_{0}^{2}-2\Delta^{2}+\sqrt{9M_{0}^{4}-12M_{0}^{2}\Delta^{2}+36\Delta^{4}}\big{)}^{2}}{32\Delta^{4}}}~{}\right)^{-1}\,,$		(21)

where $\Delta^{2}=m_{\overline{K}}^{2}-m_{\overline{\pi}}^{2}$, and $m_{\overline{\pi}}$ and $m_{\overline{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 $\pi^{0}$, the axion $a$, $\eta$ and $\eta^{\prime}$. On one hand, the isospin breaking (IB) contributions can not be ignored any more, since the mixings between the $\pi^{0}$ and the remaining fields $\eta$, $\eta^{\prime}$ 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 $\pi^{0}$, $\eta$ and $\eta^{\prime}$ are suppressed by the factor $F/f_{a}$. Additional power counting rules will be introduced to simplify the discussions, apart from the conventional $\delta$-counting that is widely employed in the $\eta$-$\eta^{\prime}$ case. Generally speaking, the magnitudes of the IB corrections are around $1\sim 2\%$. Compared to the correction in the $\delta$ 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 higher-order corrections in the $\delta$-counting scheme by keeping the leading IB and $F/f_{a}$ terms in this work.

Similarly as the $\eta$-$\eta^{\prime}$ case, it is much more convenient to use the fields $\overline{\pi}^{0}$, $\overline{\eta}$, $\overline{\eta}^{\prime}$ and $\overline{a}$ that already diagonalize LO mass term to systematically calculate the higher-order corrections. Notice that we have used $\mathring{\overline{\eta}}$ and $\mathring{\overline{\eta}}^{\prime}$ to denote the LO diagonalized $\eta$ and $\eta^{\prime}$ fields in the isospin symmetric case, while $\overline{\eta}$ and $\overline{\eta}^{\prime}$ 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 $\overline{\pi}^{0}$, $\overline{\eta}$, $\overline{\eta}^{\prime}$ and $\overline{a}$ by performing the following field redefinitions

$$\left(\begin{array}[]{c}\overline{\pi}^{0}\\ \overline{\eta}\\ \overline{\eta}^{\prime}\\ \overline{a}\end{array}\right)\,=\,\left(\begin{array}[]{cccc}1+O(v^{2})&-v_{12}&-v_{13}&-v_{14}\\ v_{12}&1+O(v^{2})&-v_{23}&-v_{24}\\ v_{13}&v_{23}&1+O(v^{2})&-v_{34}\\ v_{14}&v_{24}&v_{34}&1+O(v^{2})\end{array}\right)\,\left(\begin{array}[]{c}\pi^{0}\\ \mathring{\overline{\eta}}\\ \mathring{\overline{\eta}}^{\prime}\\ a\end{array}\right)\,,$$

(22)

where $\pi^{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

	$\displaystyle v_{12}=$	$\displaystyle-\frac{\epsilon}{\sqrt{3}}\frac{c_{\theta}-\sqrt{2}s_{\theta}}{m_{\overline{\pi}}^{2}-m_{\overline{\eta}}^{2}}\,,$			(23)
	$\displaystyle v_{13}=$	$\displaystyle-\frac{\epsilon}{\sqrt{3}}\frac{\sqrt{2}c_{\theta}+s_{\theta}}{m_{\overline{\pi}}^{2}-m_{\overline{\eta}^{\prime}}^{2}}\,,$			(24)
	$\displaystyle v_{14}=$	$\displaystyle\frac{M_{0}^{2}\epsilon}{3\sqrt{2}(m_{a,0}^{2}-m_{\overline{\pi}}^{2})}\frac{F}{f_{a}}\bigg{[}\frac{s_{\theta}(c_{\theta}-\sqrt{2}s_{\theta})(m_{a,0}^{2}+m_{\overline{\pi}}^{2}-m_{\overline{\eta}}^{2})}{(m_{a,0}^{2}-m_{\overline{\eta}}^{2})(m_{\overline{\eta}}^{2}-m_{\overline{\pi}}^{2})}-\frac{c_{\theta}(\sqrt{2}c_{\theta}+s_{\theta})(m_{a,0}^{2}+m_{\overline{\pi}}^{2}-m_{\overline{\eta}^{\prime}}^{2})}{(m_{a,0}^{2}-m_{\overline{\eta}^{\prime}}^{2})(m_{\overline{\eta}^{\prime}}^{2}-m_{\overline{\pi}}^{2})}\bigg{]}\,,$
	$\displaystyle v_{24}=$	$\displaystyle-\frac{M_{0}^{2}s_{\theta}}{\sqrt{6}(m_{a,0}^{2}-m_{\overline{\eta}}^{2})}\frac{F}{f_{a}}\,,$			(26)
	$\displaystyle v_{34}=$	$\displaystyle\frac{M_{0}^{2}c_{\theta}}{\sqrt{6}(m_{a,0}^{2}-m_{\overline{\eta}^{\prime}}^{2})}\frac{F}{f_{a}}\,,$			(27)

where $m_{\overline{\pi}}$, $m_{\overline{\eta}}$ and $m_{\overline{\eta}^{\prime}}$ correspond to the LO masses of the $\pi$, $\eta$ and $\eta^{\prime}$ mesons, and $\epsilon=B(m_{u}-m_{d})$ corresponds to the leading strong IB factor. Regarding the matrix element $v_{23}$ in Eq. (22), it describes the IB contribution to the $\eta$-$\eta^{\prime}$ 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 $\eta$-$\eta^{\prime}$ 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 [43, 36], i.e. $(m_{K^{+}}^{2}-m_{K^{0}}^{2})_{\rm EM}=m_{\pi^{+}}^{2}-m_{\pi^{0}}^{2}$, and at this order one can write $\epsilon\equiv B(m_{u}-m_{d})=m_{K^{+}}^{2}-m_{K^{0}}^{2}-(m_{\pi^{+}}^{2}-m_{\pi^{0}}^{2})$.

After the diagonalization at LO, the mass of the axion field $\overline{a}$ is found to be

$\displaystyle m_{\overline{a}}^{2}=m_{a,0}^{2}+\frac{M_{0}^{2}F^{2}}{6f_{a}^{2}}\bigg{[}1+\frac{c_{\theta}^{2}M_{0}^{2}(2m_{a,0}^{2}-m_{\overline{\eta}^{\prime}}^{2})}{(m_{a,0}^{2}-m_{\overline{\eta}^{\prime}}^{2})^{2}}+\frac{s_{\theta}^{2}M_{0}^{2}(2m_{a,0}^{2}-m_{\overline{\eta}}^{2})}{(m_{a,0}^{2}-m_{\overline{\eta}}^{2})^{2}}\bigg{]}+O(\epsilon)\,,$

(28)

which by taking $m_{a,0}=0$ reduces to the minimal setup of the QCD axion case

$\displaystyle m_{\overline{a}}^{2}=\frac{M_{0}^{2}F^{2}}{6f_{a}^{2}}\bigg{[}1-\frac{c_{\theta}^{2}M_{0}^{2}}{m_{\overline{\eta}^{\prime}}^{2}}-\frac{s_{\theta}^{2}M_{0}^{2}}{m_{\overline{\eta}}^{2}}\bigg{]}+O(\epsilon)\,.$

(29)

It should be stressed that the neat analytical expressions in Eqs. (23)-(29) are obtained by keeping the leading IB effects in $\pi^{0}$-$a$ and $\pi^{0}-\eta^{(^{\prime})}$ mixing and neglecting the IB contributions to the $a$-$\eta^{(^{\prime})}$ mixing and the masses of the axion and pNGBs. The concise formulas in Eqs. (28) and (29) give us direct access to estimate the influences from the pertinent $\eta$-$\eta^{\prime}$ mixing parameters to the axion mass in different scenarios. We take the minimal QCD axion case, i.e. Eq. (29), to carry out some exploratory phenomenological discussions.

By taking the LO expressions of the $\eta$-$\eta^{\prime}$ mixing in Eqs. (19), (20) and (21), it is easy to demonstrate that the axion mass squared can be recast as

$\displaystyle m_{\overline{a}}^{2}=\frac{F^{2}}{f_{a}^{2}}\frac{M_{0}^{2}m_{\overline{\pi}}^{2}(2m_{\overline{K}}^{2}-m_{\overline{\pi}}^{2})}{8M_{0}^{2}m_{\overline{K}}^{2}-2m_{\overline{\pi}}^{2}(M_{0}^{2}-6m_{\overline{K}}^{2})-6m_{\overline{\pi}}^{4}}\,,$

(30)

which reduces to the well celebrated result $m_{\overline{a}}^{2}=F^{2}m_{\overline{\pi}}^{2}/(4f_{a}^{2})$ [2] by keeping the leading expansions of $m_{\overline{\pi}}^{2}/m_{\overline{K}}^{2}$ and $m_{\overline{\pi}}^{2}/M_{0}^{2}$.

For the LO fit in Ref. [44], the value of $M_{0}$ is determined to be $M_{0}=820.0$ MeV, which leads to

$\displaystyle m_{\overline{\eta}}=493.6~{}{\rm MeV},\quad m_{\overline{\eta}^{\prime}}=957.7~{}{\rm MeV},\quad\theta=-19.6^{\circ}\,.$

(31)

By substituting these values to Eqs. (LABEL:eq.v14)-(28) with $m_{a,0}=0$, we obtain

$\displaystyle m_{\overline{a}}=6.1~{}{\rm\mu eV}\frac{10^{12}{\rm GeV}}{f_{a}}~{}\,,\quad v_{14}=\frac{-0.012}{f_{a}}\,,\quad v_{24}=\frac{-0.035}{f_{a}}\,,\quad v_{34}=\frac{-0.026}{f_{a}}\,,$

(32)

where the mass agrees well with the recent determinations [29, 45]. The units of the numbers in the numerators of $v_{14/24/34}$ are taken as GeV here. Nevertheless, it is clear that the LO predictions for the masses of $\eta$ and $\eta^{\prime}$ in Eq. (31) are still not satisfactory when compared to their physical values. In naively estimating $m_{\overline{\eta}}$ and $m_{\overline{\eta}^{\prime}}$ by the physical masses of $\eta$ and $\eta^{\prime}$, the axion mass can be predicted from Eq. (29) with the result

$\displaystyle m_{\overline{a}}=9.7~{}{\rm\mu eV}\frac{10^{12}{\rm GeV}}{f_{a}}~{}\,,\quad v_{14}=\frac{-0.011}{f_{a}}\,,\quad v_{24}=\frac{-0.028}{f_{a}}\,,\quad v_{34}=\frac{-0.026}{f_{a}}\,,$

(33)

with $M_{0}=820$ MeV and $\theta=-19.6^{\circ}$. In turn, if we take $M_{0}=820$ MeV and $\theta=-10.0^{\circ}$ the result turns out to be

$\displaystyle m_{\overline{a}}=14.5~{}{\rm\mu eV}\frac{10^{12}{\rm GeV}}{f_{a}}~{}\,,\quad v_{14}=\frac{-0.008}{f_{a}}\,,\quad v_{24}=\frac{-0.015}{f_{a}}\,,\quad v_{34}=\frac{-0.027}{f_{a}}\,.$

(34)

The obvious changes of the predictions in Eqs. (32)-(34) by taking different phenomenological inputs for the $m_{\bar{\eta}}$ and $m_{\bar{\eta}^{\prime}}$, 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.

Up to the NLO in the $\delta$ counting, the bilinear terms involving the $\overline{\pi}^{0},\overline{\eta},\overline{\eta}^{\prime}$ and the axion field $\overline{a}$, which are diagonalized at LO, can be generally written as

	$\displaystyle\mathcal{L}=$	$\displaystyle\dfrac{1+\delta^{\eta}_{k}}{2}\partial_{\mu}\overline{\eta}\partial^{\mu}\overline{\eta}+\dfrac{1+\delta^{\eta^{\prime}}_{k}}{2}\partial_{\mu}\overline{\eta}^{\prime}\partial^{\mu}\overline{\eta}^{\prime}+\delta_{k}^{\eta\eta^{\prime}}\partial_{\mu}\overline{\eta}\partial^{\mu}\overline{\eta}^{\prime}-\dfrac{m^{2}_{{\overline{\eta}}}+\delta_{m^{2}_{{\eta}}}}{2}\overline{\eta}\,\overline{\eta}-\dfrac{m^{2}_{{\overline{\eta}^{\prime}}}+\delta_{m^{2}_{{\eta^{\prime}}}}}{2}\overline{\eta}^{\prime}\,\overline{\eta}^{\prime}-\delta_{m^{2}}^{\eta\eta^{\prime}}\overline{\eta}\,\overline{\eta}^{\prime}$		(35)
		$\displaystyle+\dfrac{1+\delta^{\pi}_{k}}{2}\partial_{\mu}\overline{\pi}^{0}\partial^{\mu}\overline{\pi}^{0}+\delta^{\pi\eta}_{k}\partial_{\mu}\overline{\pi}^{0}\partial^{\mu}\overline{\eta}+\delta^{\pi\eta^{\prime}}_{k}\partial_{\mu}\overline{\pi}^{0}\partial^{\mu}\overline{\eta}^{\prime}-\dfrac{m_{\overline{\pi}}^{2}+\delta_{m^{2}_{\pi}}}{2}\overline{\pi}^{0}\,\overline{\pi}^{0}-\delta^{\pi\eta}_{m^{2}}\overline{\pi}^{0}\overline{\eta}-\delta^{\pi{\eta}^{\prime}}_{m^{2}}\overline{\pi}^{0}\overline{\eta}^{\prime}$
		$\displaystyle+\dfrac{1+\delta^{a}_{k}}{2}\partial_{\mu}\overline{a}\partial^{\mu}\overline{a}+\delta^{a\pi}_{k}\partial_{\mu}\overline{a}\partial^{\mu}\overline{\pi}^{0}+\delta^{a\eta}_{k}\partial_{\mu}\overline{a}\partial^{\mu}\overline{\eta}+\delta^{a{\eta}^{\prime}}_{k}\partial_{\mu}\overline{a}\partial^{\mu}\overline{\eta}^{\prime}-\dfrac{m^{2}_{\overline{a}}+\delta_{m^{2}_{a}}}{2}\overline{a}\,\overline{a}-\delta^{a\pi}_{m^{2}}\overline{a}\,\overline{\pi}^{0}$
		$\displaystyle-\delta^{a\eta}_{m^{2}}\overline{a}\,\overline{\eta}-\delta^{a{\eta}^{\prime}}_{m^{2}}\overline{a}\,\overline{\eta}^{\prime}\,,$

where the $\delta_{k}^{X}$ with subscript $k$ denotes the corrections from higher orders to the kinetic terms and other $\delta_{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 $\eta$-$\eta^{\prime}$ mixing in Ref. [42, 44], we extend the calculations by simultaneously including the $\pi^{0}$ and the axion $a$ in the above equation. All the $\delta_{i}$ parameters in Eq. (35) 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 $\delta$ 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

	$\displaystyle\left(\begin{array}[]{c}\hat{\pi}^{0}\\ \hat{\eta}\\ \hat{\eta}^{\prime}\\ \hat{a}\end{array}\right)\,=\,$	$\displaystyle\left(\begin{array}[]{cccc}1&-y_{12}&-y_{13}&-y_{14}\\ y_{12}&1&-y_{23}&-y_{24}\\ y_{13}&y_{23}&1&-y_{34}\\ y_{14}&y_{24}&y_{34}&1\end{array}\right)\times$		(44)
		$\displaystyle\left(\begin{array}[]{cccc}1-x_{11}&-x_{12}&-x_{13}&-x_{14}\\ -x_{12}&1-x_{22}&-x_{23}&-x_{24}\\ -x_{13}&-x_{23}&1-x_{33}&-x_{34}\\ -x_{14}&-x_{24}&-x_{34}&1-x_{44}\end{array}\right)\,\left(\begin{array}[]{c}\overline{\pi}^{0}\\ \overline{\eta}\\ \overline{\eta}^{\prime}\\ \overline{a}\end{array}\right)\,,$		(53)

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

	$\displaystyle x_{11}=-\frac{\delta^{\pi}_{k}}{2}\,,\quad x_{12}=-\frac{\delta_{k}^{\pi\eta}}{2}\,,\quad x_{13}=-\frac{\delta_{k}^{\pi\eta^{\prime}}}{2}\,,\quad x_{14}=-\frac{\delta_{k}^{a\pi}}{2}\,,\quad x_{22}=-\frac{\delta^{\eta}_{k}}{2}\,,$
	$\displaystyle x_{23}=-\frac{\delta_{k}^{\eta\eta^{\prime}}}{2}\,,\quad x_{24}=-\frac{\delta_{k}^{a\eta}}{2}\,,\quad x_{33}=-\frac{\delta^{\eta^{\prime}}_{k}}{2}\,,\quad x_{34}=-\frac{\delta_{k}^{a\eta^{\prime}}}{2}\,,\quad x_{44}=-\frac{\delta^{a}_{k}}{2}\,,$		(54)

and the expressions of the $y_{ij}$ read

	$\displaystyle y_{12}=\frac{\delta_{m^{2}}^{\pi\eta}+x_{12}(m_{\overline{\eta}}^{2}+m_{\overline{\pi}}^{2})}{m_{\overline{\eta}}^{2}-m_{\overline{\pi}}^{2}}\,,\quad y_{13}=\frac{\delta_{m^{2}}^{\pi\eta^{\prime}}+x_{13}(m_{\overline{\eta}^{\prime}}^{2}+m_{\overline{\pi}}^{2})}{m_{\overline{\eta}^{\prime}}^{2}-m_{\overline{\pi}}^{2}}\,,\quad y_{14}=\frac{\delta_{m^{2}}^{a\pi}+x_{14}(m_{a,0}^{2}+m_{\overline{\pi}}^{2})}{m_{a,0}^{2}-m_{\overline{\pi}}^{2}}\,,$
	$\displaystyle y_{23}=\frac{\delta_{m^{2}}^{\eta\eta^{\prime}}+x_{23}(m_{\overline{\eta}}^{2}+m_{\overline{\eta}^{\prime}}^{2})}{m_{\overline{\eta}^{\prime}}^{2}-m_{\overline{\eta}}^{2}}\,,\quad y_{24}=\frac{\delta_{m^{2}}^{a\eta}+x_{24}(m_{\overline{\eta}}^{2}+m_{a,0}^{2})}{m_{a,0}^{2}-m_{\overline{\eta}}^{2}}\,,\quad y_{34}=\frac{\delta_{m^{2}}^{a\eta^{\prime}}+x_{34}(m_{\overline{\eta}^{\prime}}^{2}+m_{a,0}^{2})}{m_{a,0}^{2}-m_{\overline{\eta}^{\prime}}^{2}}\,.$

By keeping the contributions up to NLO, the masses of the diagonalized canonical states take the form

	$\displaystyle m_{\hat{\pi}}^{2}$	$\displaystyle=$	$\displaystyle m_{\overline{\pi}}^{2}+\delta_{m_{\pi}^{2}}-m_{\overline{\pi}}^{2}\delta^{\pi}_{k}\,,$		(56)
	$\displaystyle m_{\hat{\eta}}^{2}$	$\displaystyle=$	$\displaystyle m_{\overline{\eta}}^{2}+\delta_{m_{\eta}^{2}}-m_{\overline{\eta}}^{2}\delta^{\eta}_{k}\,,$		(57)
	$\displaystyle m_{\hat{\eta}^{\prime}}^{2}$	$\displaystyle=$	$\displaystyle m_{\overline{\eta}^{\prime}}^{2}+\delta_{m_{\eta^{\prime}}^{2}}-m_{\overline{\eta}^{\prime}}^{2}\delta^{\eta^{\prime}}_{k}\,,$		(58)
	$\displaystyle m_{\hat{a}}^{2}$	$\displaystyle=$	$\displaystyle m_{\overline{a}}^{2}+\delta_{m_{a}^{2}}-m_{a,0}^{2}\delta^{a}_{k}\,,$		(59)

and the kaon mass squared up to NLO is similarly given by

$$m_{\hat{K}}^{2}=m_{\overline{K}}^{2}+\delta_{m_{K}^{2}}-m_{\overline{K}}^{2}\delta^{K}_{k}\,.$$

(60)

Up to NLO the relations between the physical states (denoted by the hatted fields) and the LO diagonalized ones reduce to

$\displaystyle\left(\begin{array}[]{c}\hat{\pi}^{0}\\ \hat{\eta}\\ \hat{\eta}^{\prime}\\ \hat{a}\end{array}\right)\,=\,$

$\displaystyle\left(\begin{array}[]{cccc}1-x_{11}&-x_{12}-y_{12}&-x_{13}-y_{13}&-x_{14}-y_{14}\\ -x_{12}+y_{12}&1-x_{22}&-x_{23}-y_{23}&-x_{24}-y_{24}\\ -x_{13}+y_{13}&-x_{23}+y_{23}&1-x_{33}&-x_{34}-y_{34}\\ -x_{14}+y_{14}&-x_{24}+y_{24}&-x_{34}+y_{34}&1-x_{44}\end{array}\right)\left(\begin{array}[]{c}\overline{\pi}^{0}\\ \overline{\eta}\\ \overline{\eta}^{\prime}\\ \overline{a}\end{array}\right)\,.$

(73)