Topological susceptibility and η ′ meson mass from N f =2 lattice QCD at the physical point
Petros Dimopoulos, Christopher Helmes, Christian Jost, Bastian Knippschild, Bartosz Kostrzewa, Liuming Liu, Konstantin Ottnad, Marcus Petschlies, Carsten Urbach, Urs Wenger, Markus Werner
TTopological susceptibility and η (cid:48) meson massfrom N f = 2 lattice QCD at the physical point P. Dimopoulos, C. Helmes, C. Jost, B. Knippschild, B. Kostrzewa, L. Liu, K. Ottnad,
2, 4
M. Petschlies, C. Urbach, U. Wenger, and M. Werner (ETM Collaboration) Dipartimento di Fisica, Universit`a di Roma Tor Vergata, 00133 Roma, Italy HISKP and BCTP, Rheinische Friedrich-Wilhelms Universit¨at Bonn, 53115 Bonn, Germany Institute of Modern Physics, Chinese Academy of Sciences, 730000 Lanzhou, China Institut f¨ur Kernphysik, Johann-Joachim-Becher-Weg 45, University of Mainz, 55099 Mainz, Germany Institute for Theoretical Physics, Albert Einstein Center for Fundamental Physics,University of Bern, 3012 Bern, Switzerland
In this paper we explore the computation of topological susceptibility and η (cid:48) meson massin N f = 2 flavor QCD using lattice techniques with physical value of the pion mass as wellas larger pion mass values. We observe that the physical point can be reached without asignificant increase in the statistical noise. The mass of the η (cid:48) meson can be obtained fromboth fermionic two point functions and topological charge density correlation functions,giving compatible results. With the pion mass dependence of the η (cid:48) mass being flat wearrive at M η (cid:48) = 772(18) MeV without an explicit continuum limit. For the topologicalsusceptibility we observe a linear dependence on M π , however, with an additional constantstemming from lattice artifacts. a r X i v : . [ h e p - l a t ] M a r I. INTRODUCTION
Due to the persisting 3 − σ deviation in the anomalous magnetic moment of the muon a µ between theory and experiment there is considerable interest in the decays η → γ (cid:63) γ (cid:63) and η (cid:48) → γ (cid:63) γ (cid:63) because a better knowledge of the corresponding transition form factors could help to reduce theuncertainty in the hadronic light-by-light contribution to a µ ; see for instance Ref. [1]. Moreover, η and η (cid:48) mesons are interesting from a theoretical point of view because the large mass of the η (cid:48) meson is explained by the anomalously broken U A (1) axial symmetry in QCD. The η, η (cid:48) mixingpattern and the aforementioned transition form factors can be computed nonperturbatively usinglattice techniques.There has been considerable progress in studying η and η (cid:48) mesons from lattice QCD. In Refs. [2,3] 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 tothe physical pion mass value excellent agreement to experiment was found. Further lattice resultsfor η, η (cid:48) can be found in Refs. [4–8].Through the anomaly, the mass of the η (cid:48) is also tightly connected to topology and in particularthe topological susceptibility χ top . The latter quantity must decrease as M π toward the chiral limit,if the η (cid:48) is not a Goldstone boson [9]. For recent lattice studies of the topological susceptibility;see for instance [10, 11]. There is now particular interest in χ top due to its connection to axiondark matter; see for instance Refs. [12–14].In this paper we attempt to study the η (cid:48) meson and the topological susceptibility directly at thephysical point, however, in a first step in N f = 2 flavor QCD. In N f = 2 flavor QCD there exist apion triplet and one flavor singlet, which is related to the aforementioned anomaly. We will denoteit as the η meson to distinguish it from the η (cid:48) meson in full QCD, which is only approximately aflavor eigenstate. The η and the η (cid:48) meson have in common that both receive significant fermionicdisconnected contributions. In an earlier study [15] their masses have been found to differ onlyby 200 MeV, with the additional strange quark introducing only a moderate shift in the mass. Inparticular, both are expected to have a similar dependence on the light quark mass. The mostrecent lattice QCD studies of the η meson can be found in Refs. [16, 17].We investigate M η using fermionic correlation functions and in addition topological chargedensity correlators. The topological susceptibility is studied using gradient flow techniques [18].Studying the η meson and the topological susceptibility at the physical point will reveal on theone hand important qualitative information on the implementation of the anomaly in QCD. On Ensemble β c sw aµ (cid:96) ( L/a ) × T /a N conf cA . .
48 2.10 1.57551 0.009 48 ×
96 615 cA . .
48 2.10 1.57551 0.030 48 ×
96 352 cA . .
24 2.10 1.57551 0.030 24 ×
48 352 cA . .
32 2.10 1.57551 0.060 32 ×
64 337 cA . .
24 2.10 1.57551 0.060 24 ×
48 424Table I: The gauge ensembles used in this study. The labeling of the ensembles follows the notations inRef. [21]. In addition to the relevant input parameters we give the lattice volume (
L/a ) × T /a and thenumber of evaluated configurations N conf . the other hand it represents a feasibility study for a later investigation of η and η (cid:48) in N f = 2 + 1 + 1QCD at the physical point [19]. The results obtained here are also important prerequisites for anexploratory study of η → γ (cid:63) γ (cid:63) .The paper is organized as follows: in the following two sections we discuss the lattice details ofour computation. In Sec. IV we present the analysis methods and in Sec. V the results. We closewith a discussion and summary. For a first account of this work we refer to Ref. [20]. II. LATTICE ACTION
The results presented in this paper are based on the gauge configurations generated by theETMC with Wilson clover twisted mass quark action at maximal twist [21]. We employ theIwasaki gauge action [22]. The measurements are performed on a set of N f = 2 ensembles withthe pion mass ranging from its physical value to 340 MeV. In Table I we list all the ensemblestogether with the relevant input parameters, the lattice volume, and the number of configurations.The lattice spacing is a = 0 . D (cid:96) = D − iγ τ (cid:20) W cr + i c sw σ µν F µν (cid:21) + µ (cid:96) , (1)which acts on a flavor doublet spinor ψ = ( u, d ) T . In Eq. (1) we have D = γ µ ( ∇ ∗ µ + ∇ µ ) / ∇ µ and ∇ ∗ µ the forward and backward lattice covariant derivatives, and the Wilson term W cr = − ra ∇ ∗ µ ∇ µ + m cr with the critical mass m cr , the Wilson parameter r = 1, and the latticespacing a . The average up/down (twisted) quark mass is denoted by µ (cid:96) , while c sw is the so-calledSheikoleslami-Wohlert improvement coefficient [24] multiplying the clover term. It is in our casenot used for O ( a ) improvement but serves to significantly reduce the effects of isospin breaking [21].The critical mass has been determined as described in Refs. [25, 26]. This guarantees automatic O ( a ) improvement [27], which is one of the main advantages of the Wilson twisted mass formulationof lattice QCD. III. OBSERVABLES
As a smearing scheme in the computation of fermionic correlation functions we use the stochasticLaplacian Heaviside (sLapH) method [28, 29]. The details of our sLapH parameter choices for aset of N f = 2 + 1 + 1 Wilson twisted mass ensembles are given in Ref. [30]. The parameters forthe ensembles used in this work are the same as those for N f = 2 + 1 + 1 ensembles with thecorresponding lattice volume. A. η and pion correlation functions In N f = 2 flavor QCD there is the neutral pion, corresponding to the neutral of the three pionsin the triplet, and the η , 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 andthe η with our action. We employ the following pseudoscalar interpolating operators projected tozero momentum, which are all local and Hermitian P ( t ) = 1 √ (cid:88) x ¯ ψiγ τ ψ ( x , t ) , P ( t ) = 1 √ (cid:88) x ¯ ψiγ f ψ ( x , t ) . (2)Here, τ is the third Pauli and f the unit matrix, both acting in flavor space. From those onebuilds the correlation functions C π ( t − t (cid:48) ) = (cid:104)P ( t ) ( P ( t (cid:48) )) † (cid:105) , (3) C η ( t − t (cid:48) ) = (cid:104)P ( t ) ( P ( t (cid:48) )) † (cid:105) , (4)which allow one to determine the masses M π and M η from their decay in Euclidean time. Bothcorrelation functions in Eqs. (3) and (4) do have a fermionic connected and a fermionic disconnectedcontribution, the latter of which vanishes exactly in case of the neutral pion in an isospin symmetrictheory. Since this is not the case for Wilson twisted mass fermions, we have to take the disconnectedcontributions into account also for the π .For the disconnected part of C η we consider the loop (cid:104) ¯ ψ u iγ ψ u ( x ) + ¯ ψ d iγ ψ d ( x ) (cid:105) F = − i Tr { γ G xxu } − i Tr { γ G xxd } = − i Tr { γ G xxu } − i Tr { ( G xxu ) † γ } = − i ReTr { γ G xxu } . (5)Here, we have used the γ hermiticity property D d = γ D † u γ . G xyu/d represents the up or down propagator. Similarly, one shows for C π (cid:104) ¯ ψ u iγ ψ u ( x ) − ¯ ψ d iγ ψ d ( x ) (cid:105) F = 2 ImTr { γ G xxu } . (6)The fermionic connected contribution is identical for the two correlation functions Eq. (3) andEq. (4). With similar arguments as for the loops one finds C conn ( t − t (cid:48) ) = ReTr { γ G tt (cid:48) u γ G t (cid:48) tu } , (7)where we have suppressed the spatial indices. From Eqs. (5), (6), and (7) we infer the expressionsfor the π and η correlation functions as follows: C π ( t − t (cid:48) ) = Tr { γ G tt (cid:48) u γ G t (cid:48) tu } + 2 ImTr { γ G ttu } · ImTr { γ G t (cid:48) t (cid:48) u } ,C η ( t − t (cid:48) ) = Tr { γ G tt (cid:48) u γ G t (cid:48) tu } − { γ G ttu } · ReTr { γ G t (cid:48) t (cid:48) u } . (8)For completeness, the correlation function of the charged pion is constructed as C π ± ( t − t (cid:48) ) = (cid:104)P + ( t ) ( P + ( t (cid:48) )) † (cid:105) (9)with P + ( t ) = (cid:88) x ¯ ψiγ τ + iτ ψ ( x , t ) (10)and τ and τ the first and second Pauli matrices, respectively. B. Topological charge density correlations and susceptibility
The naive field theoretical definition of the topological charge density given by q ( x ) = − π ε µνρσ Tr F µν ( x ) F ρσ ( x ) (11)defines the topological charge density point-to-point correlator C qq ( x − y ) = (cid:104) q ( x ) q ( y ) (cid:105) . (12)The topological susceptibility, which is a measure for the fluctuations of the topological charge, isdefined as χ top = 1 V (cid:90) d x (cid:90) d y (cid:104) q ( x ) q ( y ) (cid:105) (13)where V is the spacetime volume. Due to the pseudoscalar nature of the topological charge densitythe topological charge density correlator is strictly negative for finite separations, C qq ( x − y > <
0. On the other hand, it is clear that the susceptibility is strictly positive, because V · χ top = (cid:104) Q (cid:105) >
0, where Q = (cid:82) dx q ( x ) is the total topological charge of the gauge field. This apparentcontradiction is resolved by recalling that C qq suffers from contact-term singularities at x − y = 0which need to be renormalized in order for the susceptibility to make physical sense. Hence, thephysics 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 | > η . As a consequence, the behavior of the topological charge density point-to-point correlator is dominated by the single boson propagator for the η meson and thereforefollows the form of the scalar propagator [12, 31] C qq ( x, y ) ∼ M | x − y | K ( M · | x − y | ) (14)where K is the modified Bessel function of the first kind, and M is the mass of the lightest particlein the pseudoscalar meson sector, i.e., the mass M η of the η meson.On the lattice the topological charge density is discretized with a clover-type discretizationof the field strength tensor F µν which extends over a distance of 2 a in lattice units. Hence thecontact-term contributions to the topological charge density correlator are also distributed overthe distance | x − y | ∼ a . Moreover, since the discretized field strength tensor is constructed fromsmeared gauge links in order to remove ultraviolet fluctuations of the gauge field, the positivecontact-term contributions to the correlator are spread over a range R which depends on thedetails of the smoothing scheme. However, the behavior of the correlator for | x − y | (cid:29) R shouldbe independent of the details of the smearing scheme and hence any smoothing scheme is supposedto yield the same physics, i.e. the same mass M η .For the topological charge density correlator we use the array processor experiment (APE)smearing scheme [32] with various smearing levels ranging up to 90 iterative smearing steps. Thisis in order to check for the independence of the results from the smearing scheme. To compute thetopological susceptibility χ top we employ the gradient flow technique as introduced for lattice QCDin Ref. [18]. It has the advantage of yielding a renormalized topological susceptibility at finite flowtime t [33], in particular it renormalizes the contact term singularities in the continuum limit at anyfixed, physical value of t . Since the renormalized susceptibility is scale invariant, i.e., independentof the renormalization scale, χ 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, wenote 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 firstobtain the topological charge at flow time t from the topological charge density q t ( x ) evaluated onthe flown gauge field configuration, Q ( t ) = a (cid:88) x q t ( x ) , (15)and then the susceptibility via χ top ( t ) = (cid:104) Q ( t ) (cid:105) /V . We choose t = 3 t where t is the usual gradi-ent flow reference scale defined through the renormalized action density [18] on the correspondingensemble. In addition, we also make use of the related reference scale t in order to facilitate thecomparison of our results with those in Ref. [34]. IV. ANALYSIS METHOD
In the following Secs. IV A and IV B we will first give more details on the analysis of the fermioniccorrelation function, before we turn to the discussion for the gluonic correlators in Sec. IV C.The fermionic correlation function data are generally analyzed using the blocked bootstrapprocedure with R = 10000 bootstrap samples. Depending on the ensemble, we have chosen theblock size such that at least (cid:38)
100 blocked data points are left. The relevant masses are computedfrom correlated fits to the correlation function data.The gluonic correlation function data are analyzed using a jackknife procedure. It turns outthat the correlators at separate distances are highly correlated even for large separations, such thatthe covariance matrix cannot be taken into account reliably in the fitting procedur;e see furtherdetails below. The data for the topological charge susceptibility (and the gluonic scales t and t ) are analyzed using a blocked bootstrap procedure with R = 1000 bootstrap samples and blocksizes such that (cid:38)
30 blocked data points are left. The so obtained error is compared to the naiveone corrected by the integrated autocorrelation time τ int , and the larger of the two is always chosenas the final error. Since these calculations are inexpensive, we use at least double the number ofconfigurations indicated in Table I. A. Excited state subtraction
In particular for the η meson, the fermionic disconnected contributions are very noisy. As aconsequence, the signal is lost relatively early in Euclidean time. For this reason we have in thepast applied a method to subtract excited states [2, 3, 17, 35], originally proposed in Ref. [36]. Itactually works very well and we will apply it here again for the η meson. It consists of subtractingexcited 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) tothe full η correlation functions.The underlying assumption is that disconnected contributions are large for the ground state,i.e. the η , but not for excited states. If this assumption is correct, the effective mass M eff = − log C η ( t ) C η ( t + 1) (16)should show a plateau from very early Euclidean times on. We have found in Refs. [2, 3, 17] thatthis approach works very well for the η meson in N f = 2 flavor QCD as well as for η and η (cid:48) mesonsin N f = 2 + 1 + 1 flavor QCD. B. Shifted correlation functions
The expected time dependence of the fermionic correlation functions considered here reads asfollows: C ( t ) = |(cid:104) |O| (cid:105)| + (cid:88) n |(cid:104) |O| n (cid:105)| E n (cid:16) e − E n t + e − E n ( T − t ) (cid:17) , (17)where O = P + , P , P and n labels the states with the corresponding quantum numbers. Thetime 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 P + and P the VEV must bezero, while this is not the case for P . We deal with the VEV by building the shifted correlationfunction ˜ C ( t ) = C ( t ) − C ( t + 1) . (18)The difference cancels the constant VEV contribution, while also changing the time dependence tobe antisymmetric in time, ˜ C ( t ) ∝ (cid:16) e − E n t − e − E n ( T − t ) (cid:17) . (19)As an alternative, one can also compute the VEV |(cid:104) |O| (cid:105)| from the data and subtract it explicitly.Since the VEV has to be zero for P up to statistical fluctuations, strictly speaking we do notneed to use the shifting procedure for the η meson. However, as has been argued in Ref. [37] andfirst investigated in Ref. [38], there is an additive finite volume effect to C η constant in Euclideantime of the form ∝ a T (cid:18) χ top + Q V (cid:19) (20)proportional to the topological susceptibility χ top and the squared topological charge Q . If present,such a term will cause the η correlation function to stay finite at large Euclidean times. Dependingon the sign of the coefficient in front of the finite volume effect, the correlation function may eventurn negative at relatively small Euclidean times. Clearly, a finite volume effect of this type canbe subtracted again using the shifting procedure, which has first been proposed and applied inRef. [3]. C. Topological charge density correlators
For the computation of the topological charge density correlator we make use of the full transla-tional 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, calcu-lating the correlator in Fourier space and transforming it back to coordinate space. In this sensethe evaluation is exact, in contrast to the computation of the disconnected contributions to thefermionic 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 correlationfunction, i.e. in this case the contribution of the η state. Third, with an increasing smearing levelthe contact term is distributed over larger distances and hence distorts the correlation function upto larger and larger separations. Obviously, these effects compete with each other with respect tothe optimal fit range.0In principle, the choice of the fit ranges should be determined by the quality of the fits. Unfor-tunately, here this is not possible, because the correlators at separate distances r and r (cid:48) are highlycorrelated. We illustrate this in Figure 1 where we show the covariance Cov( C qq ( r ) , C qq ( r (cid:48) )) of thecorrelation functions as a function of ( r − r (cid:48) ) /a for different values of r and smearing level n = 90. We are essentially looking at separate columns of the covariance matrix. For ease of comparison C ov ( C qq (r) C qq (r ' )) r/a= 6.5r/a= 7.5r/a= 8.5r/a= 9.5r/a=10.5r/a=11.5r/a=12.5 Figure 1: Covariance of the binned correlation functions as a function of ( r − r (cid:48) ) /a for different values of r and smearing level n = 90 on ensemble cA2.09.48. we normalize the covariance by Cov( C qq ( r ) , C qq ( r )), and in addition bin the data into bins of size∆ r/a = 0 . C qq ’s are positively correlated for ( r (cid:48) − r ) /a (cid:46) . r (cid:48) − r ) /a ∼ . r , the columns of the covariance matrix arehighly linear dependent and the matrix itself is very ill-conditioned. As a consequence, it cannotbe taken into account for reliably estimating the quality of the χ -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 ofthe gauge fields. The data for the other smearing levels look very similar. Ensemble aM π aM π aM π c aM ferm. η aM gl. η t /a · t χ top cA . .
48 0 . . . . . . . cA . .
48 0 . . . . . . . cA . .
24 0 . . . . . . . cA . .
32 0 . . . . . . . cA . .
24 0 . . . . . . . η meson in lattice units (from fermionic and gluonic correlators), the gluonic gradient flow lattice scale t /a ,and the topological susceptibility in units of t for the five ensembles considered. π π π π ± t/a ˜ C ( t ) L = 24 aL = 32 aL = 48 aM π / MeV M π / M e V Figure 2: Overview of results for the pions. Left panel: Shifted correlation functions ˜ C ( t ) for the charged andthe neutral pion on ensemble cA2.09.48 with physical quark mass. In the case of the neutral pion we showthe full correlation function as well as individual quark-connected and quark-disconnected contributions.Right panel: Mass of the (full) neutral pion as a function of the charged pion mass. V. RESULTS
In Table II we show results for the pion and η ferm.2 masses which have been computed fromfermionic correlation functions, the η gl.2 masses obtained from the gluonic topological charge densitycorrelation functions, as well as the gluonic gradient flow lattice scale t /a and the topologicalsusceptibility χ top . For the charged pion we will always use the shorthand M π , while for the neutralpion M π and M π c refer to the full and quark-connected masses, respectively. In the following wediscuss these results in more detail.2 A. Neutral pion
In contrast to the η meson discussed later, the signal for the neutral pion can be resolved forall values of t/a . In the left panel of Figure 2 we show the shifted correlation function ˜ C ( t ) for theneutral pion as well as the individual quark-connected and quark-disconnected contributions. Notethat in this case the function shift is required to remove the offset from the vacuum expectationvalue in the quark-disconnected contribution. Clearly the signal in the quark-disconnected part iswell behaved even for the largest values of t/a . For comparison the charged pion has been includedin the plot as well.As already visible from the left panel of Figure 2, charged and neutral pions appear to havevery similar mass values. This is even more apparent from the right panel where the neutral pionmass values are plotted versus the charged pion mass values, both in physical units, for all theensembles considered here. The points fall almost on the bisecting line, which indicates no masssplitting between neutral and charged pion mass. This finding, which we pointed out already inRef. [21], is rather important: this mass splitting is basically the only large a lattice artefact thatwas found for simulations with Wilson twisted fermions at maximal twist (see also Refs. [39, 40]).Including the clover term appears to reduce its size drastically. We refer to Ref. [41] for a systematicinvestigation of this splitting for simulations without the clover term. B. η meson mass from fermionic correlators In Figure 3 we show C η and its shifted version ˜ C η as functions of Euclidean time t/a , in theleft panel for the physical point ensemble cA2.09.48 and in the right panel for cA2.60.32. For thephysical point (left panel) we observe a sign change in C η around t/a = 8. However, from evenslightly earlier values of t/a , the correlation function is compatible with zero, at least within twosigma. The point errors are large compared to the observed fluctuations between different t/a valuesindicating large correlations. The shifting has two effects. First, the error bars are dramaticallydecreased in ˜ C η compared to C η with at the same time strongly reduced correlations. Second,˜ C η turns negative only at t/a = 18 and stays compatible with zero within two sigma from thenon.In the right panel of Figure 3 for M π ≈
340 MeV the unshifted correlation function does notshow a sign change. Still, the shifted correlation function exhibits significantly smaller error barsdue to largely reduced correlations.3 ˜ C η ( t ) C η ( t ) t/a C η ( t ) ˜ C η ( t ) C η ( t ) t/a C η ( t ) Figure 3: η correlation function C η and its shifted version ˜ C η as a function of t/a . For better visibility ofthe tail some of the numerically very large data points at small values of t/a are not included in the plot.Left panel: Ensemble cA2.09.48. Right panel: Ensemble cA2.60.32. disconnectedconnectedfull t/a C η ( t ) disconnectedconnectedfull t/a ˜ C η ( t ) Figure 4: Connected, disconnected and full η correlation function versus t/a for the physical point ensemblecA2.09.48. Left panel: Original correlation function C η . Right panel: Shifted correlation function ˜ C η . In Figure 4 we focus on the physical point ensemble cA2.09.48. We show in a half logarithmicplot the connected, disconnected and full η correlation function versus t/a . We recall that thefull correlation function is obtained as the difference between the connected and disconnectedcontribution, cf. Eq. (8). In the left panel we show the unshifted correlators and in the right4 from ˜ C sub η ( t )from ˜ C η ( t )fitted mass t/a a M η C sub η ( t )from ˜ C η ( t )fitted mass t/a a M η Figure 5: Lattice data for effective masses computed from ˜ C η without and from ˜ C sub η with excited statessubtracted. The result for the mass and its error from a correlated fit to the correlator data of ˜ C sub η has beenincluded. Note that the end points of the fit ranges lie outside of the plots, because we fit to the correlationfunction and not to the effective masses. Left panel: Ensemble cA2.09.48. Right panel: Ensemble cA2.60.32. panel the corresponding shifted correlators. While the observations are the same as obtained fromFigure 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 contributionstays approximately constant until close to t = T /
2. Therefore, the connected correlation functioncan be fitted at large Euclidean times using the ansatz f ± ( t ; A, M ) = A (cid:16) e − Mt ± e − M ( T − t ) (cid:17) , (21)where the ± depends on whether the shifted or unshifted correlation function is analyzed. Addi-tionally, one learns from Figure 4 that the error on the full correlation functions mainly stems fromthe disconnected contribution.Once the connected-only part is fitted with the ansatz above, we can apply the excited statesubtraction as explained earlier. We denote the corresponding subtracted and shifted η correlationfunction as ˜ C sub η ( t ). In Figure 5 we show the effective masses computed from ˜ C sub η ( t ) as a functionof t/a . In the left panel we show the data for the physical point ensemble cA2.09.48, in the rightpanel for cA2.60.32. In both cases we observe a plateau in the effective masses from t/a = 2 oreven 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, becausewe obtain a signal in the correlation function further out in t/a . For comparison, we also showthe effective masses computed from ˜ C η without excited state subtraction, for which a plateau canclearly not be identified with confidence.The final values for M η are determined from a fit of ansatz Eq. (21) to ˜ C sub η ( t ). The corre-sponding results are compiled in Table II.5 C. η Meson mass from topological charge density correlators
When determining the fit range in fitting the form in Eq. (14) to the topological charge densitycorrelators 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 onensemble cA2.09.48 for different smearing levels n = 15 , , , , ,
90. Since the maximum ofthe correlator at distance r = 0 is suppressed with an increased smearing level and varies by anorder of magnitude between smearing levels n = 15 and n = 90, we normalize the correlatorsby C qq ( r = 0). The smearing range can be described by the two characteristic scales R and C qq (r) / C qq ( ) APE15APE30APE45APE60APE75APE90 s m ea r i ng r a ng e / a R min R Figure 6: Topological charge density correlator for various APE n smearing with levels n = 15 , , . . . , C qq ( r ) /C qq (0) as a function of theseparation r . Right panel: Scales characterizing the smearing range of the contact term as a function of thesmearing levels. R min , defined by the conditions C qq ( r = R ) = 0 and C qq ( r = R min ) where the correlator has itsminimum value. The dependence of these smearing ranges on the smearing levels is displayed inthe right panel of Figure 6 for the correlators on ensemble cA2.09.48 together with fits of the form c + c log( n ) + c log( n ) .In Figure 7 we show the long distance behavior of the correlators on a logarithmic scale for thevarious smearing levels. It is comforting to see that the correlators start to asymptotically fall ontop of each other for increasing smearing level. Smearing levels n = 75 and 90 for example arestatistically indistinguishable for r/a (cid:38)
11. Note that since the asymptotic form of the correlator6
10 12 14 16 18r/a-12-11-10-9-8-7-6 l n | C qq (r) | APE15APE30APE45APE60APE75APE90
Figure 7: Long distance behavior of the topological charge density correlator for various APE n smearing onensemble cA2.09.48. in Eq. (14) is C qq ( r ) ∼ (cid:114) Mr r e − Mr (cid:18) O (cid:18) M r (cid:19)(cid:19) for large r , (22)rather than purely exponential, the choice of the optimal fit range cannot be guided by an effectivemass plot. From Figure 7 we infer that for the lowest smearing level n = 15 the signal is essentiallylost after r/a (cid:38)
16, while for the largest smearing level n = 90 the fit range is limited to r (cid:38)
12 dueto the contamination by the smeared-out contact term. Consequently, the intermediate smearinglevels seem to provide the longest fit ranges when both restrictions are taken into account.When trying to maximize the fit range [ r min , r max ] for the different smearing levels, we noticethat the fit results are not particularly sensitive to the choice of r max as long as r max (cid:38)
16. On theother hand, the error depends strongly on the choice of r min . This is illustrated in the left panelof Fig. 8 where we show the fit results for aM as a function of r min /a while keeping r max /a = 20fixed. As we lower r min the error becomes smaller, but at some point the fit result starts to changedue to the influence of the smeared contact term, and possibly also excited state contributions.Consequently, for each smearing level we minimize r min /a while making sure that the result is stillstable under a variation of r min /a . In Fig. 9 we give an example for such a fit. The left panelshows the correlation function on ensemble cA2.09.48 at smearing level n = 45 together with the fitfunction from a fit using r min /a = 10 and r max /a = 20, while the right panel shows the differencesbetween the fit function and the data points. In this example we get aM η = 0 . min /a0.350.4 APE90 h Figure 8: Fit results for different APE n smearings on ensemble cA2.09.48. Left panel: As a function of r min /a for fixed r max /a = 20. Right panel: Variation of the fit result with the smearing level.
10 12 14 16 18 20r/a-12-11-10-9-8-7-6 l n | C qq (r) | APE45 10 12 14 16 18 20r/a-0.4-0.200.20.40.6 d l n | C qq (r) | APE45
Figure 9: Example for a fit of the topological charge density correlation function with APE45 smearing onensemble cA2.09.48. Left panel: Fit function and data. Right panel: Differences between fit function anddata. χ / dof = 0 .
61 with 299 (correlated) degrees of freedom. This result is very stable under a largevariation of the fit range.Our choice for the r min /a values are r min /a ∼ . , . , . , . , . , . n = 15 , , , , ,
90, respectively, and in the right panel of Fig. 8 we display the final fit resultfor each smearing level.8Finally, 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 fromthe result at level n = 75 which also roughly covers the systematic error from varying n . Our finalresult aM gl. η = 0 . aM gl. η compiled in TableII. We note that the values on the smaller lattice volumes have a significantly larger error. Thisis mainly due to two reasons. First and foremost, the calculations on the smaller lattices cannotbenefit from self-averaging as much as the ones on the larger lattices. Second, due to the smallerlattice extent, the fitting ranges, in particular r max , are more restricted leading to a larger variationof the fitted masses with the smearing levels and hence to a larger systematic error. D. Topological susceptibility
In Table II we have also compiled our results for the topological susceptibility evaluated at flowtime t = 3 t as discussed in Sec. III B, and the gradient flow scale t /a . The values for t /a canbe found in Ref. [21]. We express the susceptibility in units of t in order to facilitate comparisonwith Ref. [34] and display the values in Figure 10 as a function of t M π . In leading order Wilsonchiral perturbation theory one expects the following dependence of χ top on the lattice spacing andthe pion mass [34] written in units of the gradient flow scale t : t χ top = 18 t f π M π + a c t . (24)Apart from the ensemble with too small volume cA2.30.24, our data are nicely compatible withthis expectation: the solid line in Figure 10 represents a fit of the function g ( M π ) = c t M π + a c t to the data with fit parameters c and c . Note that we use the charged pion mass, becausecharged and neutral pion masses are degenerate within errors. The best fit parameter for c iscompatible with t f π /
8. Note that ensemble cA2.30.24 has a very small volume explaining theoutlier in Figure 10.9 M p t c t op L/a = 48L/a = 24L/a = 32
Figure 10: Topological susceptibility χ top as a function of the squared pion mass, both in appropriate unitsof t . The solid line with shaded error band indicates a fit to the data according to Eq. (24). The fitted value for c can be compared to the results of Ref. [34] using Wilson clover fermions.They obtain c = 5 . × − , while our value reads c = 2 . × − indicating a sizablereduction of the corresponding lattice artifact. VI. DISCUSSION
In Figure 11 we show M ferm .η in units of the Sommer parameter r as a function of ( r M π ) ,with the value of r /a = 5 . M π L .We compare the results presented in this paper determined from the fermionic correlators to otherlattice determinations available in the literature: the two UKQCD results stem from Refs. [42],the PACS-CS result from Ref. [43] and the DWF result from Ref. [44]. The twisted mass resultswithout clover are taken from Ref. [17].From this figure we conclude that there is overall very good agreement between the differentdeterminations. 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 picture0 CLQCDUKQCDDWFCP-PACSETMCthis work,
L/a = 24this work,
L/a = 32this work,
L/a = 48 ( r M π ) r M η Figure 11: Compilation of literature values for the N f = 2 η (cid:48) meson: r M η as a function of ( r M π ) . Thetwo UKQCD results stem from Ref. [42] with the filled symbol for r /a = 5 .
04 and the open symbol for r /a = 5 .
32, the PACS-CS result from Ref. [43] with r /a = 4 .
49, the DWF result from Ref. [44] with r /a = 4 .
28, the CLQCD result from Ref. [16] with r /a = 4 .
22 and the ETMC results from Ref. [17] withfilled symbols for r /a = 5 .
22 and open symbols for r /a = 6 . toward the physical point, with a value r M η = 1 . r = 0 . M η = 772(18) MeVwhere the scale setting error has been propagated into the final error estimate. While the resultis a bit lower than what is quoted in Ref. [17], the flat dependence of M η on the light quarkmass is confirmed. Interestingly, this value agrees very well with an estimate from Ref. [15],where a phenomenological analysis of the full η, η (cid:48) mixing matrix has been performed to arrive at M η ≈
776 MeV.With this determination of M η at the physical pion mass value it is almost certain that the η meson will have a finite mass in the chiral limit, agreeing with the picture that the η is nota Goldstone boson. It implies that the topological susceptibility must decrease as M π toward thechiral limit [9].1 VII. SUMMARY
In this paper we have presented results for the η meson related to the axial anomaly andthe topological susceptibility in two-flavor QCD. The results have been obtained using N f = 2lattice QCD ensembles generated by ETMC with the Wilson twisted clover discretization [21].Pion mass values reach from the physical value up to 340 MeV at a single lattice spacing value of a = 0 . η we could confirm the almost constant extrapolation in M π towardthe physical point. Errors are significantly reduced compared to previous calculations. Latticeartifacts seem to be not larger than our statistical uncertainty.Regarding a future study of the η and η (cid:48) at physical quark masses in the N f = 2+1+1 theory weconclude that such a calculation should now be feasible assuming a roughly similar signal-to-noiseratio as in the two-flavor case. Since it is known from earlier N f = 2+1+1 simulations at unphysicalquark 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 requireadditional 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 involvedanalysis because — unlike the η — the η (cid:48) is not a ground state. However, all the relevant analysismethods have been developed and successfully applied previously in Refs. [2, 3] in a study of the η , η (cid:48) at unphysical quark masses, and the analysis at physical quark masses can be done in thesame way.We complement the determination of M η at the physical point from fermionic correlationfunctions with one from the topological charge density correlator. We find that with the numberof APE smearing steps larger than or equal to 60 the estimated value of M η becomes stable.The so determined value for M η is fully compatible with the one from fermionic correlators andhas 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 η (cid:48) meson: except for the mixing withthe η , which can be taken into account by appropriately modifying the fit function, we do notexpect any additional complications.The topological susceptibility has been computed using the gradient flow. As expected, χ top is proportional to M π (for small M π ) up to an additive lattice artifact independent of M π . Evenif we are not able to finally confirm this with only a single lattice spacing at hand, this constantterm should be of O ( a ). The size of this artifact appears to be significantly smaller than what isobserved with Wilson clover fermions in Ref. [34].2 Acknowledgments
We thank the members of ETMC for the most enjoyable collaboration. We thank G. Rossi forvaluable comments on the manuscript. The computer time for this project was made available tous by the John von Neumann-Institute for Computing (NIC) on the Jureca and Juqueen systemsin J¨ulich and the HPC Cluster in Bern. This project was funded by the DFG as a project in theSino-German CRC110. U. W. acknowledges support from the Swiss National Science Foundation.The open source software packages tmLQCD [45–47], Lemon [48], DD α AMG [49], and R [50] havebeen used. [1] F. Jegerlehner, EPJ Web Conf. , 01016 (2016), arXiv:1511.04473 [hep-ph] .[2]
ETM
Collaboration, C. Michael, K. Ottnad and C. Urbach, Phys.Rev.Lett. , 181602 (2013), arXiv:1310.1207 [hep-lat] .[3]
ETM
Collaboration, K. Ottnad and C. Urbach, Phys. Rev.
D97 , 054508 (2018), arXiv:1710.07986[hep-lat] .[4] N. Christ et al. , Phys.Rev.Lett. , 241601 (2010), arXiv:1002.2999 [hep-lat] .[5] J. J. Dudek et al. , Phys.Rev.
D83 , 111502 (2011), arXiv:1102.4299 [hep-lat] .[6]
UKQCD
Collaboration, E. B. Gregory, A. C. Irving, C. M. Richards and C. McNeile, Phys.Rev.
D86 ,014504 (2012), arXiv:1112.4384 [hep-lat] .[7]
Hadron Spectrum
Collaboration, J. J. Dudek, R. G. Edwards, P. Guo and C. E. Thomas, Phys.Rev.
D88 , 094505 (2013), arXiv:1309.2608 [hep-lat] .[8]
JLQCD
Collaboration, H. Fukaya et al. , Phys. Rev.
D92 , 111501 (2015), arXiv:1509.00944[hep-lat] .[9] H. Leutwyler and A. Smilga, Phys. Rev.
D46 , 5607 (1992).[10]
JLQCD
Collaboration, S. Aoki, G. Cossu, H. Fukaya, S. Hashimoto and T. Kaneko, PTEP ,043B07 (2018), arXiv:1705.10906 [hep-lat] .[11] C. Alexandrou et al. , Phys. Rev.
D97 , 074503 (2018), arXiv:1709.06596 [hep-lat] .[12] N. J. Dowrick and N. A. McDougall, Phys. Lett.
B285 , 269 (1992).[13] G. D. Moore, EPJ Web Conf. , 01009 (2018), arXiv:1709.09466 [hep-ph] .[14] P. Di Vecchia, G. Rossi, G. Veneziano and S. Yankielowicz, JHEP , 104 (2017), arXiv:1709.00731[hep-th] .[15] UKQCD
Collaboration, C. McNeile and C. Michael, Phys.Lett.
B491 , 123 (2000), arXiv:hep-lat/0006020 [hep-lat] .[16] W. Sun et al. , Chin. Phys.
C42 , 093103 (2018), arXiv:1702.08174 [hep-lat] .[17]
ETM
Collaboration, K. Jansen, C. Michael and C. Urbach, Eur.Phys.J.
C58 , 261 (2008), arXiv:0804.3871 [hep-lat] .[18] M. L¨uscher, JHEP , 071 (2010), arXiv:1006.4518 [hep-lat] , [Erratum: JHEP03,092(2014)].[19] C. Alexandrou et al. , Phys. Rev. D98 , 054518 (2018), arXiv:1807.00495 [hep-lat] .[20]
ETM
Collaboration, C. Helmes et al. , EPJ Web Conf. , 05025 (2018), arXiv:1710.03698[hep-lat] .[21]
ETM
Collaboration, A. Abdel-Rehim et al. , Phys. Rev.
D95 , 094515 (2017), arXiv:1507.05068[hep-lat] .[22] Y. Iwasaki, Nucl. Phys.
B258 , 141 (1985).[23]
ALPHA
Collaboration, R. Frezzotti, P. A. Grassi, S. Sint and P. Weisz, JHEP , 058 (2001), hep-lat/0101001 .[24] B. Sheikholeslami and R. Wohlert, Nucl.Phys. B259 , 572 (1985).[25] T. Chiarappa et al. , Eur.Phys.J.
C50 , 373 (2007), arXiv:hep-lat/0606011 [hep-lat] .[26]
ETM
Collaboration, R. Baron et al. , JHEP , 111 (2010), arXiv:1004.5284 [hep-lat] .[27] R. Frezzotti and G. C. Rossi, JHEP , 007 (2004), hep-lat/0306014 .[28] Hadron Spectrum
Collaboration, M. Peardon et al. , Phys. Rev.
D80 , 054506 (2009), arXiv:0905.2160 [hep-lat] .[29] C. Morningstar et al. , Phys. Rev.
D83 , 114505 (2011), arXiv:1104.3870 [hep-lat] .[30]
ETM
Collaboration, C. Helmes et al. , JHEP , 109 (2015), arXiv:1506.00408 [hep-lat] .[31] E. V. Shuryak and J. J. M. Verbaarschot, Phys. Rev. D52 , 295 (1995), arXiv:hep-lat/9409020[hep-lat] .[32]
APE
Collaboration, M. Albanese et al. , Phys. Lett.
B192 , 163 (1987).[33] M. L¨uscher and P. Weisz, JHEP , 051 (2011), arXiv:1101.0963 [hep-th] .[34] ALPHA
Collaboration, M. Bruno, S. Schaefer and R. Sommer, JHEP , 150 (2014), arXiv:1406.5363 [hep-lat] .[35] L. Liu et al. , Phys. Rev. D96 , 054516 (2017), arXiv:1612.02061 [hep-lat] .[36] H. Neff, N. Eicker, T. Lippert, J. W. Negele and K. Schilling, Phys.Rev.
D64 , 114509 (2001), arXiv:hep-lat/0106016 [hep-lat] .[37] S. Aoki, H. Fukaya, S. Hashimoto and T. Onogi, Phys. Rev.
D76 , 054508 (2007), arXiv:0707.0396[hep-lat] .[38] G. S. Bali, S. Collins, S. D¨urr and I. Kanamori, Phys. Rev.
D91 , 014503 (2015), arXiv:1406.5449[hep-lat] .[39]
XLF
Collaboration, K. Jansen et al. , Phys. Lett.
B624 , 334 (2005), arXiv:hep-lat/0507032[hep-lat] .[40] P. Dimopoulos, R. Frezzotti, C. Michael, G. C. Rossi and C. Urbach, Phys. Rev.
D81 , 034509 (2010), arXiv:0908.0451 [hep-lat] .[41] G. Herdoiza, K. Jansen, C. Michael, K. Ottnad and C. Urbach, JHEP , 038 (2013), arXiv:1303.3516[hep-lat] . [42] UKQCD
Collaboration, C. R. Allton et al. , Phys. Rev.
D70 , 014501 (2004), hep-lat/0403007 .[43]
CP-PACS
Collaboration, V. I. Lesk et al. , Phys. Rev.
D67 , 074503 (2003), hep-lat/0211040 .[44] K. Hashimoto and T. Izubuchi, Prog. Theor. Phys. , 599 (2008), arXiv:0803.0186 [hep-lat] .[45] K. Jansen and C. Urbach, Comput.Phys.Commun. , 2717 (2009), arXiv:0905.3331 [hep-lat] .[46] A. Deuzeman, K. Jansen, B. Kostrzewa and C. Urbach, PoS
LATTICE2013 , 416 (2013), arXiv:1311.4521 [hep-lat] .[47] A. Abdel-Rehim et al. , PoS
LATTICE2013 , 414 (2014), arXiv:1311.5495 [hep-lat] .[48]
ETM
Collaboration, A. Deuzeman, S. Reker and C. Urbach, arXiv:1106.4177 [hep-lat] .[49] C. Alexandrou et al. , Phys. Rev.
D94 , 114509 (2016), arXiv:1610.02370 [hep-lat] .[50] R Development Core Team,