Topological susceptibility and $\eta'$ meson mass from $N_f=2$ lattice QCD at the physical point

TL;DR

This study calculates the topological susceptibility and $ ext{η'}$ meson mass in $N_f=2$ lattice QCD at the physical pion mass, demonstrating consistent results and a linear pion mass dependence of the susceptibility.

Contribution

It provides a lattice QCD analysis of topological susceptibility and $ ext{η'}$ mass at the physical point with new insights into their pion mass dependence.

Findings

$ ext{η'}$ mass is approximately 772 MeV with small uncertainty.

Topological susceptibility shows linear dependence on $M_ ext{π}^2$ with lattice artifacts.

Physical point can be reached without increased statistical noise.

Abstract

In this paper we explore the computation of topological susceptibility and $\eta'$ meson mass in $N_f=2$ flavor QCD using lattice techniques with physical value of the pion mass as well as larger pion mass values. We observe that the physical point can be reached without a significant increase in the statistical noise. The mass of the $\eta'$ meson can be obtained from both fermionic two point functions and topological charge density correlation functions, giving compatible results. With the pion mass dependence of the $\eta'$ mass being flat we arrive at $M_{\eta'}= 772(18)\ \mathrm{MeV}$ without an explicit continuum limit. For the topological susceptibility we observe a linear dependence on $M_\pi^2$, however, with an additional constant stemming from lattice artifacts.

Tables

Table 1. Table 1: The gauge ensembles used in this study. The labeling of the ensembles follows the notations in Ref.  Abdel-Rehim:2015pwa . In addition to the relevant input parameters we give the lattice volume ( L / a ) 3 × T / a superscript 𝐿 𝑎 3 𝑇 𝑎 (L/a)^{3}\times T/a and the number of evaluated configurations N conf subscript 𝑁 conf N_{\mathrm{conf}} .

Ensemble	$\beta$	$c_{\mathrm{sw}}$	$a{\mu }_{\ell }$	${\left(L/a\right)}^3\times T/a$	$N_{\mathrm{conf}}$
$cA\mathrm{2.09.48}$	2.10	1.57551	0.009	${48}^3\times 96$	$615$
$cA\mathrm{2.30.48}$	2.10	1.57551	0.030	${48}^3\times 96$	$352$
$cA\mathrm{2.30.24}$	2.10	1.57551	0.030	${24}^3\times 48$	$352$
$cA\mathrm{2.60.32}$	2.10	1.57551	0.060	${32}^3\times 64$	$337$
$cA\mathrm{2.60.24}$	2.10	1.57551	0.060	${24}^3\times 48$	$424$

Table 2. Table 2: Results for the masses of the charged and the neutral pion (full and quark-connected only), the η 2 subscript 𝜂 2 \eta_{2} meson in lattice units (from fermionic and gluonic correlators), the gluonic gradient flow lattice scale t 1 / a 2 subscript 𝑡 1 superscript 𝑎 2 t_{1}/a^{2} , and the topological susceptibility in units of t 1 subscript 𝑡 1 t_{1} for the five ensembles considered.

Ensemble	$a{M}_{\pi }$	$a{M}_{{\pi }^0}$	$a{M}_{{\pi }_c^0}$	$a{M}_{{\eta }_2}^{\text{ferm.}}$	$a{M}_{{\eta }_2}^{\text{gl.}}$	$t_1/a^2$	${10}^3\cdot t_1^2{\chi }_{\text{top}}$
$cA\mathrm{2.09.48}$	$0.06211\left(06\right)$	$0.0576\left(25\right)$	$0.1196\left(02\right)$	$0.361\left(14\right)$	$0.369\left(10\right)$	$6.890\left(08\right)$	$0.48\left(3\right)$
$cA\mathrm{2.30.48}$	$0.11199\left(06\right)$	$0.0976\left(35\right)$	$0.1521\left(01\right)$	$0.376\left(11\right)$	$0.356\left(17\right)$	$6.761\left(08\right)$	$0.56\left(5\right)$
$cA\mathrm{2.30.24}$	$0.11461\left(37\right)$	$0.1110\left(95\right)$	$0.1519\left(29\right)$	$0.425\left(22\right)$	$0.386\left(45\right)$	$6.828\left(30\right)$	$0.46\left(4\right)$
$cA\mathrm{2.60.32}$	$0.15783\left(12\right)$	$0.1555\left(64\right)$	$0.1883\left(03\right)$	$0.396\left(10\right)$	$0.379\left(14\right)$	$6.562\left(08\right)$	$0.78\left(4\right)$
$cA\mathrm{2.60.24}$	$0.15908\left(28\right)$	$0.1347\left(86\right)$	$0.1883\left(11\right)$	$0.399\left(12\right)$	$0.345\left(48\right)$	$6.550\left(17\right)$	$0.79\left(7\right)$

Equations

$$D_{\ell }=D-i{\gamma }_5{\tau }_3\left[W_{\mathrm{cr}}+\frac{i}{4}{c}_{\mathrm{sw}}{\sigma }^{\mu \nu }{\mathcal{F}}^{\mu \nu }\right]+{\mu }_{\ell },$$

$${\mathcal{P}}^3(t)=\frac{1}{\sqrt{2}}\sum _{\mathbf{x}}\bar{\psi }i{\gamma }_5{\tau }^3\psi (\mathbf{x},t),{\mathcal{P}}^0(t)=\frac{1}{\sqrt{2}}\sum _{\mathbf{x}}\bar{\psi }i{\gamma }_5\text{mathds}{1}_f\psi (\mathbf{x},t).$$

$$\displaystyle C_{\pi^{0}}(t-t^{\prime})\$$

$$\displaystyle C_{\eta_{2}}(t-t^{\prime})\$$

$$\begin{array}{rl}⟨{\bar{\psi }}_{u}i{\gamma }_5{\psi }_u(x)+{\bar{\psi }}_{d}i{\gamma }_5{\psi }_d(x){⟩}_F& =-i\mathrm{Tr}\{{\gamma }_5{G}_u^{xx}\}-i\mathrm{Tr}\{{\gamma }_5{G}_d^{xx}\}\\ & =-i\mathrm{Tr}\{{\gamma }_5{G}_u^{xx}\}-i\mathrm{Tr}\{(G_u^{xx}{)}^{†}{\gamma }_5\}\\ & =-2i\mathrm{Re}\mathrm{Tr}\{{\gamma }_5{G}_u^{xx}\}.\end{array}$$

$$⟨{\bar{\psi }}_{u}i{\gamma }_5{\psi }_u(x)-{\bar{\psi }}_{d}i{\gamma }_5{\psi }_d(x){⟩}_F=2\mathrm{Im}\mathrm{Tr}\{{\gamma }_5{G}_u^{xx}\}.$$

$$C^{\mathrm{conn}}(t-t^{\prime })=\mathrm{Re}\mathrm{Tr}\{{\gamma }_5{G}_u^{t{t}^{\prime }}{\gamma }_5{G}_u^{t^{\prime }t}\},$$

$$\begin{array}{rl}{C}_{{\pi }^0}(t-t^{\prime })& =\mathrm{Tr}\{{\gamma }_5{G}_u^{t{t}^{\prime }}{\gamma }_5{G}_u^{t^{\prime }t}\}+2\mathrm{Im}\mathrm{Tr}\{{\gamma }_5{G}_u^{tt}\}\cdot \mathrm{Im}\mathrm{Tr}\{{\gamma }_5{G}_u^{t^{\prime }{t}^{\prime }}\},\\ C_{{\eta }_2}(t-t^{\prime })& =\mathrm{Tr}\{{\gamma }_5{G}_u^{t{t}^{\prime }}{\gamma }_5{G}_u^{t^{\prime }t}\}-2\mathrm{Re}\mathrm{Tr}\{{\gamma }_5{G}_u^{tt}\}\cdot \mathrm{Re}\mathrm{Tr}\{{\gamma }_5{G}_u^{t^{\prime }{t}^{\prime }}\}.\end{array}$$

$$C_{{\pi }^{\pm }}(t-t^{\prime })=⟨{\mathcal{P}}^{+}(t)({\mathcal{P}}^{+}(t^{\prime }){)}^{†}⟩$$

$${\mathcal{P}}^{+}(t)=\sum _{\mathbf{x}}\bar{\psi }i{\gamma }_5\frac{{\tau }^1+i{\tau }^2}{2}\psi (\mathbf{x},t)$$

$$q(x)=-\frac{1}{32{\pi }^2}{\epsilon }_{\mu \nu \rho \sigma }\mathrm{Tr}{F}_{\mu \nu }(x)F_{\rho \sigma }(x)$$

$$C_{qq}(x-y)=⟨q(x)q(y)⟩.$$

$${\chi }_{\text{top}}=\frac{1}{V}\int \mathrm{d}x\int \mathrm{d}y⟨q(x)q(y)⟩$$

$$C_{qq}(x,y)\sim \frac{M}{|x-y|}{K}_1(M\cdot |x-y|)$$

$$Q(t)=a^4\sum _x{q}_t(x),$$

$$M_{\mathrm{eff}}=-\log \frac{C_{{\eta }_2}(t)}{C_{{\eta }_2}(t+1)}$$

$$C(t)=|⟨0|\mathcal{O}|0⟩{|}^2+\sum _n\frac{|⟨0|\mathcal{O}|n⟩{|}^2}{2E_n}\left(e^{-E_{n}t}+e^{-E_n(T-t)}\right),$$

$$\tilde{C}(t)=C(t)-C(t+1).$$

$$\tilde{C}(t)\propto \left(e^{-E_{n}t}-e^{-E_n(T-t)}\right).$$

$$\propto \frac{a^5}{T}\left({\chi }_{\text{top}}+\frac{Q^2}{V}\right)$$

$$f^{\pm }(t;A,M)=A\left(e^{-Mt}\pm e^{-M(T-t)}\right),$$

$$C_{qq}(r)\sim \sqrt{\frac{M}{r}}\frac{1}{r}{e}^{-Mr}\left(1+\mathcal{O}\left(\frac{1}{Mr}\right)\right)\text{for large r,}$$

$$a{M}_{{\eta }_2}^{\text{gl.}}=0.3687(98)$$

$$t_1^2{\chi }_{\text{top}}=\frac{1}{8}{t}_1^2{f}_{\pi }^2{M}_{{\pi }^0}^2+a^2\frac{c_2}{t_1}.$$

$$g(M_{\pi }^2)=c_1{t}_1{M}_{\pi }^2+a^2\frac{c_2}{t_1}$$

$$r_0{M}_{{\eta }_2}=1.92(8)$$

$$M_{{\eta }_2}=772(18)\mathrm{MeV}$$

Full text

I Introduction

Due to the persisting $3-5\,\sigma$ deviation in the anomalous magnetic moment of the muon $a_{\mu}$ between theory and experiment there is considerable interest in the decays $\eta\to\gamma^{\star}\gamma^{\star}$ and $\eta^{\prime}\to\gamma^{\star}\gamma^{\star}$ because a better knowledge of the corresponding transition form factors could help to reduce the uncertainty in the hadronic light-by-light contribution to $a_{\mu}$; see for instance Ref. Jegerlehner:2015stw . Moreover, $\eta$ and $\eta^{\prime}$ mesons are interesting from a theoretical point of view because the large mass of the $\eta^{\prime}$ meson is explained by the anomalously broken $U_{A}(1)$ axial symmetry in QCD. The $\eta,\eta^{\prime}$ mixing pattern and the aforementioned transition form factors can be computed nonperturbatively using lattice techniques.

There has been considerable progress in studying $\eta$ and $\eta^{\prime}$ mesons from lattice QCD. In Refs. Michael:2013gka ; Ottnad:2017bjt the corresponding mixing has been studied for three values of the lattice spacing and a large, but still unphysical range of pion mass values in $N_{f}=2+1+1$ flavor QCD. After extrapolation to the physical pion mass value excellent agreement to experiment was found. Further lattice results for $\eta,\eta^{\prime}$ can be found in Refs. Christ:2010dd ; Dudek:2011tt ; Gregory:2011sg ; Dudek:2013yja ; Fukaya:2015ara .

Through the anomaly, the mass of the $\eta^{\prime}$ is also tightly connected to topology and in particular the topological susceptibility $\chi_{\text{top}}$. The latter quantity must decrease as $M_{\pi}^{2}$ toward the chiral limit, if the $\eta^{\prime}$ is not a Goldstone boson Leutwyler:1992yt . For recent lattice studies of the topological susceptibility; see for instance Aoki:2017paw ; Alexandrou:2017bzk . There is now particular interest in $\chi_{\text{top}}$ due to its connection to axion dark matter; see for instance Refs. Dowrick:1991sj ; Moore:2017ond ; DiVecchia:2017xpu .

In this paper we attempt to study the $\eta^{\prime}$ meson and the topological susceptibility directly at the physical point, however, in a first step in $N_{f}=2$ flavor QCD. In $N_{f}=2$ flavor QCD there exist a pion triplet and one flavor singlet, which is related to the aforementioned anomaly. We will denote it as the $\eta_{2}$ meson to distinguish it from the $\eta^{\prime}$ meson in full QCD, which is only approximately a flavor eigenstate. The $\eta_{2}$ and the $\eta^{\prime}$ meson have in common that both receive significant fermionic disconnected contributions. In an earlier study McNeile:2000hf their masses have been found to differ only by $200\ \mathrm{MeV}$, with the additional strange quark introducing only a moderate shift in the mass. In particular, both are expected to have a similar dependence on the light quark mass. The most recent lattice QCD studies of the $\eta_{2}$ meson can be found in Refs. Sun:2017ipk ; Jansen:2008wv .

We investigate $M_{\eta_{2}}$ using fermionic correlation functions and in addition topological charge density correlators. The topological susceptibility is studied using gradient flow techniques Luscher:2010iy . Studying the $\eta_{2}$ meson and the topological susceptibility at the physical point will reveal on the one hand important qualitative information on the implementation of the anomaly in QCD. On the other hand it represents a feasibility study for a later investigation of $\eta$ and $\eta^{\prime}$ in $N_{f}=2+1+1$ QCD at the physical point Alexandrou:2018egz . The results obtained here are also important prerequisites for an exploratory study of $\eta_{2}\to\gamma^{\star}\gamma^{\star}$.

The paper is organized as follows: in the following two sections we discuss the lattice details of our computation. In Sec. IV we present the analysis methods and in Sec. V the results. We close with a discussion and summary. For a first account of this work we refer to Ref. Helmes:2017ccf .

II Lattice action

The results presented in this paper are based on the gauge configurations generated by the ETMC with Wilson clover twisted mass quark action at maximal twist Abdel-Rehim:2015pwa . We employ the Iwasaki gauge action Iwasaki:1985we . The measurements are performed on a set of $N_{f}=2$ ensembles with the pion mass ranging from its physical value to 340 MeV. In Table 1 we list all the ensembles together with the relevant input parameters, the lattice volume, and the number of configurations. The lattice spacing is $a=0.0931(2)\ \mathrm{fm}$ for all five ensembles. More details about the ensembles are presented in Ref. Abdel-Rehim:2015pwa .

The sea quarks are described by the Wilson clover twisted mass action. The Dirac operator for the light quark doublet consists of the Wilson twisted mass Dirac operator Frezzotti:2000nk combined with the clover term

[TABLE]

which acts on a flavor doublet spinor $\psi=(u,d)^{T}$. In Eq. (1) we have $D=\gamma_{\mu}(\nabla^{\ast}_{\mu}+\nabla_{\mu})/2$ with $\nabla_{\mu}$ and $\nabla_{\mu}^{\ast}$ the forward and backward lattice covariant derivatives, and the Wilson term $W_{\mathrm{cr}}=-ra\nabla_{\mu}^{\ast}\nabla_{\mu}+m_{\mathrm{cr}}$ with the critical mass $m_{\mathrm{cr}}$, the Wilson parameter $r=1$, and the lattice spacing $a$. The average up/down (twisted) quark mass is denoted by $\mu_{\ell}$, while $c_{\mathrm{sw}}$ is the so-called Sheikoleslami-Wohlert improvement coefficient Sheikholeslami:1985ij multiplying the clover term. It is in our case not used for $\mathcal{O}(a)$ improvement but serves to significantly reduce the effects of isospin breaking Abdel-Rehim:2015pwa .

The critical mass has been determined as described in Refs. Chiarappa:2006ae ; Baron:2010bv . This guarantees automatic $\mathcal{O}\left(a\right)$ improvement Frezzotti:2003ni , which is one of the main advantages of the Wilson twisted mass formulation of lattice QCD.

III Observables

As a smearing scheme in the computation of fermionic correlation functions we use the stochastic Laplacian Heaviside (sLapH) method Peardon:2009gh ; Morningstar:2011ka . The details of our sLapH parameter choices for a set of $N_{f}=2+1+1$ Wilson twisted mass ensembles are given in Ref. Helmes:2015gla . The parameters for the ensembles used in this work are the same as those for $N_{f}=2+1+1$ ensembles with the corresponding lattice volume.

III.1 $\eta_{2}$ and pion correlation functions

In $N_{f}=2$ flavor QCD there is the neutral pion, corresponding to the neutral of the three pions in the triplet, and the $\eta_{2}$, the flavor singlet pseudoscalar meson related to the axial $U_{A}(1)$ anomaly. Since up and down quarks are mass degenerate, there is no mixing among the neutral pion and the $\eta_{2}$ with our action. We employ the following pseudoscalar interpolating operators projected to zero momentum, which are all local and Hermitian

[TABLE]

Here, $\tau^{3}$ is the third Pauli and $\mathds{1}_{f}$ the unit matrix, both acting in flavor space. From those one builds the correlation functions

[TABLE]

which allow one to determine the masses $M_{\pi^{0}}$ and $M_{\eta_{2}}$ from their decay in Euclidean time. Both correlation functions in Eqs. (3) and (4) do have a fermionic connected and a fermionic disconnected contribution, the latter of which vanishes exactly in case of the neutral pion in an isospin symmetric theory. Since this is not the case for Wilson twisted mass fermions, we have to take the disconnected contributions into account also for the $\pi^{0}$.

For the disconnected part of $C_{\eta_{2}}$ we consider the loop

[TABLE]

Here, we have used the $\gamma_{5}$ hermiticity property $D_{d}=\gamma_{5}D_{u}^{\dagger}\gamma_{5}$. $G_{u/d}^{xy}$ represents the up or down propagator. Similarly, one shows for $C_{\pi^{0}}$

[TABLE]

The fermionic connected contribution is identical for the two correlation functions Eq. (3) and Eq. (4). With similar arguments as for the loops one finds

[TABLE]

where we have suppressed the spatial indices. From Eqs. (5), (6), and (7) we infer the expressions for the $\pi^{0}$ and $\eta_{2}$ correlation functions as follows:

[TABLE]

For completeness, the correlation function of the charged pion is constructed as

[TABLE]

with

[TABLE]

and $\tau^{1}$ and $\tau^{2}$ the first and second Pauli matrices, respectively.

III.2 Topological charge density correlations and

susceptibility

The naive field theoretical definition of the topological charge density given by

[TABLE]

defines the topological charge density point-to-point correlator

[TABLE]

The topological susceptibility, which is a measure for the fluctuations of the topological charge, is defined as

[TABLE]

where $V$ is the spacetime volume. Due to the pseudoscalar nature of the topological charge density the topological charge density correlator is strictly negative for finite separations, $C_{qq}(x-y>0)<0$. On the other hand, it is clear that the susceptibility is strictly positive, because $V\cdot\chi_{\text{top}}=\langle Q^{2}\rangle>0$, where $Q=\int dx\,q(x)$ is the total topological charge of the gauge field. This apparent contradiction is resolved by recalling that $C_{qq}$ suffers from contact-term singularities at $x-y=0$ which need to be renormalized in order for the susceptibility to make physical sense. Hence, the physics of the topological susceptibility is intricately hidden in the difference between the contact-term contribution of the correlator at $|x-y|=0$ and the contributions at $|x-y|>0$.

Another interesting property of the topological charge density correlator is that it couples to the flavor singlet pseudoscalar mesons. It is in fact this coupling which is thought to be responsible for the large mass of the $\eta_{2}$. As a consequence, the behavior of the topological charge density point-to-point correlator is dominated by the single boson propagator for the $\eta_{2}$ meson and therefore follows the form of the scalar propagator Dowrick:1991sj ; Shuryak:1994rr

[TABLE]

where $K_{1}$ is the modified Bessel function of the first kind, and $M$ is the mass of the lightest particle in the pseudoscalar meson sector, i.e., the mass $M_{\eta_{2}}$ of the $\eta_{2}$ meson.

On the lattice the topological charge density is discretized with a clover-type discretization of the field strength tensor $F_{\mu\nu}$ which extends over a distance of $2a$ in lattice units. Hence the contact-term contributions to the topological charge density correlator are also distributed over the distance $|x-y|\sim 2a$. Moreover, since the discretized field strength tensor is constructed from smeared gauge links in order to remove ultraviolet fluctuations of the gauge field, the positive contact-term contributions to the correlator are spread over a range $R_{0}$ which depends on the details of the smoothing scheme. However, the behavior of the correlator for $|x-y|\gg R_{0}$ should be independent of the details of the smearing scheme and hence any smoothing scheme is supposed to yield the same physics, i.e. the same mass $M_{\eta_{2}}$.

For the topological charge density correlator we use the array processor experiment (APE) smearing scheme Albanese:1987ds with various smearing levels ranging up to 90 iterative smearing steps. This is in order to check for the independence of the results from the smearing scheme. To compute the topological susceptibility $\chi_{\text{top}}$ we employ the gradient flow technique as introduced for lattice QCD in Ref. Luscher:2010iy . It has the advantage of yielding a renormalized topological susceptibility at finite flow time $t$ Luscher:2011bx , in particular it renormalizes the contact term singularities in the continuum limit at any fixed, physical value of $t$. Since the renormalized susceptibility is scale invariant, i.e., independent of the renormalization scale, $\chi_{\text{top}}$ becomes independent of the flow time $t$ at sufficiently large $t$ toward the continuum limit. This is indeed what we observe in our calculation. However, we note that lattice artifacts might well be very different for the susceptibility at different values of $t$. Instead of calculating the topological susceptibility via the lattice version of Eq. (13), we first obtain the topological charge at flow time $t$ from the topological charge density $q_{t}(x)$ evaluated on the flown gauge field configuration,

[TABLE]

and then the susceptibility via $\chi_{\text{top}}(t)=\langle Q(t)^{2}\rangle/V$. We choose $t=3t_{0}$ where $t_{0}$ is the usual gradient flow reference scale defined through the renormalized action density Luscher:2010iy on the corresponding ensemble. In addition, we also make use of the related reference scale $t_{1}$ in order to facilitate the comparison of our results with those in Ref. Bruno:2014ova .

IV Analysis Method

In the following Secs. IV.1 and IV.2 we will first give more details on the analysis of the fermionic correlation function, before we turn to the discussion for the gluonic correlators in Sec. IV.3.

The fermionic correlation function data are generally analyzed using the blocked bootstrap procedure with $R=10000$ bootstrap samples. Depending on the ensemble, we have chosen the block size such that at least $\gtrsim 100$ blocked data points are left. The relevant masses are computed from correlated fits to the correlation function data.

The gluonic correlation function data are analyzed using a jackknife procedure. It turns out that the correlators at separate distances are highly correlated even for large separations, such that the covariance matrix cannot be taken into account reliably in the fitting procedur;e see further details below. The data for the topological charge susceptibility (and the gluonic scales $t_{0}$ and $t_{1}$) are analyzed using a blocked bootstrap procedure with $R=1000$ bootstrap samples and block sizes such that $\gtrsim 30$ blocked data points are left. The so obtained error is compared to the naive one corrected by the integrated autocorrelation time $\tau_{\text{int}}$, and the larger of the two is always chosen as the final error. Since these calculations are inexpensive, we use at least double the number of configurations indicated in Table 1.

IV.1 Excited state subtraction

In particular for the $\eta_{2}$ meson, the fermionic disconnected contributions are very noisy. As a consequence, the signal is lost relatively early in Euclidean time. For this reason we have in the past applied a method to subtract excited states Jansen:2008wv ; Michael:2013gka ; Liu:2016cba ; Ottnad:2017bjt , originally proposed in Ref.  Neff:2001zr . It actually works very well and we will apply it here again for the $\eta_{2}$ meson. It consists of subtracting excited states from the connected contribution only. This is feasible, because the connected part — representing a pion correlation function — has a signal for all Euclidean time values. Therefore, we can fit to it at large enough Euclidean times such that excited states have decayed sufficiently. Next, we replace the connected correlation function at small times by the fitted (ground state) function. Thereafter, the so subtracted connected contribution is summed according to Eq. (8) to the full $\eta_{2}$ correlation functions.

The underlying assumption is that disconnected contributions are large for the ground state, i.e. the $\eta_{2}$, but not for excited states. If this assumption is correct, the effective mass

[TABLE]

should show a plateau from very early Euclidean times on. We have found in Refs. Jansen:2008wv ; Michael:2013gka ; Ottnad:2017bjt that this approach works very well for the $\eta_{2}$ meson in $N_{f}=2$ flavor QCD as well as for $\eta$ and $\eta^{\prime}$ mesons in $N_{f}=2+1+1$ flavor QCD.

IV.2 Shifted correlation functions

The expected time dependence of the fermionic correlation functions considered here reads as follows:

[TABLE]

where $\mathcal{O}=\mathcal{P}^{+},\mathcal{P}^{3},\mathcal{P}^{0}$ and $n$ labels the states with the corresponding quantum numbers. The time independent first term on the right-hand side corresponds to the vacuum expectation value (VEV). Using the symmetries of our action one can show that for $\mathcal{P}^{+}$ and $\mathcal{P}^{0}$ the VEV must be zero, while this is not the case for $\mathcal{P}^{3}$. We deal with the VEV by building the shifted correlation function

[TABLE]

The difference cancels the constant VEV contribution, while also changing the time dependence to be antisymmetric in time,

[TABLE]

As an alternative, one can also compute the VEV $|\langle 0|\mathcal{O}|0\rangle|^{2}$ from the data and subtract it explicitly.

Since the VEV has to be zero for $\mathcal{P}^{0}$ up to statistical fluctuations, strictly speaking we do not need to use the shifting procedure for the $\eta_{2}$ meson. However, as has been argued in Ref. Aoki:2007ka and first investigated in Ref. Bali:2014pva , there is an additive finite volume effect to $C_{\eta_{2}}$ constant in Euclidean time of the form

[TABLE]

proportional to the topological susceptibility $\chi_{\text{top}}$ and the squared topological charge $Q^{2}$. If present, such a term will cause the $\eta_{2}$ correlation function to stay finite at large Euclidean times. Depending on the sign of the coefficient in front of the finite volume effect, the correlation function may even turn negative at relatively small Euclidean times. Clearly, a finite volume effect of this type can be subtracted again using the shifting procedure, which has first been proposed and applied in Ref. Ottnad:2017bjt .

IV.3 Topological charge density correlators

For the computation of the topological charge density correlator we make use of the full translational invariance. In order to do so, we obtain the topological charge density correlator in Eq. (12) by Fourier transforming the topological charge density on each gauge field configuration, calculating the correlator in Fourier space and transforming it back to coordinate space. In this sense the evaluation is exact, in contrast to the computation of the disconnected contributions to the fermionic correlators in Eqs. (3) and (4), which can only be evaluated stochastically.

The employed smearing level has several effects on the correlation function $C_{qq}(x-y)$. First, it reduces the statistical errors because the smearing suppresses ultraviolet fluctuations. Hence, with increasing smearing levels the signal can be followed over larger and larger separations $x-y$. Second, the increased smearing enhances the contribution of the ground state in the correlation function, i.e. in this case the contribution of the $\eta_{2}$ state. Third, with an increasing smearing level the contact term is distributed over larger distances and hence distorts the correlation function up to larger and larger separations. Obviously, these effects compete with each other with respect to the optimal fit range.

In principle, the choice of the fit ranges should be determined by the quality of the fits. Unfortunately, here this is not possible, because the correlators at separate distances $r$ and $r^{\prime}$ are highly correlated. We illustrate this in Figure 1 where we show the covariance Cov($C_{qq}(r),C_{qq}(r^{\prime})$) of the correlation functions as a function of $(r-r^{\prime})/a$ for different values of $r$ and smearing level $n=90$.111The data for the other smearing levels look very similar.

We are essentially looking at separate columns of the covariance matrix. For ease of comparison we normalize the covariance by Cov($C_{qq}(r),C_{qq}(r)$), and in addition bin the data into bins of size $\Delta r/a=0.125$.

We see that the $C_{qq}$’s are positively correlated for $(r^{\prime}-r)/a\lesssim 4.5$ and become more and more strongly anticorrelated until a maximum of anticorrelation is reached at around $(r^{\prime}-r)/a\sim 7.5$. Since this correlation is essentially independent of $r$, the columns of the covariance matrix are highly linear dependent and the matrix itself is very ill-conditioned. As a consequence, it cannot be taken into account for reliably estimating the quality of the $\chi^{2}$-fits.

We note that all the above conclusions hold independently of the bin size and the smearing level, and we suspect that the peculiar behavior is due to some underlying structure in the topology of the gauge fields.

V Results

In Table 2 we show results for the pion and $\eta_{2}^{\text{ferm.}}$ masses which have been computed from fermionic correlation functions, the $\eta_{2}^{\text{gl.}}$ masses obtained from the gluonic topological charge density correlation functions, as well as the gluonic gradient flow lattice scale $t_{1}/a^{2}$ and the topological susceptibility $\chi_{\text{top}}$. For the charged pion we will always use the shorthand $M_{\pi}$, while for the neutral pion $M_{\pi^{0}}$ and $M_{\pi^{0}_{c}}$ refer to the full and quark-connected masses, respectively. In the following we discuss these results in more detail.

V.1 Neutral pion

In contrast to the $\eta_{2}$ meson discussed later, the signal for the neutral pion can be resolved for all values of $t/a$. In the left panel of Figure 2 we show the shifted correlation function $\tilde{C}(t)$ for the neutral pion as well as the individual quark-connected and quark-disconnected contributions. Note that in this case the function shift is required to remove the offset from the vacuum expectation value in the quark-disconnected contribution. Clearly the signal in the quark-disconnected part is well behaved even for the largest values of $t/a$. For comparison the charged pion has been included in the plot as well.

As already visible from the left panel of Figure 2, charged and neutral pions appear to have very similar mass values. This is even more apparent from the right panel where the neutral pion mass values are plotted versus the charged pion mass values, both in physical units, for all the ensembles considered here. The points fall almost on the bisecting line, which indicates no mass splitting between neutral and charged pion mass. This finding, which we pointed out already in Ref. Abdel-Rehim:2015pwa , is rather important: this mass splitting is basically the only large $a^{2}$ lattice artefact that was found for simulations with Wilson twisted fermions at maximal twist (see also Refs. Jansen:2005cg ; Dimopoulos:2009qv ). Including the clover term appears to reduce its size drastically. We refer to Ref. Herdoiza:2013sla for a systematic investigation of this splitting for simulations without the clover term.

V.2 $\eta_{2}$ meson mass from fermionic correlators

In Figure 3 we show $C_{\eta_{2}}$ and its shifted version $\tilde{C}_{\eta_{2}}$ as functions of Euclidean time $t/a$, in the left panel for the physical point ensemble cA2.09.48 and in the right panel for cA2.60.32. For the physical point (left panel) we observe a sign change in $C_{\eta_{2}}$ around $t/a=8$. However, from even slightly earlier values of $t/a$, the correlation function is compatible with zero, at least within two sigma. The point errors are large compared to the observed fluctuations between different $t/a$ values indicating large correlations. The shifting has two effects. First, the error bars are dramatically decreased in $\tilde{C}_{\eta_{2}}$ compared to $C_{\eta_{2}}$ with at the same time strongly reduced correlations. Second, $\tilde{C}_{\eta_{2}}$ turns negative only at $t/a=18$ and stays compatible with zero within two sigma from then on.

In the right panel of Figure 3 for $M_{\pi}\approx 340\ \mathrm{MeV}$ the unshifted correlation function does not show a sign change. Still, the shifted correlation function exhibits significantly smaller error bars due to largely reduced correlations.

In Figure 4 we focus on the physical point ensemble cA2.09.48. We show in a half logarithmic plot the connected, disconnected and full $\eta_{2}$ correlation function versus $t/a$. We recall that the full correlation function is obtained as the difference between the connected and disconnected contribution, cf. Eq. (8). In the left panel we show the unshifted correlators and in the right panel the corresponding shifted correlators. While the observations are the same as obtained from Figure 3, the effect of the shift is better visible due to the logarithmic scale on the $y$ axis.

Moreover, one sees from Figure 4 that the signal-to-noise ratio of the connected only contribution stays approximately constant until close to $t=T/2$. Therefore, the connected correlation function can be fitted at large Euclidean times using the ansatz

[TABLE]

where the $\pm$ depends on whether the shifted or unshifted correlation function is analyzed. Additionally, one learns from Figure 4 that the error on the full correlation functions mainly stems from the disconnected contribution.

Once the connected-only part is fitted with the ansatz above, we can apply the excited state subtraction as explained earlier. We denote the corresponding subtracted and shifted $\eta_{2}$ correlation function as $\tilde{C}^{\mathrm{sub}}_{\eta_{2}}(t)$. In Figure 5 we show the effective masses computed from $\tilde{C}^{\mathrm{sub}}_{\eta_{2}}(t)$ as a function of $t/a$. In the left panel we show the data for the physical point ensemble cA2.09.48, in the right panel for cA2.60.32. In both cases we observe a plateau in the effective masses from $t/a=2$ or even $t/a=1$ on. The result of a fit to the correlation function is indicated by the horizontal lines, indicating also the fit range. The end points of the fit ranges lie outside the plotted region, because we obtain a signal in the correlation function further out in $t/a$. For comparison, we also show the effective masses computed from $\tilde{C}_{\eta_{2}}$ without excited state subtraction, for which a plateau can clearly not be identified with confidence.

The final values for $M_{\eta_{2}}$ are determined from a fit of ansatz Eq. (21) to $\tilde{C}^{\mathrm{sub}}_{\eta_{2}}(t)$. The corresponding results are compiled in Table 2.

V.3 $\eta_{2}$ Meson mass from topological charge density

correlators

When determining the fit range in fitting the form in Eq. (14) to the topological charge density correlators $C_{qq}$, one needs to take into account the range over which the contact term is smeared, as discussed above. For this reason, we show in the left panel of Figure 6 the correlators on ensemble cA2.09.48 for different smearing levels $n=15,30,45,60,75,90$. Since the maximum of the correlator at distance $r=0$ is suppressed with an increased smearing level and varies by an order of magnitude between smearing levels $n=15$ and $n=90$, we normalize the correlators by $C_{qq}(r=0)$.

The smearing range can be described by the two characteristic scales $R_{0}$ and $R_{\text{min}}$, defined by the conditions $C_{qq}(r=R_{0})=0$ and $C_{qq}(r=R_{\text{min}})$ where the correlator has its minimum value. The dependence of these smearing ranges on the smearing levels is displayed in the right panel of Figure 6 for the correlators on ensemble cA2.09.48 together with fits of the form $c_{0}+c_{1}\log(n)+c_{2}\log(n)^{2}$.

In Figure 7 we show the long distance

behavior of the correlators on a logarithmic scale for the various smearing levels. It is comforting to see that the correlators start to asymptotically fall on top of each other for increasing smearing level. Smearing levels $n=75$ and 90 for example are statistically indistinguishable for $r/a\gtrsim 11$. Note that since the asymptotic form of the correlator in Eq. (14) is

[TABLE]

rather than purely exponential, the choice of the optimal fit range cannot be guided by an effective mass plot. From Figure 7 we infer that for the lowest smearing level $n=15$ the signal is essentially lost after $r/a\gtrsim 16$, while for the largest smearing level $n=90$ the fit range is limited to $r\gtrsim 12$ due to the contamination by the smeared-out contact term. Consequently, the intermediate smearing levels seem to provide the longest fit ranges when both restrictions are taken into account.

When trying to maximize the fit range $[r_{\text{min}},r_{\text{max}}]$ for the different smearing levels, we notice that the fit results are not particularly sensitive to the choice of $r_{\text{max}}$ as long as $r_{\text{max}}\gtrsim 16$. On the other hand, the error depends strongly on the choice of $r_{\text{min}}$. This is illustrated in the left panel of Fig. 8 where we show the

fit results for $aM$ as a function of $r_{\text{min}}/a$ while keeping $r_{\text{max}}/a=20$ fixed. As we lower $r_{\text{min}}$ the error becomes smaller, but at some point the fit result starts to change due to the influence of the smeared contact term, and possibly also excited state contributions. Consequently, for each smearing level we minimize $r_{\text{min}}/a$ while making sure that the result is still stable under a variation of $r_{\text{min}}/a$. In Fig. 9 we give an

example for such a fit. The left panel shows the correlation function on ensemble cA2.09.48 at smearing level $n=45$ together with the fit function from a fit using $r_{\text{min}}/a=10$ and $r_{\text{max}}/a=20$, while the right panel shows the differences between the fit function and the data points. In this example we get $aM_{\eta_{2}}=0.3791(71)$ and $\chi^{2}/\text{dof}=0.61$ with 299 (correlated) degrees of freedom. This result is very stable under a large variation of the fit range.

Our choice for the $r_{\text{min}}/a$ values are $r_{\text{min}}/a\sim 9.5,10.0,10.0,11.0,12.0,13.0$ for smearing level $n=15,30,45,60,75,90$, respectively, and in the right panel of Fig. 8 we display the final fit result for each smearing level.

Finally, we choose as our final value the weighted average between the three smearing levels $n=60,75$ and 90, at which the fit results seem to stabilize, and we use the statistical error from the result at level $n=75$ which also roughly covers the systematic error from varying $n$. Our final result

[TABLE]

is displayed in the right panel of Fig. 8 as the vertical orange band. We note that this is well compatible with the result from the fermionic correlators in Sec. V.2, but it is here obtained from smeared topological charge density correlators which are significantly cheaper to calculate.

Repeating this procedure for the other ensembles yields the results for $aM_{\eta_{2}}^{\text{gl.}}$ compiled in Table 2. We note that the values on the smaller lattice volumes have a significantly larger error. This is mainly due to two reasons. First and foremost, the calculations on the smaller lattices cannot benefit from self-averaging as much as the ones on the larger lattices. Second, due to the smaller lattice extent, the fitting ranges, in particular $r_{\text{max}}$, are more restricted leading to a larger variation of the fitted masses with the smearing levels and hence to a larger systematic error.

V.4 Topological susceptibility

In Table 2 we have also compiled our results for the topological susceptibility evaluated at flow time $t=3t_{0}$ as discussed in Sec. III.2, and the gradient flow scale $t_{1}/a^{2}$. The values for $t_{0}/a^{2}$ can be found in Ref. Abdel-Rehim:2015pwa . We express the susceptibility in units of $t_{1}$ in order to facilitate comparison with Ref. Bruno:2014ova and display the values in Figure 10 as a function of $t_{1}M_{\pi}^{2}$. In leading order Wilson chiral perturbation theory one expects the following dependence of $\chi_{\text{top}}$ on the lattice spacing and the pion mass Bruno:2014ova written in units of the gradient flow scale $t_{1}$:

[TABLE]

Apart from the ensemble with too small volume cA2.30.24, our data are nicely compatible with this expectation: the solid line in Figure 10 represents a fit of the function

[TABLE]

to the data with fit parameters $c_{1}$ and $c_{2}$. Note that we use the charged pion mass, because charged and neutral pion masses are degenerate within errors. The best fit parameter for $c_{1}$ is compatible with $t_{1}f_{\pi}^{2}/8$. Note that ensemble cA2.30.24 has a very small volume explaining the outlier in Figure 10.

The fitted value for $c_{2}$ can be compared to the results of Ref. Bruno:2014ova using Wilson clover fermions. They obtain $c_{2}=5.1(7)\times 10^{-3}$, while our value reads $c_{2}=2.86(26)\times 10^{-3}$ indicating a sizable reduction of the corresponding lattice artifact.

VI Discussion

In Figure 11 we show $M_{\eta_{2}}^{\mathrm{ferm.}}$ in units of the Sommer parameter $r_{0}$ as a function of $(r_{0}M_{\pi})^{2}$, with the value of $r_{0}/a=5.317(48)$ from ensemble cA2.09.48 taken from Ref. Abdel-Rehim:2015pwa . The outlier in our data points stems again from the ensemble cA2.30.24, which has a very small value of $M_{\pi}L$. We compare the results presented in this paper determined from the fermionic correlators to other lattice determinations available in the literature: the two UKQCD results stem from Refs. Allton:2004qq , the PACS-CS result from Ref. Lesk:2002gd and the DWF result from Ref. Hashimoto:2008xg . The twisted mass results without clover are taken from Ref. Jansen:2008wv .

From this figure we conclude that there is overall very good agreement between the different determinations. Even if the different investigations do not cover a wide range in the lattice spacing, there is no room for sizable lattice artifacts. The results presented in this work complete the picture toward the physical point, with a value

[TABLE]

at the physical point. Using $r_{0}=0.4907(86)$ from Ref. Abdel-Rehim:2015pwa we arrive at

[TABLE]

where the scale setting error has been propagated into the final error estimate. While the result is a bit lower than what is quoted in Ref. Jansen:2008wv , the flat dependence of $M_{\eta_{2}}$ on the light quark mass is confirmed. Interestingly, this value agrees very well with an estimate from Ref. McNeile:2000hf , where a phenomenological analysis of the full $\eta,\eta^{\prime}$ mixing matrix has been performed to arrive at $M_{\eta_{2}}\approx 776\ \mathrm{MeV}$.

With this determination of $M_{\eta_{2}}$ at the physical pion mass value it is almost certain that the $\eta_{2}$ meson will have a finite mass in the chiral limit, agreeing with the picture that the $\eta_{2}$ is not a Goldstone boson. It implies that the topological susceptibility must decrease as $M_{\pi}^{2}$ toward the chiral limit Leutwyler:1992yt .

VII Summary

In this paper we have presented results for the $\eta_{2}$ meson related to the axial anomaly and the topological susceptibility in two-flavor QCD. The results have been obtained using $N_{f}=2$ lattice QCD ensembles generated by ETMC with the Wilson twisted clover discretization Abdel-Rehim:2015pwa . Pion mass values reach from the physical value up to $340\ \mathrm{MeV}$ at a single lattice spacing value of $a=0.0931(2)\ \mathrm{fm}$. For the $\eta_{2}$ we could confirm the almost constant extrapolation in $M_{\pi}^{2}$ toward the physical point. Errors are significantly reduced compared to previous calculations. Lattice artifacts seem to be not larger than our statistical uncertainty.

Regarding a future study of the $\eta$ and $\eta^{\prime}$ at physical quark masses in the $N_{f}=2+1+1$ theory we conclude that such a calculation should now be feasible assuming a roughly similar signal-to-noise ratio as in the two-flavor case. Since it is known from earlier $N_{f}=2+1+1$ simulations at unphysical quark masses that the total error is dominated by the error on the light quark disconnected loops, such an assumption seems reasonable. While the nondegenerate heavy quark doublet will require additional inversions, it should only lead to a moderate increase in the total computational cost. An additional complication in the $N_{f}=2+1+1$ case arises from the technically more involved analysis because — unlike the $\eta_{2}$ — the $\eta^{\prime}$ is not a ground state. However, all the relevant analysis methods have been developed and successfully applied previously in Refs. Michael:2013gka ; Ottnad:2017bjt in a study of the $\eta$, $\eta^{\prime}$ at unphysical quark masses, and the analysis at physical quark masses can be done in the same way.

We complement the determination of $M_{\eta_{2}}$ at the physical point from fermionic correlation functions with one from the topological charge density correlator. We find that with the number of APE smearing steps larger than or equal to $60$ the estimated value of $M_{\eta_{2}}$ becomes stable. The so determined value for $M_{\eta_{2}}$ is fully compatible with the one from fermionic correlators and has an even smaller statistical uncertainty. It is straightforward to apply this methodology in the $N_{f}=2+1+1$ theory in order to determine the mass of the $\eta^{\prime}$ meson: except for the mixing with the $\eta$, which can be taken into account by appropriately modifying the fit function, we do not expect any additional complications.

The topological susceptibility has been computed using the gradient flow. As expected, $\chi_{\text{top}}$ is proportional to $M_{\pi}^{2}$ (for small $M_{\pi}^{2}$) up to an additive lattice artifact independent of $M_{\pi}$. Even if we are not able to finally confirm this with only a single lattice spacing at hand, this constant term should be of $\mathcal{O}(a^{2})$. The size of this artifact appears to be significantly smaller than what is observed with Wilson clover fermions in Ref. Bruno:2014ova .

Acknowledgements.

We thank the members of ETMC for the most enjoyable collaboration. We thank G. Rossi for valuable comments on the manuscript. The computer time for this project was made available to us by the John von Neumann-Institute for Computing (NIC) on the Jureca and Juqueen systems in Jülich and the HPC Cluster in Bern. This project was funded by the DFG as a project in the Sino-German CRC110. U. W. acknowledges support from the Swiss National Science Foundation. The open source software packages tmLQCD Jansen:2009xp ; Deuzeman:2013xaa ; Abdel-Rehim:2013wba , Lemon Deuzeman:2011wz , DD$\alpha$AMG Alexandrou:2016izb , and R R:2005 have been used.