Topological Susceptibility and QCD Axion Mass: QED and NNLO corrections

TL;DR

This paper refines the calculation of the QCD axion mass by incorporating electromagnetic and NNLO corrections into the topological susceptibility, demonstrating good convergence and potential for high-precision determination.

Contribution

It introduces the first inclusion of $O( ext{em})$ and NNLO corrections in the chiral expansion for topological susceptibility, enhancing the precision of the QCD axion mass estimate.

Findings

Electromagnetic correction is 0.65(21)%

NNLO correction is -0.71(29)%

Chiral expansion shows good convergence

Abstract

We improve the precision of the topological susceptibility of QCD, and therefore of the QCD axion mass, by including $O(\alpha_{\rm em})$ and NNLO corrections in the chiral expansion, which amount to 0.65(21)% and -0.71(29)% respectively. Both corrections are one order of magnitude smaller than the known NLO ones, confirming the very good convergence of the chiral expansion and its reliability. Using the latest estimates for the light quark masses the current uncertainty is dominated by the one of the low-energy constant $\ell_7$. When combined with possible improvements on the light quark mass ratio and $\ell_7$ from lattice QCD, our computation could allow to determine the QCD axion mass with per-mille accuracy.

Figures

Tables

Table 1. Table 1: Numerical values of the n f = 2 subscript 𝑛 𝑓 2 n_{f}=2 EM LECs k i r superscript subscript 𝑘 𝑖 𝑟 k_{i}^{r} at the scale μ ¯ = 770 ¯ 𝜇 770 \bar{\mu}=770  MeV extracted using their relation to K i r superscript subscript 𝐾 𝑖 𝑟 K_{i}^{r} . To the k i r superscript subscript 𝑘 𝑖 𝑟 k_{i}^{r} is assigned a conservative 100 % percent 100 100\% uncertainty.

$k_1^r$	$k_2^r$	$k_3^r$	$k_4^r$	$k_7^r$
$8.4$	$3.4$	$2.7$	$1.4$	$2.2$	$\times {10}^{-3}$

Table 2. Table 2: Top: Numerical values of the c ¯ i subscript ¯ 𝑐 𝑖 \bar{c}_{i} constants appearing in δ Γ loc superscript subscript 𝛿 Γ loc \delta_{\Gamma}^{\rm loc} from  [ 34 ] . Bottom: Numerical values of the n f = 3 subscript 𝑛 𝑓 3 n_{f}=3 radiative and leptonic LECs from  [ 32 , 33 , 35 ] at the scale μ ¯ = 770 ¯ 𝜇 770 \bar{\mu}=770  MeV. All constants are assigned a conservative 100 % percent 100 100\% uncertainty.

${\overline{c}}_2$	${\overline{c}}_3$	${\overline{c}}_4$	${\overline{c}}_{2t}$
$5.2$	$-10.5$	$1.69$	$0$

Table 3. Table 3: Numerical values of the n f = 2 subscript 𝑛 𝑓 2 n_{f}=2 LECs in δ 2 subscript 𝛿 2 \delta_{2} at the scale μ ¯ = 770 ¯ 𝜇 770 \bar{\mu}=770  MeV extracted by combining the lattice results of  [ 27 ] and the matching with the n f = 3 subscript 𝑛 𝑓 3 n_{f}=3 LECs of  [ 45 ] (see app.  B for more details).

${10}^6{c}_6^r$	${10}^6{c}_7^r$	${10}^6{c}_8^r$	${10}^6{c}_9^r$	${10}^6{c}_{10}^r$	${10}^6{c}_{11}^r$	${10}^6{c}_{17}^r$	${10}^6{c}_{18}^r$	${10}^6{c}_{19}^r$
$1.0\left(3.8\right)$	$-2.2\left(2.0\right)$	$2.1\left(1.1\right)$	$-0.3\left(1.2\right)$	$-1.0\left(1.1\right)$	$-0.1\left(0.7\right)$	$6.4\left(3.8\right)$	$-2.2\left(5.9\right)$	$-0.1\left(9.5\right)$

Table 4. Table 4: Numerical values of the combinations of SU(2) LECs at the scale μ ¯ = 770 ¯ 𝜇 770 \bar{\mu}=770  MeV extracted from the 450 450 450 MeV cut-fit of the partially quenched simulations in  [ 45 ] (tab. 6).

$c_6^r-c_{17}^r$	$-5.33\left(77\right)\cdot {10}^{-6}$
$2c_6^r-12c_{10}^r+18c_{11}^r-c_{18}^r$	$14.5\left(3.9\right)\cdot {10}^{-6}$
$c_7^r$	$-3.9\left(2.3\right)\cdot {10}^{-6}$
$c_8^r$	$0.0\left(1.8\right)\cdot {10}^{-6}$
$2c_7^r+4c_8^r$	$6.2\left(3.2\right)\cdot {10}^{-6}$
$c_9^r$	$-0.2\left(1.2\right)\cdot {10}^{-6}$
$c_{10}^r$	$-1.0\left(1.1\right)\cdot {10}^{-6}$
$19c_{11}^r-12c_{10}^r$	$10.1\left(3.1\right)\cdot {10}^{-6}$

Table 5. Table 5: Numerical value of the NLO and NNLO couplings L i r superscript subscript 𝐿 𝑖 𝑟 L_{i}^{r} and C i r superscript subscript 𝐶 𝑖 𝑟 C_{i}^{r} at the scale μ ¯ = 770 ¯ 𝜇 770 \bar{\mu}=770  MeV. We associated 100 % percent 100 100\% uncertainty to the C i r superscript subscript 𝐶 𝑖 𝑟 C_{i}^{r} .

$L_1^r$	$L_2^r$	$L_3^r$	$L_4^r$	$L_5^r$	$L_6^r$	$L_7^r$	$L_8^r$
$0.5\left(2\right)$	$0.8\left(3\right)$	$-3.1\left(1.0\right)$	$0.09\left(34\right)$	$1.19\left(40\right)$	$0.16\left(20\right)$	$-0.34\left(11\right)$	$0.55\left(18\right)$	$\times {10}^{-3}$

Table 6. Table 6: Numerical input values used in the computations. Dimensionful quantities are given in MeV.

$z$	0.472(11)	eq. (4)	$w^{-1}$	27	[26]
$f_{\pi }$	92.3	[28]	${\overline{\ell }}_3$	2.81(49)	[27]
$m_{{\pi }^0}$	134.98	[28]	${\overline{\ell }}_4$	4.02(25)	[27]
$m_{{\pi }^{+}}$	139.57	[28]	$h_1^r-h_3^r-{\ell }_4^r$	-0.0049(12)	eq. (5)
$m_K$	$495$	[28]	${\ell }_7$	0.0065(38)	[27]
$m_{\rho }$	775	[28]	${\alpha }_{\mathrm{em}}^{-1}$	$137$	[28]
$m_{\mu }$	105.658	[28]	${\Gamma }_{{\pi }^{+}\to \mu \nu \left(\gamma \right)}$	$2.5281\cdot {10}^{-14}$	[28]
$G_F$	$1.16638\cdot {10}^{-11}$	[28]	$V_{ud}$	$0.9742$	[28]

Equations

$${\chi }_{\mathrm{top}}^{\mathrm{LO}}=\frac{z}{(1+z{)}^2}{m}_{{\pi }^0}^2{f}_{\pi }^2,z\equiv \frac{m_u}{m_d},$$

$${\chi }_{\mathrm{top}}^{\mathrm{NLO}}=$$

$${\delta }_1=$$

$${\chi }_{\mathrm{top}}^{1/4}=75.5(5)\mathrm{MeV},m_a=5.70(6{)}_z(4{)}_{{\ell }_i^r}\mu \mathrm{eV}\frac{10^{12}\mathrm{GeV}}{f_a},$$

$$z={\left(\frac{m_u}{m_d}\right)}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=0.472(11),$$

$$h_1^r-h_3^r-{\ell }_4^r=-0.0049(12),{\ell }_7=0.0065(38).$$

$${\delta }_1=-0.042(13),$$

$${\chi }_{\mathrm{top}}^{1/4}=75.46(29)\mathrm{MeV},m_a=5.69(2{)}_z(4{)}_{{\ell }_i^r}\mu \mathrm{eV}\frac{10^{12}\mathrm{GeV}}{f_a}.$$

$$\frac{\partial \log z}{\partial \log \mu }=\frac{6{\alpha }_{\mathrm{em}}}{4\pi }\left[{\left(\frac{2}{3}\right)}^2-{\left(-\frac{1}{3}\right)}^2\right]=\frac{{\alpha }_{\mathrm{em}}}{2\pi }$$

$${\chi }_{\mathrm{top}}$$

$${\delta }_e$$

$$Z=\frac{m_{{\pi }^{+}}^2-m_{{\pi }^0}^2}{2e^2{f}_{{\pi }^{+}}^2}+\cdots \simeq 0.81.$$

$${\delta }_e=0.0065(21).$$

$$\mu \frac{\partial }{\partial \mu }{k}_5^r=-\frac{3}{5}\frac{1}{(4\pi {)}^2},\mu \frac{\partial }{\partial \mu }{k}_7^r=-\frac{3}{4}\frac{1}{(4\pi {)}^2}.$$

$$\Gamma_{\pi^{+}\to\mu\nu(\gamma)}=\frac{G_{F}^{2}|V_{ud}|^{2}m_{\pi^{+}}m_{\mu}^{2}f_{\pi^{+}}^{2}}{4\pi}\Bigl{(}1-\frac{m_{\mu}^{2}}{m_{\pi^{+}}^{2}}\Bigr{)}^{2}\,\left[1+\delta_{\Gamma}^{\rm loc}+\delta_{\Gamma}^{\rm had}\right]$$

$${\delta }_{\Gamma }^{\mathrm{loc}}$$

$${\delta }_{\Gamma }^{\mathrm{had}}$$

$$+\frac{1}{(4\pi {)}^2}\left(2-3Z-Z\log \left(\frac{m_K^2}{{\bar{\mu }}^2}\right)+(3-2Z)\log \left(\frac{m_{\pi }^2}{{\bar{\mu }}^2}\right)\right)],$$

$$F(x)$$

$${\delta }_{\Gamma }^{\mathrm{loc}}+{\delta }_{\Gamma }^{\mathrm{had}}=0.0177(38),$$

$$\chi_{\rm top}=\frac{z}{(1+z)^{2}}m_{\pi^{0}}^{2}\frac{\Gamma_{\pi^{+}\to\mu\nu(\gamma)}}{\frac{G_{F}^{2}|V_{ud}|^{2}m_{\pi^{+}}m_{\mu}^{2}}{4\pi}\Bigl{(}1-\frac{m_{\mu}^{2}}{m_{\pi^{+}}^{2}}\Bigr{)}^{2}}\left[1+\delta_{e}-\delta_{\Gamma}^{\rm loc}-\delta_{\Gamma}^{\rm had}+\dots\right]\,.$$

$${\delta }_e-{\delta }_{\Gamma }^{\mathrm{had}}$$

$$=e^2[2K_4^r-4K_3^r+\frac{8}{3}\Delta (K_9^r+K_{10}^r)+\frac{4}{3}{X}_1^r+4(X_2^r-X_3^r)+{\tilde{X}}_6^{r,\mathrm{eff}}$$

$$-\frac{1}{(4\pi {)}^2}\left(2(1+Z+\Delta Z)+(3+2Z)\log \frac{m_{\pi }^2}{{\bar{\mu }}^2}+2\Delta Z\log \frac{m_K^2}{{\bar{\mu }}^2}\right)].$$

$${\delta }_e-{\delta }_{\Gamma }^{\mathrm{had}}-{\delta }_{\Gamma }^{\mathrm{loc}}=0.024(6),$$

$${\delta }_{\Gamma }^{\mathrm{loc}}+{\delta }_{\Gamma }^{\mathrm{had}}=0.0169(15),$$

$${\chi }_{\mathrm{top}}^{\mathrm{NNLO}}$$

$${\delta }_2$$

$$+\frac{h_1^r-h_3^r-{\ell }_4^r-{\ell }_7}{8{\pi }^2}{\bar{\ell }}_3-\frac{1}{(4\pi {)}^4}\left({\bar{\ell }}_4^2-{\bar{\ell }}_3{\bar{\ell }}_4+\frac{25}{16}\right)$$

$$\displaystyle\left.\qquad\qquad+\Biggl{(}32\left(c_{9}^{r}-6c_{10}^{r}-4\left(c_{11}^{r}+c_{17}^{r}+c_{18}^{r}\right)-6c_{19}^{r}\right)-\frac{7\ell_{7}}{8\pi^{2}}\log\frac{m_{\pi}^{2}}{\bar{\mu}^{2}}\right.$$

$$\displaystyle\left.\qquad\qquad+8(h_{1}^{r}-h_{3}^{r}-\ell_{4}^{r})\ell_{7}-12\ell_{7}^{2}+\frac{2\bar{\ell}_{3}+2\bar{\ell}_{4}-1}{8\pi^{2}}\ell_{7}\Biggr{)}\Delta^{2}+20\ell_{7}^{2}\,\Delta^{4}\right]\,,$$

$${\delta }_2=-0.0071(01{)}_z(23{)}_{{\ell }_i^r}(19{)}_{c_i^r}=-0.0071(29).$$

$${\chi }_{\mathrm{top}}=\frac{z}{(1+z{)}^2}{m}_{{\pi }^0}^2{f}_{{\pi }^{+}}^2\left[1+{\delta }_1+{\delta }_2+{\delta }_e\right],$$

$${\chi }_{\mathrm{top}}^{1/4}=75.44(34)\mathrm{MeV},m_a=5.691(51)\mu \mathrm{eV}\frac{10^{12}\mathrm{GeV}}{f_a}.$$

$$m_{a}=\Bigl{[}\underbrace{5.815(22)_{z}(04)_{f_{\pi}}}_{\text{LO}}\underbrace{-0.121(38)_{\ell^{r}_{i}}}_{\text{NLO}}\underbrace{-0.022(07)_{\ell^{r}_{i}}(05)_{c_{i}^{r}}}_{\text{NNLO}}\underbrace{+0.019(06)_{k_{i}^{r}}}_{\text{EM}}\Bigr{]}\,\mu{\rm eV}\,\frac{10^{12}~{}{\rm GeV}}{f_{a}}\,,$$

$$m_{{\pi }^0}^2=$$

Full text

1. 1 Topological Susceptibility at LO and NLO

2. 2 QED corrections

3. 3 NNLO corrections

4. 4 Final Results and Axion Mass

5. A Results in terms of the quark masses

6. A.1 Two-flavor results

7. A.2 Three-flavor results

8. B Extraction of the NNLO LECs and input parameters

1 Topological Susceptibility at LO and NLO

In QCD the topological susceptibility $\chi_{\rm top}$ is one of the fundamental observables describing the non-trivial properties of the QCD vacuum. Defined as the second derivative of the free energy with respect to the $\theta$-angle at $\theta=0$, it determines how the QCD vacuum energy depends on the CP-violating $\theta$ parameter. While $\chi_{\rm top}$ vanishes at all order in perturbation theory, at high temperatures its value is expected to be well reproduced by semi-classical small-instanton configurations.111 Several lattice simulations in pure Yang-Mills [1, 2] and in QCD [3, 4, 5, 6] indicate that the temperature behavior predicted by the instanton gas approximation might be valid already at small temperatures, just above the QCD transition, although the overall size of $\chi_{\rm top}$ is not yet well reproduced. At zero temperatures such a description is not valid, and indeed in the chiral limit current algebra relates $\chi_{\rm top}$ to the chiral condensate, which is notoriously associated to renormalons rather than to a semi-classical field configuration.

Recently the determination of $\chi_{\rm top}$ has seen a wave of renewed interest. Indeed, the most plausible known solution to the strong-CP problem [7] involves the presence of a light pseudo-Goldstone boson, the QCD axion [8, 9, 10, 11, 12, 13], whose mass is determined by $\chi_{\rm top}$ through the relation $m_{a}^{2}=\chi_{\rm top}/f_{a}^{2}$ (where $f_{a}$ is the Peccei-Quinn scale controlling the axion coupling to the Standard Model). Since the existence of the QCD axion could also explain the dark matter abundance in our Universe [14, 15, 16], a multitude of experiments are being pursued to search for this particle (see e.g. [17, 18]). Using various forms of resonance effects to amplify the otherwise too feeble signal, several of these experiments would be able to measure the axion mass with very high precision, even down to ${\cal O}(10^{-6})$. When combined with measurements of the axion couplings and possibly the information of the axion relic abundance, such precision could be used to learn about the dynamics of new physics at much higher scales, as well as physics of the early universe including inflation, reheating and pre-BBN evolution.

At a first look a high precision determination of $\chi_{\rm top}$ seems hampered by its non-perturbative nature. On the contrary, chiral effective theories are particularly powerful in this case. By exploiting the freedom to rotate the whole $\theta$-dependence of the QCD functional into the phases of the lightest quark masses, the $\theta$-dependence of low-energy QCD observables (e.g. the vacuum energy) can be computed analytically using chiral Lagrangians, which are expansions in the light quark masses. In particular using the lightest two quarks, the effective expansion parameter, $m_{u,d}/m_{s}\sim$ few percent, is rather small, suggesting a fast convergence.

Indeed, as shown explicitly in [19], the leading order formula [8]

[TABLE]

where $m_{\pi^{0}}$ is the physical neutral pion mass and $f_{\pi}$ its decay constant (normalized as $f_{\pi}\sim 92.3$ MeV), is accurate at the few percent level. Next to leading order (NLO) corrections in the chiral expansion have been computed in [19], and result in

[TABLE]

where the coefficients $h_{i}^{r}$ and $\ell_{i}^{r}$ are the low-energy constants (LECs) of the chiral Lagrangian defined in ref. [20]. Using the available estimates for the quark mass ratio $z=0.48(3)$ and for the LECs, the topological susceptibility and the corresponding axion mass were estimated to be

[TABLE]

where the uncertainties in $m_{a}$ come respectively from $z$ and the LECs. This estimate represents the current state of the art for the topological susceptibility and the QCD axion mass.

Since the results in ref. [19] new lattice simulations became available. Direct measurements of the topological susceptibility were performed in the isospin limit $z=1$ in refs. [21] and [4], which, once corrected for the leading isospin-breaking effects (i.e. the factor $z^{1/4}\sqrt{2/(1+z)}$), give respectively $\chi^{1/4}_{\rm top}=\,$75(3) and 75(2) MeV, in nice agreement with eq. (3). While the current precision of these estimates is still roughly a factor of four worse than (3), the situation may change in the near future as systematics are further reduced on the lattice side.

At the same time new lattice estimates of the quark mass ratio $z$ appeared and improved the previous ones, namely $z=0.485(20)$ [22] with three dynamical quarks and $z=0.513(31)$ [23] and $z=0.453(16)$ [24], with four dynamical quarks. By combining them with the older $z=0.470(56)$ [25] (also with four dynamical quarks), we get the following improved estimate

[TABLE]

which agrees with the previous estimate ($z=0.48(3)$) used in ref. [19] and improves its precision by a factor of three. We should warn the reader that here the error has been computed by simply propagating the uncertainties quoted by each collaboration, since a proper combination is not yet available. In the remainder of the paper we will use the value in eq. (4) as reference, however we will always report separately the uncertainty originating from $z$ and the total one so that it can easily be rescaled if needed.

For the LECs appearing in eq. (1) we now use the values

[TABLE]

The first combination is computed using the matching to $L_{8}^{r}$ as described in ref. [19] and using the latest FLAG estimate $L_{8}^{r}=0.00055(15)$ [26], while the second value is taken from the direct lattice simulation of ref. [27]. These values give

[TABLE]

where the error is dominated by the one from $\ell_{7}$. Combining everything together we get the updated values for the topological susceptibility and the axion mass at NLO

[TABLE]

Given the improved value for the light quark mass ratio, the dominant error now became the one from the NLO LECs, in particular $\ell_{7}$, which also controls the strong isospin breaking effect in the pion mass splitting, indeed poorly known. An improvement on this quantity would directly translate into an equivalent improvement in our knowledge of $\chi_{\rm top}$ and thus $m_{a}$. Conversely, improvements in the direct computation of $\chi_{\rm top}$ on the lattice could be used to better determine both $z$ and $\ell_{7}$.

A natural question to ask is how much an advance in our knowledge of the light quark masses and the NLO LECs can increase the precision of $\chi_{\rm top}$, before other unknown corrections need to be considered. Among the latter, the most relevant are the NNLO corrections of the chiral expansion and ${\cal O}(\alpha_{\rm em})$ electromagnetic (EM) corrections. The firsts do not only determine the ultimate precision reachable with eq. (1) but also measure the convergence and reliability of the chiral expansion. Of course the size of the NNLO corrections is only relevant in the chiral expansion approach and does not represent a source of uncertainty for lattice simulations222On the other hand, lattice simulations have to face a number of systematic uncertainties which are not present in the chiral expansion such as finite volume, finite lattice spacing effects, explicit chiral symmetry breaking, etc., some of which require delicate and careful analyses., which contain the full non-perturbative result. The EM corrections, on the other hand, are common to both approaches and so far have never been considered. As we will show in the next section, their size is smaller with the choice made in eq. (1) of using the value of the neutral pion mass in the LO formula. Even with this choice, however, the value of the EM corrections is just below the size of the present uncertainties for $\chi_{\rm top}$, which means that further improvements cannot ignore them.

In the rest of the paper we will present first the analysis of the EM corrections to the topological susceptibility in sec. 2, and then the NNLO ones in sec. 3. We will combine all of them in sec. 4 where we will also give our final estimate for the axion mass with a discussion of the various sources of uncertainties. In appendix A we report the formulas for our results with the explicit quark mass dependence, suitable to be used in lattice simulation fits. Finally, in appendix B we give all the details of the numerical extraction of the values of the LECs used in the results.

2 QED corrections

While the QCD axion has a vanishing electric charge, its mass can receive ${\cal O}(\alpha_{\rm em})$ corrections from several sources. Indeed, the leading order formula (1) involves a number of quantities that can introduce potentially large EM corrections depending on the way they are defined and extracted by experiments.

• •

The pion masses for the neutral and the charged states are degenerate at leading order, but differ at higher orders due to isospin and EM effects. The latter largely dominate this difference, which amounts to $m_{\pi^{+}}-m_{\pi^{0}}=4.5936(5)$ MeV, i.e. around 4% of the total mass. The main effect comes from the charged pion mass, whose corrections are ${\cal O}(e^{2})$, while those in the neutral pion mass start at ${\cal O}(e^{2}p^{2})$. Therefore, depending on which pion mass is used in eq. (1), the axion mass can vary by 4%, which is more than the quoted uncertainties of the previous section. As we will see below, the naive expectation that the neutral pion mass should be used to minimize EM effects is the correct one. Indeed the pion mass entering in the leading order formula can be understood as arising from the mixing between the axion and the neutral pion state.

• •

In QCD the pion decay constant $f_{\pi}$ is not unambiguously defined when EM interactions are turned on. In chiral perturbation theory, on the other hand, $\alpha_{\rm em}$ can be controlled analytically and it is possible to define $f_{\pi}$ unambiguously. The best determination of $f_{\pi}$ at the moment comes from (radiative) leptonic pion decays $\pi^{+}\to\mu{\nu_{\mu}}(\gamma)$ where both experimental and theoretical uncertainties are small [28]. As we will discuss in more detail below, the EM corrections to $\Gamma_{\pi^{+}\to\mu\nu(\gamma)}$ are dominated by a calculable short distance contribution. The long distance hadronic contribution (which is of the same order of the EM corrections we want to compute for the axion mass) is subleading but dominates the current error of $f_{\pi}$. Given the importance of such corrections for our computation, we revisit their estimate and analyze their interplay with the genuine corrections to the axion mass. An alternative determination of $f_{\pi}$ could be obtained from the neutral pion decay $\pi^{0}\to\gamma\gamma$, however both the theoretical and experimental uncertainties are not competitive with the charged pion channel [28].

• •

While the light quark mass ratio $z=m_{u}/m_{d}$ at leading order is renormalization group (RG) invariant with respect to QCD corrections, it is not with respect to the QED ones [29]. This introduces an ${\cal O}(\alpha_{\rm em})$ ambiguity in the tree-level formula of the axion mass that should be removed by the sub-leading EM corrections:

[TABLE]

A change of ${\cal O}(1)$ in the renormalization scale introduces a shift of ${\cal O}(10^{-3})$ in $z$ that can be taken as a lower bound to the order of magnitude of the expected EM corrections to the axion mass.

We start by reporting the result for the computation of the leading EM corrections to the topological susceptibility, which begin at ${\cal O}(e^{2}p^{2})$ in the chiral expansion once the leading order term is written in terms of the physical333Whenever $m_{\pi}$ appears in the following formulas, it can be equivalently understood as $m_{\pi^{0}}$ or $m_{\pi^{+}}$ because the difference will be accounted by higher orders in either $e^{2}$ or $p^{2}$. For the numerical estimates we used $m_{\pi}=m_{\pi^{0}}$. neutral pion mass $m_{\pi^{0}}$ (including EM corrections) and the physical charged pion decay constant $f_{\pi^{+}}$ (defined in pure QCD, i.e. at $\alpha_{\rm em}=0$):

[TABLE]

where dots in eq. (9) represent the non-EM corrections discussed in secs. 1 and 3, the coefficients $Z$ and $k_{i}^{r}$ are the $n_{f}=2$ EM low-energy constants from [30], and $\bar{\mu}$ is the renormalization scale of the chiral Lagrangian, whose dependence cancels against that from the $k^{r}_{i}$ coefficients. As anticipated before, once the LO formula is written in terms of $m_{\pi^{0}}^{2}$, the EM corrections starts at ${\cal O}(e^{2}p^{2})$ (the $\delta_{e}$ term). In particular the EM pion mass splitting effects parametrized by

[TABLE]

are loop suppressed. Although the value for the couplings $k_{i}^{r}$ is not known directly, it can been inferred, as in [31], using their relation to the $n_{f}=3$ constants $K_{i}^{r}$, which have been estimated in refs. [32, 33] using various techniques including sum rules and vector meson dominance. The values for the $k_{i}^{r}$ we use are taken from [31] (with $K_{9}^{r}$ from [33]) and reported in tab. 1.

Because of the model dependence of such estimates we decided to assign a conservative 100% uncertainty on each LEC, i.e. we use the mentioned values as an order of magnitude estimate of their size. Substituting the numerical values we find

[TABLE]

While we have assigned 100% uncertainties to the LECs $k_{i}^{r}$, the uncertainty on $\delta_{e}$ only amounts to 30% because the dominant contribution comes from the last term in eq. (10).

As discussed before, the QED RG scale dependence from the quark mass ratio $z$ in the leading order formula (1) must be reabsorbed by ${\cal O}(\alpha_{\rm em})$ corrections. Indeed the EM LECs $k_{5,7}^{r}$ have the non-trivial UV-scale $\mu$ dependence444This can be derived by computing the operators generated in the chiral Lagrangian by an RG transformation of the quark mass matrix in terms of the EM charge spurions $Q_{L,R}$.:

[TABLE]

It is easy to check that the variation of $k^{r}_{7}$ reabsorbs the dependence induced by the variation of $z$ in the leading order formula (in the $n_{f}=3$ case the RG scale dependence is reabsorbed by $K_{9}^{r}$). In fact, the light quark mass ratio $z$ and the constants $k^{r}_{5,7}$ cannot be determined independently and only the RG invariant combination enters physical quantities. The numerical value of $k_{7}^{r}$ in tab. 1 is of the same order of the scale dependence in eq. (13), which therefore dominates its determination. In any case, the current uncertainties on the quark mass ratio $z$ are still bigger than the effects from the scale dependence in $z$, and therefore bigger than the effects from $k^{r}_{7}$.

To complete the computation of $\chi_{\rm top}$ we need the value of the pion decay constant $f_{\pi^{+}}$ at $\alpha_{\rm em}=0$. Currently the best determination comes from the charged pion leptonic decay, which according to the PDG [28] provides $f_{\pi^{+}}=92.28(9)$. This estimate however involves EM corrections of the same order of $\delta_{e}$, so that a consistent calculation of $\chi_{\rm top}$ within the chiral expansion should consider the two sources of EM corrections together. In more details $f_{\pi^{+}}$ is related to the EM inclusive pion decay rate via

[TABLE]

where the $\delta_{\Gamma}$ terms computed in [34] are the ${\cal O}(\alpha_{\rm em})$ corrections, which we split into two terms: the local contribution $\delta_{\Gamma}^{\rm loc}$ and the IR one $\delta_{\Gamma}^{\rm had}$, which parametrizes the hadronic form factors and depends on the chiral LECs. Explicitly they read:

[TABLE]

The constants $\bar{c}_{i}$ have been estimated in a model dependent way in ref. [34], and for this reason we assign a conservative 100% uncertainty to them. The corresponding numerical values from [34] are reported in tab. 2.

The constants $K_{i}^{r}$ and $X_{i}^{r}$ are the $n_{f}=3$ radiative and leptonic LECs respectively, defined in refs. [36, 37] (except for $\tilde{X}_{6}^{r,{\rm eff}}$ defined in ref. [35]). We use the values estimated in ref. [32] for $K_{1,\dots,6}^{r}$, in ref. [33] for $K_{9,10}^{r}$, and in ref. [35] for $X_{i}^{r}$, which we report in tab. 2 and to which we associate conservatively a 100% uncertainty.

Numerically the size of the EM corrections to $\Gamma_{\pi^{+}\to\mu\nu(\gamma)}$ amounts to

[TABLE]

very close to the PDG estimate (0.0176(21)) but with larger error (here we have been more conservative). Note that, although the uncertainties of all LECs have been taken ${\cal O}(1)$, the result has a 20% accuracy, since the first term in $\delta_{\Gamma}^{\rm loc}$ largely dominates over all the others. Combining eqs. (14) and (16) we get $f_{\pi^{+}}=92.26(18)$.

Since some of the LECs appearing in $\delta_{\Gamma}^{\rm had}$ are common with some of those appearing in $\delta_{e}$, the topological susceptibility $\chi_{\rm top}$ should be written directly in terms of the $\Gamma_{\pi^{+}\to\mu\nu(\gamma)}$ rather than $f_{\pi^{+}}$, i.e.

[TABLE]

The combination $\delta_{e}-\delta_{\Gamma}^{\rm had}$ can be written either all entirely in terms of $n_{f}=3$ LECs, or in a hybrid way in terms of $n_{f}=2$ $k_{i}^{r}$ and $n_{f}=3$ $X_{i}^{r}$ (because the $n_{f}=2$ leptonic LECs are not available):

[TABLE]

In this way we get

[TABLE]

which in combination with eq. (17) can be used to evaluate the EM contribution to $\chi_{\rm top}$.

The direct extraction of $f_{\pi^{+}}$ from lattice QCD simulations is not competitive with the estimate above. However, recently the EM corrections to $\Gamma_{\pi^{+}\to\mu\nu(\gamma)}$ have been computed in a preliminary study on the lattice [38], giving555Note however that ref. [39] alerts about upcoming results which slightly deviate from this quoted value.

[TABLE]

which is in very good agreement with eq. (16). Accidentally this value is very close to the PDG one, both in size and uncertainty. Eq. (20) implies $f_{\pi^{+}}=92.30(7)$. Given the compatibility of the chiral and the lattice results, and the fact that the latter has better precision and less model dependence, we will use eq. (12) and this lattice estimation for $f_{\pi^{+}}$, bearing in mind that numerically this choice is also equivalent to using the PDG determination.

3 NNLO corrections

Given the smallness of the expansion parameter, the $n_{f}=2$ chiral expansion is expected to converge very fast and the NNLO corrections to be only a few percent with respect to the NLO ones. Nevertheless, depending on the magnitude of the low-energy coefficients, they might be competitive with the EM corrections discussed in the previous section. Their estimation is therefore essential for a precise calculation of $\chi_{\rm top}$.

The NNLO corrections to the topological susceptibility receive contributions from the diagrams in fig. 1, which correspond to: 1) the two-loop diagrams constructed from $\mathcal{O}(p^{2})$ vertices, 2) the one-loop diagrams generated by $\mathcal{O}(p^{4})$ vertices and 3) the tree-level graphs from the $\mathcal{O}(p^{6})$ Lagrangian.666Note that even with a quark field redefinition that avoids the tree-level mixing between the axion and the neutral pion, this mixing arises at one loop, producing effects at NNLO in $\chi_{\rm top}$. At $\alpha_{\rm em}=0$, the full NNLO result is

[TABLE]

with $\delta_{1}$ given in eq. (1) and the $c_{i}^{r}$ being the $\mathcal{O}(p^{6})$ $n_{f}=2$ LECs introduced in [40, 41]. Note that the charged and neutral pion decay constants (defined, as mentioned, at $\alpha_{\rm em}=0$) differ only at two loops due to isospin breaking effects, so $f_{\pi}$ in $\delta_{1}$ and $\delta_{2}$ can be understood either as $f_{\pi^{+}}$ or $f_{\pi^{0}}$, being the difference accounted by $\mathcal{O}(p^{8})$ terms. The scale dependence of the combinations of $c_{i}^{r}$ in eq. (22) is fully reabsorbed by the $\log\frac{m_{\pi}^{2}}{\bar{\mu}^{2}}$ and $\log^{2}\frac{m_{\pi}^{2}}{\bar{\mu}^{2}}$ terms, and since $h_{1}^{r}-h_{3}^{r}-\ell_{4}^{r}$, $\ell_{7}$ and $\bar{\ell}_{3,4}$ are scale invariant, the scale dependence of $\delta_{2}$ cancels separately in each line of eq. (22).

While the numerical value of most of the $\mathcal{O}(p^{4})$ LECs is reasonably well known, the determination of the $c_{i}^{r}$ is in much worse shape. In fact, only few combinations of $c_{i}^{r}$ can be extracted directly because there are not enough experimental observables to fit all the $n_{f}=2$ Lagrangian parameters.777In fact, of the $c_{i}^{r}$ that appear in eq. (22), only $c_{6}^{r}$ has been extracted semi-directly from experiments, in particular from the pion scalar form factor [42], with some phenomenological modeling. As explained in app. B, since on its numerical value there is still disagreement [43, 42, 44], we will not use it in our analysis. Recent partially-quenched lattice QCD simulations [27] provided results for some of the combinations of $c_{i}^{r}$ appearing in eq. (22). For the remaining ones we matched the relevant combinations to the $n_{f}=3$ LECs, for which numerical fits exist [45] (taken with a conservative $100\%$ error). In this way we have been able to extract an order of magnitude estimate for all the $c_{i}^{r}$ appearing in eq. (22), which we report in tab. 3 (see app. B for more details).

The LECs in tab. 3 and eq. (5), combined with the values $\bar{\ell}_{3}=2.81(49)$ and $\bar{\ell}_{4}=4.02(25)$ from [27], lead to the following numerical result for the NNLO corrections:

[TABLE]

While the uncertainty from $z$ is very small, those from $\ell_{i}^{r}$ (of which $\ell_{7}$ provides the largest contribution) and $c_{i}^{r}$ have similar size. Notice that although the relative uncertainties of the $c_{i}^{r}$ are large, they only have a milder impact on the final uncertainty of $\delta_{2}$, because numerically $\delta_{2}$ receives bigger contributions from the ${\cal O}(p^{4})$ LECs and the non-local contributions. Moreover, the isospin-breaking terms in $\delta_{2}$ (the last two lines in eq. (22)), which are suppressed by powers of $\Delta^{2}\approx 0.1$, contribute less than 20% to the final result and are within the uncertainty of $\delta_{2}$. As a consequence, the precision on the LECs is still not enough for the result to be sensitive to isospin breaking corrections. Finally notice that $\delta_{2}$ is numerically of the same order of the EM corrections in eq. (12), but with opposite sign. Therefore, both have to be considered for a sub-percent estimate of $\chi_{\rm top}$.

4 Final Results and Axion Mass

We can now combine the analysis of secs. 1, 2 and 3 and estimate of the topological susceptibility to $\mathcal{O}(p^{6},e^{2}p^{2})$. The final result reads

[TABLE]

where the $\mathcal{O}(p^{4})$ contribution $\delta_{1}$ is given in eq. (1), the $\mathcal{O}(p^{6})$ contribution $\delta_{2}$ in eq. (22) and the $\mathcal{O}(e^{2}p^{2})$ contribution $\delta_{e}$ in eq. (10). For completeness, in app. A we also report $\chi_{\rm top}$ expressed in terms of quark masses and bare chiral Lagrangian parameters.

Substituting our numerical estimates, the final results for the topological susceptibility and the axion mass read

[TABLE]

Notice how these values almost coincide with the updated NLO ones in eq. (7), since both NNLO and EM corrections are comparable but smaller than the present uncertainties of the NLO estimate, and, having opposite sign, they tend to cancel each other. This result confirms the reliability of the NLO estimate in [19]. It is instructive to deconstruct the contributions at each order with the various uncertainties: for the axion mass case they read

[TABLE]

where the reported uncertainties on each contributions come from those of $z$, from the EM corrections in the extraction of $f_{\pi^{+}}$, and from those of the various LECs in the NLO ($\ell^{r}_{i}$ and $h_{i}^{r}$), in the NNLO ($c_{i}^{r}$) and in the EM ($k_{i}^{r}$) chiral Lagrangians.

Several comments are in order. First notice how, while NLO corrections are almost two orders of magnitude smaller than the LO result, NNLO are barely one order of magnitude below the NLO ones. On one side this means that the chiral expansion is nicely converging and, given the current uncertainties on $z$ and the LECs, the NLO result is enough. On the other side, the size of the NNLO corrections is such that they cannot be ignored in future improvements of $m_{a}$.

EM corrections are of similar size, slightly less than 0.5% and within present uncertainties. The numerical estimate of the EM corrections has been carried out using the lattice QCD results for $f_{\pi^{+}}$ extracted from eqs. (14) and (20) with $\delta_{e}$ in eq. (12), since these values are more model-independent.888A very similar result would follow using the PDG value for $f_{\pi^{+}}$, which, as mentioned, is very close to the lattice estimate both in size and error. However, one could have also used eqs. (17) and (19) obtaining essentially the same central value although with an error twice as large. As for NNLO corrections, they must be considered should the uncertainties coming from $z$ and the NLO LECs decrease. As commented before, the size of these corrections also represents the ultimate precision that can be reached in lattice estimates which do not include EM corrections.

We conclude by noticing that, if the uncertainties in $m_{u}/m_{d}$ and the NLO LECs (in particular $\ell_{7}$) are reduced by a factor of few (which is not unreasonable) our results could be used to determine the axion mass (and $\chi_{\rm top}^{1/4}$) with per-mille accuracy.

Acknowledgements

We are grateful to Johan Bijnens for correspondence and clarifications on the formulas for the pion mass and the decay constant in refs. [46, 47].

Appendix A Results in terms of the quark masses

We provide here the results for pion mass, pion decay constant and topological susceptibility in both $n_{f}=2$ and $n_{f}=3$ chiral perturbation theory (in this last case in the unbroken isospin limit $m_{u}=m_{d}$) in terms of the bare chiral Lagrangian parameters and quark masses. These are intermediate results used to obtain the formulas in the main sections and the matching of app. B.

A.1 Two-flavor results

We start with the neutral pion mass $m_{\pi^{0}}$, calculated at $\mathcal{O}(p^{6},p^{2}e^{2})$, and the charged pion decay constant $f_{\pi^{+}}$, defined at $\alpha_{\rm em}=0$ and calculated at $\mathcal{O}(p^{6})$. These have been calculated in the unbroken isospin limit in [48, 49, 50, 51, 30], and read:

[TABLE]

where $M^{2}\equiv B(m_{u}+m_{d})$, $B$ and $F$ are the LECs of the leading order $n_{f}=2$ Lagrangian [20], and

[TABLE]

where $\kappa\equiv(4\pi)^{-2}$ and $\lambda_{M}\equiv\log\frac{M^{2}}{\bar{\mu}^{2}}$. The $\mathcal{O}(p^{6})$ contributions are

[TABLE]

At the same order, the topological susceptibility reads

[TABLE]

where

[TABLE]

It is a nontrivial consistency check that the dependence on the scale $\bar{\mu}$ cancels separately in any of the previous equations. Moreover, the QED running of the quark masses is compensated by the shift of $k_{i}^{r}$ as explained in section 2, in such a way that both $m_{\pi^{0}}$ and $\chi_{\rm top}$ are independent of the QED RG scale $\mu$. Inverting eqs. (27-28) for $M$ and $F$ and plugging the result into eq. (34), we obtain the topological susceptibility $\chi_{\rm top}$ expressed as a function of the physical $\pi^{0}$ mass and $f_{\pi^{+}}$ only, as in eq. (24).

A.2 Three-flavor results

In the unbroken isospin limit $m_{u}=m_{d}\equiv m$ and at $\alpha_{\rm em}=0$, the pion mass and decay constant at NNLO in $n_{f}=3$ chiral perturbation theory are

[TABLE]

where $M_{0}^{2}\equiv 2B_{0}m$, $B_{0}$ and $F_{0}$ are the LECs of the leading order $n_{f}=3$ Lagrangian of [52] and

[TABLE]

Here $m_{s}$ is the strange quark mass, $w\equiv m/m_{s}$, $M_{K}^{2}\equiv B_{0}(m+m_{s})$ and $M_{\eta}^{2}\equiv(4/3)M_{K}^{2}-(1/3)M_{0}^{2}$ are the tree-level kaon and eta masses, $\lambda_{P}\equiv\log\frac{M_{P}^{2}}{\bar{\mu}^{2}}$ for $P=0,K,\eta$, and $L_{i}^{r}$ are the NLO LECs of the $n_{f}=3$ chiral Lagrangian of [52]. The terms $\epsilon_{2}^{M}$ and $\epsilon_{2}^{F}$ are $\mathcal{O}(p^{6})$ and depend on the LECs of the $n_{f}=3$ NNLO Lagrangian, $C_{i}^{r}$. Both $\epsilon_{2}^{M}$ and $\epsilon_{2}^{F}$ involve the calculation of a two-mass scale sunset integral at non-zero external momentum, and therefore do not admit a closed analytic expression (see [47], where a two-integral representation is provided). However, by employing the analytic form of the two-mass scale sunset integral at vanishing external momentum [53, 47] and the recursion relations for sunset integrals [56, 55, 54], $\epsilon_{2}^{M}$ and $\epsilon_{2}^{F}$ can be expanded in power series of $m$. Such an expansion has been calculated for the first time in [54] for $\epsilon_{2}^{M}$ and [46] for $\epsilon_{2}^{F}$ up to $\mathcal{O}(m^{3})$ and $\mathcal{O}(m^{2})$ respectively and is sufficient for the matching of the $n_{f}=2$ and $n_{f}=3$ LECs discussed in app. B. Since the result of $\epsilon_{2}^{M}$ and $\epsilon_{2}^{F}$ turns out to be quite involved, we will avoid reporting here the explicit expressions, for which we refer to [54, 46].

Finally, the topological susceptibility at $\alpha_{\rm em}=0$ reads

[TABLE]

where in the unbroken isospin limit $m_{u}=m_{d}\equiv m$:

[TABLE]

The result of $\epsilon_{2}^{\chi}$ has been conveniently organized into different contributions: $\epsilon_{2,C}^{\chi}$ are terms containing the LECs $C_{i}^{r}$, $\epsilon_{2,\log}^{\chi}$ and $\epsilon_{2,\log\times\log}^{\chi}$ are terms respectively linear and quadratic in the chiral logs $\lambda_{0,K,\eta}$ (but without $L_{i}^{r}$), $\epsilon_{2,L}^{\chi}$ and $\epsilon_{2,L\times L}^{\chi}$ are terms linear and quadratic in $L_{i}^{r}$ (but independent of chiral logs), and $\epsilon_{2,\log\times L}^{\chi}$ contains products of chiral logs and LECs. Finally, $\epsilon_{2,{\rm finite}}^{\chi}$ is the remaining constant piece and is automatically scale independent. In particular:

[TABLE]

In this last equation we defined

[TABLE]

which arises in the evaluation of a two-mass scale sunset integral at vanishing external momentum [46].

Appendix B Extraction of the NNLO LECs and input parameters

As mentioned in sec. 3, the $\mathcal{O}(p^{6})$ LECs $c_{i}^{r}$ of the $n_{f}=2$ chiral Lagrangian [40, 41] are still poorly known from experimental data (see [45] for a review). To estimate the value of the $c_{i}^{r}$’s appearing in eq. (22), we combined the results from recent lattice QCD simulations [27] with the information from the matching of pion mass, decay constant and topological susceptibility for $n_{f}=2$ and $n_{f}=3$ chiral perturbation theory of app. A and the fits for the $n_{f}=3$ LECs $C_{i}^{r}$ provided in [45].

• •

The SU(2) partially quenched simulation of ref. [27] provide fits of $8$ combinations of SU(2) LECs. In this analysis we consider the $450$ MeV cut-fit for such combinations, reported in the last column of tab. 6 of [27]. While this estimate is less conservative than the one extracted from the $370$ MeV cut-fit, the two are compatible for all reported combinations of LECs. In tab. 4 we quote the results expressed in terms of $c_{i}^{r}$ using the relations between SU$(N)$ and SU($2$) LECs of [40].

• •

Ref. [45] provides fits of $34$ combinations of $\mathcal{O}(p^{6})$ three-flavor LECs, $C_{i}^{r}$. The information of $C_{i}^{r}$ can be translated into a value for three combinations of $c_{i}^{r}$ by equating in the large $m_{s}/m$ limit (and for $\alpha_{\rm em}=0$ and $m\equiv m_{u}=m_{d}$) the $n_{f}=2$ and $n_{f}=3$ formulas for the pion mass in eqs. (27) and (38), the pion decay constant in eqs. (28) and (39) and the topological susceptibility in eqs. (34) and (42). Such a matching leads respectively to the three following relations:

[TABLE]

The notation in eqs. (53-55) is as in app. A.2 (except that $m/m_{\pi}^{2}$, $f_{\pi}$ and $M^{2}_{K,\eta}$ in $\lambda_{K,\eta}$ are computed at $m_{u}=m_{d}=m=0$), and the contributions on the r.h.s. have been ordered as in eq. (44). The numerical value of the NLO couplings $L_{i}^{r}$ is reported in tab. 5: in particular, $L_{4}^{r},L_{5}^{r}$ and $L_{6}^{r}$ are taken from lattice QCD studies [26], while the others from the fit of [45]. Moreover, to all $L_{i}^{r}$ a 30% intrinsic uncertainty from higher order 3-flavor corrections has been added (this is not present for 2-flavor where higher order corrections are much smaller). The value of the NNLO couplings $C_{i}^{r}$ appearing in the r.h.s. of eqs. (53-55) taken from tab. 4 of ref. [45] is also reported in tab. 5. Since the fit in ref. [45] did not provide uncertainties for the $C_{i}^{r}$ coefficients we assume that they reproduce at least the right orders of magnitude and conservatively assign to them a $100\%$ uncertainty.

Eqs. (53-55) then lead to:

[TABLE]

The final value of the $9$ couplings $c_{i}^{r}$ in tab. 3 has been extracted by combining the lattice results in tab. 4 with the 2-3 flavor matching result in eq. (• ‣ B) through a $\chi^{2}$ fit, whose quality ($\chi^{2}\sim 3$) turns out to be good. In principle, an estimate of $c_{6}^{r}$ could be directly extracted from the pion scalar form factor, as in ref. [45] where $c_{6}^{r}\approx-1.9\times 10^{-5}$. However, since there is still a factor $\sim 3$ uncertainty on how to theoretically model this last quantity [43, 42, 44], we chose not to use this estimate of $c_{6}^{r}$ in our numerical analysis. In any case, the NNLO corrections to $\chi_{\rm top}$ in eq. (22) taking into account also $c_{6}^{r}=-1.9(1.9)\times 10^{-5}$ result in $\delta_{2}=-0.006(3)$, still compatible with eq. (23), but with an overall lower quality fit of the $c_{i}^{r}$.

Finally, for convenience in tab. 6 we summarize the values of the parameters used in this work, which should be considered together with the LECs in tabs. 1, 2, 3 and 5. When uncertainties are not quoted it means that their effect was negligible and they have not been used.