Ratio of flavour non-singlet and singlet scalar density renormalisation parameters in N_\mathrm{f}=3 QCD with Wilson quarks
Jochen Heitger, Fabian Joswig, Pia L. J. Petrak, Anastassios Vladikas
MMS-TP-21-02
Ratio of flavour non-singlet and singlet scalar densityrenormalisation parameters in N f = 3 QCDwith Wilson quarks
LPHA A Collaboration
Jochen Heitger a , Fabian Joswig a , Pia L. J. Petrak a , Anastassios Vladikas b a Westfälische Wilhelms-Universität Münster, Institut für Theoretische Physik,Wilhelm-Klemm-Straße 9, 48149 Münster, Germany b INFN, “Rome Tor Vergata” Division, c/o Dipartimento di Fisica,Via della Ricerca Scientifica 1, 00133 Rome, Italy
Abstract
We determine non-perturbatively the normalisation factor r m ≡ Z S /Z , where Z S and Z are the renormalisation parameters of the flavour non-singlet and singlet scalar densities,respectively. This quantity is required in the computation of quark masses with Wilsonfermions and for instance the renormalisation of nucleon matrix elements of scalar densities.Our calculation involves simulations of finite-volume lattice QCD with the tree-levelSymanzik-improved gauge action, N f = 3 mass-degenerate O( a ) improved Wilson fermionsand Schrödinger functional boundary conditions. The slope of the current quark mass, as afunction of the subtracted Wilson quark mass is extracted both in a unitary setup (wherenearly chiral valence and sea quark masses are degenerate) and in a non-unitary setup(where all valence flavours are chiral and the sea quark masses are small). These slopes arethen combined with Z ≡ Z P / ( Z S Z A ) in order to obtain r m . A novel chiral Ward identityis employed for the calculation of the normalisation factor Z . Our results cover the rangeof gauge couplings corresponding to lattice spacings below 0 . N f = 2 + 1QCD simulations in large volumes with the same lattice action are typically performed. Keywords:
Lattice QCD, quark mass renormalisation, Ward identities, Schrödinger functional,chiral Symmetry restoration with Wilson fermions
PACS: a r X i v : . [ h e p - l a t ] J a n Introduction
Scalar and pseudoscalar flavour singlet and non-singlet dimension-3 bilinear operators havethe same anomalous dimension, since they belong to the same chiral multiplet. The sameis true for their renormalisation parameters, provided that the regularisation does notbreak chiral symmetry. Otherwise, the renormalisation parameters of the chiral multipletcomponents differ by finite terms. This is the case for the lattice regularisation withWilson fermions. For example, the renormalisation parameters of the non-singlet scalarand pseudoscalar densities (denoted as Z S and Z P , respectively) have a finite ratio which isa polynomial of the bare gauge coupling g . This ratio can be determined by chiral Wardidentities; see Refs. [1,2]. Since Z P and Z S are scale dependent, imposing a renormalisationscheme is necessary to fix one of them, and the other can be obtained using the schemeindependent ratio Z S /Z P . In this way the renormalised scalar and pseudoscalar densitiesare defined consistently in the same scheme, with the same anomalous dimension andrenormalisation group (RG) running, and chiral symmetry is restored in the continuumlimit. The ratio Z S /Z P has been computed for several gauge and Wilson fermion actions(standard, improved etc.) in the quenched approximation [2, 7–11], with two dynamicalquarks ( N f = 2 QCD) [12], and with three dynamical quarks ( N f = 3 QCD) [13–16].Far less progress has been made on the computation of the ratio of the renormalisationparameters of the non-singlet and singlet scalar densities, r m ≡ Z S /Z . For chirallysymmetric regularisations r m = 1 holds, while for Wilson fermions r m is a (finite) polynomialof the gauge coupling, arising from the sea fermion loops of the quark propagator. In thequenched approximation, r m = 1. As explained in Ref. [17], the lowest-order non-trivialperturbative contribution to this quantity is a two-loop effect; i.e., r m = 1 + O( g ). InRef. [18] the O( g ) perturbative term has been calculated for several lattice actions. Non-perturbative estimates of this quantity have been reported in Ref. [13] at two values ofthe gauge coupling for N f = 2 + 1 QCD with the tree-level Symanzik-improved gaugeaction [19] and the non-perturbatively improved Wilson-clover fermion action [20]. This isthe regularisation chosen by the CLS ( Coordinated Lattice Simulations ) initiative whichcarries out QCD simulations with N f = 2 + 1 flavours, on large physical volumes, for a rangeof bare couplings corresponding to a hadronic regime [13, 21–23]. These CLS ensemblesare suitable for the computation of correlation functions, from which low-energy hadronicquantities can be evaluated. In parallel, our group is performing N f = 3 simulations inthe same range of bare gauge couplings, but for small-volume lattices with Schrödingerfunctional boundary conditions and nearly-chiral quark masses. These ensembles areused for the numerical determination of the necessary renormalisation parameters andSymanzik improvement coefficients, see Refs. [14,15,24–28] that have various applications inlattice QCD when using this discretisation of Wilson fermions. The present work provideshigh-precision estimates of r m obtained in the same computational framework.As seen from eq. (2.2) below, ( r m −
1) contributes an O( g ) term to the renormalisationof the quark masses [17]. This is expected to be a small effect. Symanzik O( a ) countertermscontaining r m are often neglected in light quark mass determinations; cf. Ref. [29]. In In practice, distinct chiral Ward identities are used for the computation of the ratio Z S / ( Z P Z A ) and Z A ; the two results are subsequently multiplied to give Z S /Z P . Examples of renormalisation schemes are MS, RI/(S)MOM [3, 4], the Schrödinger functional (SF) [5]and the chirally rotated Schrödinger functional ( χ SF) [6]. r m can be relevant at O( a ), especially when dealing withheavy flavours, and should be taken into account in order to achieve full O( a ) improvement;see, for example, eq. (2.13) in Ref. [30]. Another application where r m plays a prominentrôle is the nucleon sigma-term, which is defined in terms of nucleon matrix elements offlavour singlet scalar densities; see Refs. [31, 32] for example and [33–35] for more recentworks. A direct determination of Z is not as straightforward as that of Z S , the formeralso requiring the computation of two-boundary (“disconnected”) quark diagrams. Thisproblem is circumvented by extracting Z as the product of Z S and r m .Our computation of r m is based on the relation between the current (PCAC) mass m and the subtracted quark mass m q . Close to the chiral limit, m ( m q ) is a linear functionwith a slope that depends on the details of the QCD model being simulated. In a unitarytheory with degenerate sea and valence quark masses, the slope of m ( m q ) is Zr m , where Z ≡ Z P / ( Z S Z A ) and Z A is the non-singlet axial current normalisation. On the other hand,in a non-unitary theory with chiral valence subtracted quark masses ( m valq = 0) and smalldegenerate sea quark masses m seaq = 0, the slope of m ( m seaq ) is Z ( r m − g . The results are combined to give estimates of r m ( g ). This approach isdescribed in Section 2.Alternatively, each of the two slopes Zr m and Z ( r m −
1) may be combined with anindependent estimate of Z , such as the results of Refs. [14, 15]. In the present work weprefer to use a novel determination of Z , relying on a chiral Ward identity which differsfrom the one of Ref. [15]. This identity is derived in Section 3.In Section 4 we present our simulation setup for N f = 3 QCD with lattices of smallphysical volumes and Schrödinger functional boundary conditions; these serve to numericallyimplement the strategies outlined in the foregoing section. Most of our gauge field ensembleswere already generated in the context of previous works; cf. Refs. [14, 15, 25–28]. Somenew ensembles have also been generated, in order to cover the region close to the origin ofthe function m ( m q ) more evenly and asses its slope reliably.Our results for r m , based on various combinations of Zr m , Z ( r m − Z arediscussed in Section 5. Different determinations of r m are compared, allowing us to settlefor a conservative final estimate with reliable systematic errors. Our final result is that ofeq. (5.5). In Table 5 we also list r m ( g ) for the g -values at which CLS simulations arebeing performed for the computations of hadronic quantities in N f = 2 + 1 QCD.In the final section we sum up our results and their uses in lattice QCD. Moredetailed calculations and definitions of the correlation functions employed can be found inAppendix A and B. Comparison of Z determinations and corresponding scaling tests canbe found in Appendix C. In this section we recapitulate the basic quark mass definitions, namely subtracted andcurrent (PCAC) quark masses, and discuss how to obtain the products Zr m and Z ( r m − i = 1 , . . . , N f , m q ,i ≡ m ,i − m crit = 12 a (cid:16) κ i − κ crit (cid:17) , (2.1)3here κ i is the hopping parameter for flavour i , κ crit its value in the chiral limit, and a isthe lattice spacing. In terms of the subtracted masses m q ,i , the corresponding renormalisedquark masses are given by m i, R = Z m " m q ,i + ( r m −
1) Tr M q N f + O( a ) , (2.2)where M q = diag( m q , , . . . , m q ,N f ) is the N f × N f bare quark mass matrix.We recall in passing that the renormalisation parameter Z m ( g , aµ ) depends on therenormalisation scale µ and diverges logarithmically in the ultraviolet. It is the inverse of Z S ( g , aµ ), the renormalisation parameter of the flavour non-singlet scalar density operator.A mass independent renormalisation scheme is implied throughout this work. In such ascheme operator renormalisation parameters (e.g. Z P , Z m , Z S ), current normalisations (i.e. Z A , Z V ) and r m are functions of the squared bare gauge coupling g . In a non-perturbativedetermination at non-zero quark mass, they are affected by O( am q ,i ), O( a Tr M q ), andO( a Λ QCD ) discretisation effects, which are part of their operational definition. As pointedout in Ref. [17], the term ( r m −
1) multiplies Tr M q , so it arises from a mass insertionin a quark loop. In perturbation theory it is a two-loop effect, contributing at O( g ).Its non-perturbative determination is the main purpose of this paper. An importantconsequence of eq. (2.2) is that a renormalised mass m i, R goes to the chiral limit only when all subtracted masses m q , , . . . , m q ,N f vanish.Alternatively, a bare current (PCAC) quark mass m ij can be defined through thefollowing relation: ( e ∂ µ ) x (cid:10) ( A I ) ijµ ( x ) O ji (cid:11) = 2 m ij (cid:10) P ij ( x ) O ji (cid:11) . (2.3)The quantity m ij is distinct from the subtracted bare quark masses, but it is related tothe mass average ( m q ,i + m q ,j ) /
2; see eq. (2.9) below. The flavour non-singlet bare axialcurrent and the pseudoscalar density are given by A ijµ ( x ) ≡ ¯ ψ i ( x ) γ µ γ ψ j ( x ) , P ij ( x ) ≡ ¯ ψ i ( x ) γ ψ j ( x ) , (2.4)with indices i, j denoting two distinct flavours ( i = j ). The pseudoscalar density P ij andthe current ( A I ) ijµ ≡ A ijµ + ac A e ∂ µ P ij are Symanzik-improved in the chiral limit, with theimprovement coefficient c A ( g ) being in principle only a function of the gauge coupling. Inthese definitions, e ∂ µ denotes the average of the usual forward and backward derivatives. The source operator O ji is defined in a region of space-time that does not include thepoint x , so as to avoid contact terms. In the O( a ) improved theory, the renormalised axialcurrent and pseudoscalar density are( A R ) ij ( x ) = Z A ( g )( A I ) ijµ ( x ) + O( am q , a ) , (2.5)( P R ) ij ( x ) = Z P ( g , aµ ) P ij ( x ) + O( am q , a ) . (2.6)The normalisation of the axial current Z A ( g ) is scale independent, depending only on thesquared gauge coupling g . The renormalisation parameter Z P ( g , aµ ) (determined, say inthe Schrödiger functional scheme of Ref. [5]) additionally depends on the renormalisation The forward derivative is defined as a∂ µ f ( x ) ≡ f ( x + a ˆ µ ) − f ( x ) and the backward derivative as a∂ ∗ µ f ( x ) ≡ f ( x ) − f ( x − a ˆ µ ). µ and diverges logarithmically in the ultraviolet. The PCAC relation, expressed byrenormalised fields,( e ∂ µ ) x (cid:10) ( A R ) ijµ ( x ) O ji (cid:11) = ( m R ,i + m R ,j ) (cid:10) ( P R ) ij ( x ) O ji (cid:11) , (2.7)valid up to discretisation effects in the continuum, combined with eqs. (2.3)–(2.6), impliesthat m i R + m j R Z A Z P m ij + O( am q , a ) . (2.8)If we calculate the average mass ( m i R + m j R ) / m ij = Z " ( m q ,i + m q ,j )2 + ( r m −
1) Tr M q N f + O( am q , a ) , (2.9)where the product of the renormalisation parameters Z ( g ) ≡ Z P ( g , µ ) / ( Z S ( g , µ ) Z A ( g ))is scale independent. We now exploit eq. (2.9) in two ways:( ) In a theory with mass-degenerate quarks ( m q ,i = m q ,j = Tr M q /N f ), it reduces to m = Zr m m q + O( am q , a ) (2.10)= Zr m a (cid:16) κ − κ crit (cid:17) + O( am q , a ) . (2.11)In the above equation, flavour indices have been dropped from the quark masses m ij , m q ,i and the hopping parameter κ i . This simplification of notation will be adopted on mostoccasions below. Thus, modelling the current quark mass am as a function of 1 /κ forvalues of κ close to κ crit , we obtain the latter as the root of the function am (1 /κ ) and thecombination Zr m as the slope of the same curve.( ) Once the critical hopping parameter κ crit is available from the previous step ( ),we use a non-unitary setup where valence and sea quarks of the same flavour have differentbare subtracted masses m valq ,i = m seaq ,i . In eq. (2.9), masses m q ,i and m q ,j on the r.h.s.are valence quark contributions, while Tr M q stands for the trace of sea quark masses;see Refs. [14, 17] for detailed explanations. In particular, we set κ val = κ crit , so as toensure m valq = 0 for all valence flavours. Moreover sea quark masses are taken to be small,degenerate, and non-zero (i.e. κ sea = κ crit , ensuring m seaq = 0 for all sea flavours). Withthese conditions, the current quark mass of eq. (2.9) reduces to m = Z ( r m − m seaq + O( am q , a ) . (2.12)It is remarkable that with non-zero bare subtracted sea quark masses (i.e. m seaq = 0), allcurrent quark masses in this setup are not chiral (i.e. m ij = 0 , ∀ i, j ), even if all subtractedvalence quark masses vanish (i.e. m valq ,i = 0 , ∀ i ). From eq. (2.12) we see that, if we compute am as a function of am seaq for several sea masses, the slope of the functions gives an estimateof Z ( r m − Zr m and Z ( r m − g , can be combined yielding estimates of r m ( g );see Subsection 5.3 for details. We stress that the above discussion concerns relations whichsuffer from O( a ) discretisation effects. For the quark masses, such effects may be removed5y introducing Symanzik counterterms, leaving us with O( a ) discretisation errors. Thesecounterterms have been worked out in Refs. [17, 36]. In Ref [17] (see also eq. (2.10) ofRef. [14]) the full O( am q ) contributions, omitted in eq. (2.9) above, are written downexplicitly. Such contributions are complicated and taking them all into account couldcompromise the numerical stability of our procedure to extract the quantities in question.We prefer a simpler and more robust strategy, consisting of working with small quarkmasses so that O( a )-effects in eq. (2.9) may be safely dropped. This must of course bechecked a posteriori , by ensuring that the function m ( m q ) is linear close to the origin,where our simulations are performed. The only improvement coefficients used in this workare c sw of the clover action and c A , of the axial current (entering the PCAC mass). Z In the previous section we have shown how the quantities Zr m and Z ( r m −
1) can beestimated from relations between suitably chosen current and subtracted Wilson quarkmasses. They may then straightforwardly be combined to give r m and Z . The latterquantity has already been measured in our setup ( N f = 3 lattice QCD with Schrödingerfunctional boundary conditions) in two ways: either by using appropriate combinationsof current and subtracted quark masses with different flavours [14], or from chiral Wardidentities [15] via Z ≡ Z P / ( Z A Z S ). Here we will describe yet another direct method, basedon a new Ward identity, very similar to the one of Ref. [15]. The reader is referred to thatwork for details, notation etc.We consider a product of two composite operators O ≡ S b ( y ) O c , defined as S b ( y ) ≡ i ¯ ψ ( y ) T b ψ ( y ) O c ≡ i a L X u , v ¯ ζ ( u ) γ T c ζ ( v ) , (3.1)where T b and T c are generators of SU ( N f ). The former operator is the flavour non-singletscalar density, located in the bulk of space-time, while the latter resides at the x = 0Dirichlet time boundary of the Schrödinger functional. The Ward identity of interest isobtained by performing axial variations on O in a region R , chosen to be the space-timevolume between the hyper-planes at t and t where t < t . With O c lying outside R , wehave δ A O = [ δ A S b ( y )] O c and δ A S b ( x ) = (cid:15) a h d abe P e ( x ) + δ ab N f ¯ ψ ( x ) ψ ( x ) i . (3.2)In what follows we simplify matters by always working with a = b , so as to eliminate thesecond contribution on the r.h.s. of the above expression. In analogy to the derivation For reasons of convenience, we have adopted a slightly different notation in this section: the flavourcontent of operators like S b or O c is determined by a single flavour index b or c , corresponding to its flavourmatrix T b or T c . The fermion fields of these operators ψ and ¯ ψ are columns in flavour space. This is to becontrasted to the notation of Section 2, where we have introduced operators like P ij and O ji , which haveexplicit indices, referring to the flavour of fields ψ j , ¯ ψ i etc. Z d y Z d x Dh A a ( t ; x ) − A a ( t ; x ) i S b ( y ; y ) O c E − m Z d y Z d x Z t t d x h P a ( x ; x ) S b ( y ; y ) O c i (3.3)= − d abe Z d y h P e ( y ) O c i . Next we adapt the previous formal manipulations to the lattice regularisation withSchrödinger functional boundary conditions. The pseudoscalar operator O c is defined onthe x = 0 time boundary. Ward identity (3.3) then becomes: Z A Z S a ( X x , y Dh ( A I ) a ( t ; x ) − ( A I ) a ( t ; x ) i S b ( y ; y ) O c E − am X x , y t X x = t w ( x ) h P a ( x ; x ) S b ( y ; y ) O c i ) (3.4)= − d abe Z P a X y h P e ( y ) O c i + O( am, a ) . In this expression repeated flavour indices e are summed, as usual. The weight factoris w ( x ) = 1 / x ∈ { t , t } and w ( x ) = 1 otherwise. It is introduced in order toimplement the trapezoidal rule for discretising integrals. Quark masses are degenerate and m is the current quark mass.The last step is to perform the Wick contractions in Ward identity (3.4). How this isdone is explained in Appendix B; eventually, flavour factors drop out and we are left with aWard identity that translates into traces of products of quark propagators and γ -matrices,graphically depicted in Fig. 1. Solving for Z we get Z ≡ Z P Z A Z S = − f IAS ( t , y ) − f IAS ( t , y ) − am ˜ f PS ( t , t , y ) f P ( y ) + O( am, a ) , (3.5)where dependencies are suppressed on the l.h.s. Assuming that we work in the chiral limit(or with nearly-vanishing quark masses, so that O( am ) effects may be safely neglected),the above Ward identity is valid up to O( a ) discretisation errors in lattice QCD withWilson quarks. In this spirit, terms proportional to Symanzik b -coefficients may also besafely ignored. The renormalisation factor of the external source O c is not taken intoconsideration, as it cancels out in the ratio (3.5). The term proportional to the currentquark mass m may also be dropped close to the chiral limit, but since we are working withmasses which are not strictly zero, it could be advantageous to keep it in practice. In fact,it was found in Refs. [15, 26] that this term stabilizes the chiral extrapolation leading tosmaller errors. This turns out to be true also in our case, as we will show in Subsection 5.2and Fig. 5.It is interesting to compare Ward identity (3.3) with those of Ref. [15]: This is even true for light (up/down, strange) non-chiral quark masses, as explicitly demonstrated inRef. [29], using the b -coefficients of Ref. [14]. a) Diagram f P (b) Diagram f AS;1 (c)
Diagram f AS;2
Figure 1:
The trace diagrams contributing to the expectation values of f P , defined in eq. (B.1)(diagram (a)) and f AS , defined in eq. (B.3) (diagrams (b) and (c)). The wall represents the timeslice x = 0 with a γ Dirac matrix between circles. The squares in the bulk represent either theinsertions of a pseudoscalar operator P ( y ) (diagram (a)) or a scalar operator S ( y ) (diagrams (b)and (c)). The diamonds stand for an axial operator A ( x ). The open circles correspond to theboundary fields ζ , while the filled circles denote ¯ ζ . The diagrams schematically represent traces,formed by starting from any point and following the lines (quark propagators) until we close theloop. The time ordering of points x and y is left unspecified in these diagrams. • In Ref. [15] the flavour factors gave rise to a multitude of identities, which werecombined in order to increase the signal-to-noise ratio, while here we only have oneidentity. On these grounds one could expect that the numerical results of Ref. [15]are more precise than the ones from the Ward identity introduced here.• On the other hand, the identities of Ref. [15] involved: (i) correlation functions withone operator insertion in the bulk of the lattice and one wall source at each timeslice; cf. Fig. 1 in that work; (ii) correlation functions with two operator insertionsin the bulk and one wall source at each time slice; cf. Fig. 2 in that work. Here wehave: (i) a correlation function with one operator insertion in the bulk and one wallsource; (ii) correlation functions with two operator insertions in the bulk and onewall source. These somewhat simpler correlation functions illustrated in Fig. 1 aboveare expected to have less statistical fluctuations. From this point of view, the resultsof the present work are expected to gain in accuracy.Thus, one of our aims is to establish which of the two approaches leads to more accurateresults. This is discussed in Subsection 5.2 and Appendix C. We employ the tree-level Symanzik-improved gauge action and N f = 3 mass-degenerate O( a )improved Wilson fermions. For the corresponding improvement coefficient c sw we use thenon-perturbative determination of Ref. [37]. As already indicated, we impose Schrödingerfunctional boundary conditions at the temporal boundaries of the lattice. For the boundaryimprovement coefficients we employ the respective tree-level values c t = ˜ c t = 1. TheSchrödinger functional setup is highly suitable for massless renormalisation schemes, sincenearly-vanishing quark masses are accessible in numerical calculations due to the spectralgap of the Dirac operator. This gap is imposed by the boundaries, so that the quark massdependence can be mapped out reliably in the vicinity of the chiral point. The generation We explicitly checked that already at the coarsest lattice spacing, a noticeable variation of the boundaryimprovement coefficients has a negligible effect on the PCAC mass plateaux. L/a ) × T /a β κ a (in fm)12 ×
17 3.3 0.13652 20 20480 A1k1 0 . ×
18 0 . A3k1 . A3k2 . A3k3 . A3k4 . A3k5 . A3k6 ×
21 3.414 0.13690 32 38400 E1k1 0 . ×
20 0.13656 18 60480
E2k1
E2k2 ×
23 3.512 0.13700 2 20480 B1k1 0 . ×
24 0.13677 1 25904
B3k1 ×
29 3.676 0.13680 1 7848 C1k1 0 . ×
35 3.810 0.13711875582 5 8416 D1k1 ∗ . Table 1:
Simulation parameters L , T , β , κ , the number of replica italics were newly generated for this study while the remaining ones were already used in previousinvestigations (see, for example Ref. [14]). The ensemble D1k1 marked by an asterisk is only usedfor the determination of the PCAC masses. The lattice spacings a are obtained by interpolatingthe results of Ref. [22] with a polynomial fit. All configurations are separated by 8 MDU’s exceptfor the ensembles A1k3 (4 MDU’s) and D1k4 (16 MDU’s). of the gauge field configurations is performed with the openQCD code [38] which employsthe RHMC algorithm [39, 40] for the third quark.9ll gauge field ensembles used in this study are summarized in Table 1 and lie on aline of constant physics, defined by a fixed spatial extent of L ≈ . T /L ≈ / r m and Z become smooth functions of the lattice spacing, with relevanthigher-order ambiguities vanishing monotonically. The gauge ensembles highlighted in italics were newly generated for this study, while the remaining ones were already usedin previous investigations; see Refs. [14, 15, 24–28]. These additional ensembles allowfor a more even and wider spread of bare quark masses around the chiral point for eachvalue of β , which enables a more precise extraction of the slopes corresponding to Zr m and Z ( r m −
1) as explained in Section 2. Since a newer version of the openQCD code wasutilised for the generation of the ensembles, the time extent T , which was odd in thepre-existing ensembles, is even for the new ones. This slight violation of the condition of aline of constant physics in T only affects the current quark masses at the level of negligiblehigher-order ambiguities, as demonstrated in Ref. [27].All Schrödinger functional correlation functions required for our numerical investiga-tions are O( a ) improved. In this context we only require the improvement coefficient c A ,non-perturbatively known from Ref. [26]. Since the Markov chain Monte Carlo sampling ofthe gauge field configurations suffers from critical slowing down of the topological charge forsmaller lattice spacings (see Ref. [42]), we project our data to the trivial topological sectoras suggested in Ref. [43], in order to account for the insufficient sampling of all topologicalsectors. For the analysis of the statistical errors we employ the Γ-method [44]. We accountfor the remaining critical slowing down of the Monte Carlo algorithm by attaching a tail tothe autocorrelation function, as suggested in Ref. [45]. The corresponding slowest mode isestimated from the autocorrelation time of the boundary-to-boundary correlation function F ij , defined in Appendix A. The error analysis is carried out with a python implementationof the Γ-method, using automatic differentiation for the error propagation as proposed inRef. [46]. In the following we present our analysis which eventually leads to several estimates forthe ratio of the renormalisation parameters of the non-singlet and singlet scalar densities, r m . We will first describe how we obtain Zr m , Z ( r m − Z individually and thendiscuss several ways of combining the three into r m . As a final result we provide aninterpolation formula for r m and extract its value at the bare couplings of large-volume CLS simulations [13, 21, 23].
As described in Section 2, the quantities Zr m and Z ( r m −
1) can be extracted from quarkmass slopes. Our results are based on the determination of the O( a ) improved PCAC In a setup with heavy sea quarks and very light valence quarks we approach a quenched-like situation inwhich exceptional configurations are to be expected; cf. Ref. [41] where a similar situation is discussed. In acareful analysis we identified only one gauge field configuration in the ensemble E2k1, with an exceptionallysmall eigenvalue of the massless Dirac operator. This leads to very large values of the correlation functions f P and f A . We have discarded this exceptional configuration. x /a − . . . . . a m B3k1 κ val = κ sea κ val = κ crit x /a − . − . . . . . a m D1k2 κ val = κ sea κ val = κ crit Figure 2:
PCAC mass am as a function of time x /a for ensembles B3k1 (left) and D1k2 (right).Squares are results obtained in the unitary setup, while diamonds are results obtained in thenon-unitary setup. The final estimate for m is obtained by averaging results in the time interval[ T / , T / masses via m ( x ) = e ∂ f ij A ( x ) + ac A ∂ ∗ ∂ f ij P ( x )2 f ij P ( x ) , (5.1)where f ij A and f ij P are Schrödinger functional correlation functions. In order to improve thesignal, these correlation functions are symmetrised with their T -symmetric counterparts g ij A ( T − x ) and g ij P ( T − x ), which are constructed from the same operators ( A I ) ij ( x ) and P ij ( x ) in the bulk but the pseudoscalar wall with operator O ji positioned at the timeboundary x = T . For exact definitions see Appendix A.We first determine the required correlation functions in a unitary setup, κ val = κ sea .From these we can obtain κ crit as will be detailed below. In a second step we compute thesame correlation functions in a non-unitary setup where κ sea = κ val = κ crit . In Fig. 2 weshow the temporal dependence of the current quark mass m ( x ) for both of these setups forthe representative ensembles B3k1 and D1k2 and demonstrate that they form well-definedplateaux as a function of time, away from the Dirichlet boundaries. Our final estimate forthe PCAC masses is obtained by averaging m ( x ) over the central third of the temporalextent of the lattice. This choice is motivated by the coarsest lattices; the plateaux for thefiner ones also extend closer to the boundary before lattice artefacts become relevant ascan be seen in Fig. 2. The plateau range is adapted according to the time extent for eachvalue of β , so as to preserve the line of constant physics. Our PCAC mass estimates inboth setups are listed in Table 2 for all ensembles.In order to extract Zr m and Z ( r m −
1) from the slopes of the current quark masses withrespect to the bare quark masses, we plot am against the inverse hopping parameter 1 / (2 κ )for both the unitary and the non-unitary setup, as demonstrated in Fig. 3. We generallyobserve that m behaves linearly as a function of 1 / (2 κ ) in the range − . (cid:46) Lm (cid:46) . am Z { T/ } Z { T/ } κ val i = κ sea i κ val i = κ crit A3k3 0 . . . . . . . . . . . . − . − . − . − . − . − . . . . . . . . . − . − . . . . . . . . . . . − . − . . . . . . . − . − . . . . . − . − . − . − . . . Table 2:
For each ensemble, identified in the first column by an ID label, we list our results for thePCAC mass am for simulations with κ val = κ sea (second column) and κ sea = κ val = κ crit (thirdcolumn). The last two columns contain Z results obtained from the Ward identity (3.5). The finalresults are those extrapolated to the chiral limit at each β = 6 /g (last line of each data grouping).The labels Z { T/ } and Z { T/ } refer to different choices of time slices with operator insertions inthe correlation functions (see text for details). .
645 3 .
650 3 .
655 3 .
660 3 .
665 3 . / (2 κ ) − . − . . . . a m β = 3.3 χ / d . o . f . = 0.71 χ / d . o . f . = 0.84 κ val = κ sea κ val = κ crit .
645 3 .
650 3 .
655 3 .
660 3 .
665 3 . / (2 κ ) − . − . . . . a m β = 3.414 χ / d . o . f . = 0.4 χ / d . o . f . = 0.85 κ val = κ sea κ val = κ crit .
645 3 .
650 3 .
655 3 .
660 3 .
665 3 . / (2 κ ) − . − . . . . a m β = 3.512 χ / d . o . f . = 1.43 χ / d . o . f . = 1.72 κ val = κ sea κ val = κ crit .
645 3 .
650 3 .
655 3 .
660 3 .
665 3 . / (2 κ ) − . − . . . . a m β = 3.676 χ / d . o . f . = 0.89 χ / d . o . f . = 0.81 κ val = κ sea κ val = κ crit .
645 3 .
650 3 .
655 3 .
660 3 .
665 3 . / (2 κ ) − . − . . . . a m β = 3.81 χ / d . o . f . = 0.03 χ / d . o . f . = 0.03 κ val = κ sea κ val = κ crit Figure 3:
PCAC masses am fitted linearly in 1 / (2 κ ), for all simulated β values (i.e. for decreasinglattice spacings from top to bottom). Open squares and filled diamonds are results in the unitaryand non-unitary setups, respectively. Note that horizontal and vertical axes are identical for allvalues of β , so as to highlight the different ranges of κ and the change of κ crit marked by the verticaldashed lines. . . . . . . . . / (2 κ ) . . . . . . a m χ / d . o . f . = 1.1, Zr m = 4.21(10) χ / d . o . f . = 1.86, Z ( r m −
1) = 3.64(10) κ val = κ sea κ val = κ crit Figure 4:
PCAC masses am , at β = 3 .
3, fitted with a quadratic polynomial in 1 / (2 κ ). Squaresand diamonds are results in the unitary and non-unitary setups, respectively. Note that the tworightmost points (A3k1, A3k3) are not included in the fit, while A3k2 (at 1 / (2 κ ) = 3 . κ crit from the linear fit; see Table 3. The Zr m and Z ( r m − For the ensembles A3k1, A3k2, and A3k3 (not displayed in Fig. 3), which correspond to Lm (cid:38) .
3, linearity is lost. Results from these ensembles have thus not been included inthe linear fits. The good linear behaviour of the data from the remaining ensembles isjustified a posteriori , by the small χ / d . o . f. of our fits, as shown in Fig. 3.We also probe the non-linear regime in both setups for β = 3 . am q ) ) effects in this case. Thetwo rightmost points (A3k1, A3k3) have not been included in these fits. Including themwould result in a very large value of χ / d . o . f. This may also be related to the fact thatno clear-cut plateaux are seen in the current quark mass data for these ensembles. Thiscould be explained by the fact that (boundary) cut-off effects for these comparatively largemasses (in lattice units) are substantial. Estimates of Zr m and Z ( r m − κ crit , are compatible with those from linearfits. The influence of the quadratic term on our final result is therefore negligible. Thisensures that our results are not affected by O(( am q ) ) systematic errors at β = 3 .
3, whichis our coarsest lattice. The same conclusion holds for the finer lattices, since also for them am q is small and linear fits have small χ / d . o . f . As implied by eq. (2.11), κ crit and Zr m are assessed as the intercept and the slopeof the linear fit to the unitary data. Similarly, eq. (2.12) tells us that Z ( r m −
1) can beestimated from the slope of the linear fit to the non-unitary data. Our final findings for Zr m , Z ( r m − κ crit are listed in Table 3. Z As the next step in our analysis, we extract the renormalisation constant Z ≡ ( Z P /Z S Z A )from the ratio (3.5), using the subset of gauge field ensembles listed in Table 1 which are14 Zr m Z ( r m − κ crit Table 3:
Results from the PCAC mass analyses. The second and fourth column show resultsobtained in a unitary setup; the third column refers to the non-unitary setup. not emphasised in italic font. The correlation functions in eq. (3.5) are computed for twochoices of t and t . Our first choice is t ≈ T / t ≈ T /
3, and the results obtainedin this fashion are denoted as Z { T/ } . Alternatively, choosing t ≈ T / t ≈ T / Z estimate denoted as Z { T/ } . When T /
T / t and t are rounded up/down to the nearest integer.In the left part of Fig. 5, we depict Z { T/ } as a function of y /T for several representativeensembles (we remind the reader that T is approximately constant in physical units).Contrary to the PCAC masses in Fig. 2, these local estimators of Z do not exhibit plateau-like behaviour; this was also observed for a similar Ward identity adopted to compute theimprovement coefficient of the vector current in Ref. [25]. Note, however, that this is notproblematic; since Z is obtained from a Ward identity, its value at any time slice qualifiesas a well-defined estimate. We prefer to err on the side of caution and quote the averageof the two central time slices as our best Z estimate. Results for the two determinationsof Z are collected in Table 2, where we see that Z { T/ } and Z { T/ } are not compatible,indicating the presence of lattice artefacts that also differ noticeably. We consider Z { T/ } the more reliable estimate because the operator insertions in this case, being further from x = 0 and x = T , are expected to lead to less contamination through cut-off effectsinduced by the boundaries.Since the Ward identity (3.5) is only valid up to lattice artefacts of O( am, a ), we haveto interpolate our data to the chiral point, in order to eliminate the O( am )-effects and beleft with O( a ) only. As an additional cross-check we also compute Z without the “massterm” 2 am ˜ f PS ( t , t ) in the Ward identity (3.5), where am is the PCAC mass from theunitary setup discussed in the previous section. This chiral interpolation is demonstratedfor β = 3 .
676 in the right part of Fig. 5. While the data including the “mass term” showsa very flat behaviour with respect to the current quark mass (where the associated fitparameter even vanishes within its uncertainty except for the coarsest lattice spacing), thetruncated Ward identity results in a considerably larger slope. If we exclude the rightmostdata point for the identity without the “mass term”, linear fits to both datasets still agreein the chiral limit. This situation resembles closely what was observed in Ref. [15], where Z was measured employing a different Ward identity. We note that the linear fit is based As explained in Section 4, the ensembles in italics have been generated for the purpose of performingreliable fits of the data in Fig. 3 and 4, in order to accurately measure their slopes. These extra ensembleshave not been used for the computation of Z , as they do not increase the accuracy of the result. D1k1(marked by an asterisk) is also not taken into account. . . . . . y /T . . . . . . Z { T / } D1k4C1k3 B1k4E1k2 . . . . . . am . . . . . β = 3.676 Z { T/ } massive Z { T/ } massless Figure 5:
Left:
Ward identity estimates of Z , plotted against time y /T , for one representativeensemble for each lattice spacing (except for β = 3 .
3, corresponding to the coarsest lattice). Thedashed vertical lines bracket the two central time slices that determine the final value of Z . Right:
Chiral extrapolation of Z at fixed β obtained from the Ward identity with the mass term (squares)and without it (diamonds). In the massless case, a possible linear range in am is illustrated bythe dashed line joining the two leftmost points. In the massive case, no significant quark massdependence is observed; the dashed line through the squares is a linear fit where the slope vanisheswithin its uncertainty. Note that the errors of the PCAC masses are also displayed and taken intoaccount in the fits via orthogonal distance regression [47]. r { u , nu } m r { u; Z,T/ } m r { nu; Z,T/ } m r { u; Z,T/ } m r { nu; Z,T/ } m Table 4:
Results for r m , obtained via eqs. (5.2) to (5.4). on the orthogonal distance regression method [47], taking into account both the error ofdependent and independent variables. The final results for Z { T/ } and Z { T/ } at the chiralpoint are also listed in Table 2. Compared to the indirect Ward identity determination ofRef. [15], they have considerably smaller errors. This confirms the expectation that thesimpler structure of the correlation functions building the Ward identity (3.4) is preferablefrom a numerical perspective; see the discussion at the end of Section 3. On the otherhand, compared to the so-called ’LCP-0’ determination of Ref. [14], our results are ofsimilar accuracy across the bare couplings investigated. We will use our results (Table 2)for a precise estimation of r m in the following. More details on the relative cut-off effectsbetween the present determination of Z and the results obtained in Refs. [13–15] can befound in Appendix C. r m In the final step of our analysis we combine the values of Zr m obtained in a unitary setup, Z ( r m −
1) in a non-unitary setup, and Z from a chiral Ward identity, in order to arrive atdifferent estimates for r m . Combining the first two, we construct r { u , nu } m , defined as r { u , nu } m = (cid:18) − (cid:20) Z ( r m − Zr m (cid:21) (cid:19) − , (5.2)where the superscripts “u” and “nu” stand for “unitary” and “non-unitary”, respectively.Combining Zr m and Z , results in r { u; Z } m , defined as r { u; Z } m = Zr m Z . (5.3)As mentioned above, this comes in two versions, r { u; Z,T/ } m and r { u; Z,T/ } m . Moreover, fromthe second and third result we gain r { nu; Z } m given by r { nu; Z } m = Z ( r m − Z + 1 , (5.4)which is again worked out for two cases, r { nu; Z,T/ } m and r { nu; Z,T/ } m . All our results for r m from these different determinations just outlined are gathered in Table 4.In principle, the different estimates can differ by O( a ) ambiguities. In Fig. 6 (left)the three determinations r { u , nu } m , r { u; Z,T/ } m , and r { nu; Z,T/ } m are plotted against the bare17 .
55 1 .
60 1 .
65 1 .
70 1 . g . . . . . . r { u , nu } m r { u;Z , T / } m r { nu;Z , T / } m .
000 0 .
002 0 .
004 0 .
006 0 .
008 0 .
010 0 . a in fm . . . . . . . . r { u , nu } m / r { nu;Z , T / } m r { nu;Z , T / } m / r { u;Z , T / } m r { u , nu } m /r { u , nu } m impr Figure 6:
Left:
Results for different r m estimates as reported in Table 4. The results for β = 3 . Right:
Ratio of different r m determinations as a function of the squared latticespacing. The dashed horizontal line indicates the expected continuum result. coupling squared; to be able to distinguish between the different estimates, the data pointscorresponding to the coarsest lattice spacing ( β = 3 .
3) are omitted as they exhibit largecut-off effects and are thus well out of the range displayed here. Results are compatiblewithin their respective 1 σ -errors. In Fig. 6 (right) we take a closer look at this behaviourby plotting ratios of different r m estimates as functions of the lattice spacing squared;the corresponding lattice spacings can be found in Table 1. Since the ratios have beencomputed on a line of constant physics, and assuming that we are in a scaling region whereSymazik’s effective theory of cut-off effects applies, they are expected to be polynomialsin the lattice spacing, tending to 1 in the continuum limit. In this context we introducean additional determination, r { u , nu;impr } m , which only differs from r { u , nu } m by an improvedversion of the derivative e ∂ in eq. (5.1). These ratios are very close to one except for oneof the data points at β = 3 .
3, for which the ratio is significantly larger. Even though itwould be sufficient to demonstrate that these ratios of r m approach unity with a rate ∝ a or higher in our particular line of constant physics framework, such ambiguities appear tobe nearly absent for a < . β = 3 . χ / d . o . f . >
3. We thus conclude that ourresults are compatible with the theoretical expectation of O( a ) lattice artefacts or higher(see also Appendix C).As our preferred determination of r m we advocate r { nu; Z,T/ } m because of its smallstatistical errors in our range of bare couplings and the poorer scaling behaviour of theother estimators at the coarsest lattice spacing. In Fig. 7 we show this result including the The improved derivative is defined as a∂ µ f ( x ) ≡ [ − f ( x +2 a ˆ µ )+8 f ( x + a ˆ µ ) − f ( x − a ˆ µ )+ f ( x − a ˆ µ )]and its corresponding second derivative by a ∂ ∗ µ ∂ µ f ( x ) ≡ [ − f ( x + 2 a ˆ µ ) + 16 f ( x + a ˆ µ ) − f ( x ) + 16 f ( x − a ˆ µ ) − f ( x − a ˆ µ )] as shown by eq. (B.4) in Ref. [14]. .
55 1 .
60 1 .
65 1 .
70 1 .
75 1 . g β → χ / d . o . f . = 1.372-loop perturbation theory r { nu;Z , T / } m r m , Bali et al. (2016) Figure 7:
Non-perturbative determination of r { nu; Z,T/ } m (open circles), compared to the resultsof Ref. [13] (filled diamonds) and those of two-loop perturbation theory [18] (horizontal dottedline). The dashed line is the interpolation (5.5) and the vertical dotted lines correspond to the barecouplings used in CLS simulations. two-loop perturbative prediction of Ref. [18]. An important observation is that the non-perturbative estimates strongly deviate from the perturbative prediction in this region ofstrong couplings. A similar behaviour was also observed in several studies of renormalisationfactors for which one-loop perturbative predictions are available (see, e.g. [25, 48]). Here,we confirm this finding also for two-loop perturbation theory. We also compare our resultswith those of other works. In Ref. [13], r m was determined for two values of the barecoupling, from an alternative renormalization condition. As inferred by Fig. 7 this resultagrees with ours at the smaller coupling, while it deviates notably at the larger coupling,most likely due to O( a ) ambiguities (or higher).Our final result consists of a continuous interpolation formula for r m = r { nu; Z,T/ } m .Our data is best described by a Padé ansatz, constrained to the two-loop prediction ofRef. [18] for small couplings, of the form r m ( g ) = 1 . . g × ( c g + c g c g ) , (5.5a)19here c i = ( − . , . , − . , (5.5b)andcov( c i , c j ) = . − . − . × − − . . . × − − . × − . × − . × − , (5.5c)which is also displayed in Fig. 7. The fit function describes our data with χ / d . o . f . = 1 . r m at the couplingsused in CLS simulations for the computation of hadronic quantities [13, 21, 23]. Sincethe
CLS coupling β = 3 .
85 lies outside the range of our r m computations, we performa short extrapolation in order to provide a value for r m ( β = 3 . β = 3 .
81 andthe extrapolated value at β = 3 .
85, is added to the statistical error in quadrature. Ourfinal r m results at the CLS couplings are collected in Table 5. β r m Table 5:
Values for r m at the couplings used in CLS simulations, obtained from the interpolationformula (5.5). As mentioned in the text, an additional systematic error was added to the β = 3 . With the non-perturbative computation of the ratio of the renormalisation constants ofnon-singlet and singlet scalar densities, r m ≡ Z S /Z , presented in this paper we haveaddressed a quantity, which not only enters the renormalisation pattern of quark massesin lattice QCD with Wilson fermions, but also constitutes an important ingredient incalculations of renormalised nucleon (and other baryon) matrix elements of singlet scalardensities, known as sigma terms.Our strategy to calculate r m merges the functional dependences of the PCAC quarkmass in terms of the subtracted quark mass, evaluated in a unitary as well as a non-unitarysetting with respect to the choice of sea and valence quark masses. In the vicinity ofthe chiral limit, these dependences are found to be linear, so that r m can be obtainedthrough the associated quark mass slopes with confidence and superior control of statisticaland systematic errors. The finite-volume numerical simulations of O( a ) improved QCDwith Schrödinger functional boundary conditions that enter the analysis realise a line ofconstant physics by working in a volume of spatial extent L ≈ . r m becomes a smoothfunction of the bare gauge coupling as the lattice spacing is varied, where any potentially20emaining intrinsic ambiguities disappear monotonically towards the continuum limit at arate that stays beyond the sensitivity of the O( a ) improved theory.Our central results, which hold for a lattice discretisation of QCD with three flavoursof non-perturbatively O( a ) improved Wilson-clover sea quarks and tree-level Symanzik-improved gluons, are the continuous parameterisation of r m as a function of the squaredbare gauge coupling g = 6 /β in eq. (5.5), as well as its values in Table 5 at the specificstrong-coupling β values of large-volume CLS simulations [13, 21–23].Along with the numerical implementation of our strategy to extract r m , we have alsodeveloped a new method to determine the scale independent combination Z = Z P / ( Z S Z A )of renormalisation parameters of quark bilinears in the pseudoscalar, (non-singlet) scalarand axial vector channel, respectively. It relies upon a Ward identity that, according toour knowledge, has not yet appeared explicitly in the literature. Since, as explained inSections 2 and 3, the renormalisation factor Z is actually required to isolate r m from theunitary and non-unitary quark mass slopes, we have employed the estimates on Z from thisapproach in our final results of r m . However, this was primarily done for practical reasonsand served the purpose of demonstrating the feasibility of the Ward identity method for Z .In fact, it is apparent from the discussion in Appendix C and Figure 8 that these new valuesfor Z are fully compatible with the earlier determinations available from Refs. [14, 15] andare neither superior in statistical precision nor in systematics regarding lattice artefacts.Thus we refrain from quoting yet another interpolating formula on Z as a function of g . The interested reader may straightforwardly obtain it from the results of Table 6 inAppendix C. Acknowledgements.
This work is supported by the Deutsche Forschungsgemeinschaft (DFG)through the Research Training Group “GRK 2149: Strong and Weak Interactions – from Hadronsto Dark Matter” (J. H., F. J. and P. L. J. P.). We acknowledge the computer resources provided bythe WWU IT, formerly ‘Zentrum für Informationsverarbeitung (ZIV)’, of the University of Münster(PALMA-II HPC cluster) and thank its staff for support. Schrödinger functional correlation functions
The Schrödinger functional correlation functions employed in this work are defined as f ij P = − a L X x , u , v D ¯ ψ i ( x ) γ ψ j ( x ) · ¯ ζ j ( v ) γ ζ i ( u ) E , (A.1) g ij P = − a L X x , u , v D ¯ ψ i ( x ) γ ψ j ( x ) · ¯ ζ j ( u ) γ ζ i ( v ) E , (A.2) f ij A = − a L X x , u , v D ¯ ψ i ( x ) γ γ ψ j ( x ) · ¯ ζ j ( u ) γ ζ i ( v ) E , (A.3) g ij A = − a L X x , u , v D ¯ ψ i ( x ) γ γ ψ j ( x ) · ¯ ζ j ( u ) γ ζ i ( v ) E , (A.4) F ij = − a L X u , v , u , v D ¯ ζ i ( u ) γ ζ j ( v ) · ¯ ζ j ( u ) γ ζ i ( v ) E . (A.5)They refer to the general case of two distinct, i.e. not necessarily mass-degenerate quarkflavours i, j . Summation over the indices i and j is not implied. The space-time point x lies in the lattice bulk; i.e. 0 < x < T . The Dirichlet boundary fields ¯ ζ j ( u ) and ζ i ( v ) liveon time slice x = 0, while ¯ ζ j ( u ) and ζ i ( v ) live on time slice x = T ; the boundary fieldsare introduced in Ref. [36]. B Wick contractions of correlation functions
In this appendix we briefly explain how to obtain eq. (3.5) from eq. (3.4). The idea is toperform the Wick contractions of the correlation functions, arriving at expressions whichare traces of flavour matrices, multiplying traces of products of quark propagators and γ -matrices. This procedure has been described in full detail in Ref. [15], which dealswith more complicated Ward identities; we refer the reader to that work for unexplainednotation. Here we will only present the main features of the proof.We start with the r.h.s. of eq. (3.4). The Wick contractions result in − d abe a X y h P e ( y ) O c i = − d abe Tr[ T e T c ] a L X y X u , v * tr (cid:26) [ ψ ( y )¯ ζ ( u )] F γ [ ζ ( v ) ¯ ψ ( y )] F γ (cid:27)+ = d abc f P ( y ) , (B.1)where the second equality implicitly defines f P (see also eq. (A.1) and Appendix B ofRef. [14]). The left-hand-side consists of correlation functions with one boundary operatorand two insertions in the bulk. So the Wick contractions of such a correlation functiongive: 22 X x , y h A a ( x ) S b ( y ) O c i = i a L Tr[ T a T b T c ] X x , y X u , v * tr (cid:26) γ γ [ ψ ( x ) ¯ ψ ( y )] F [ ψ ( y )¯ ζ ( u )] F γ [ ζ ( v ) ¯ ψ ( x )] F (cid:27)+ + i a L Tr[ T c T b T a ] X x , y X u , v * tr (cid:26) γ γ [ ψ ( x )¯ ζ ( u )] F γ [ ζ ( v ) ¯ ψ ( y )] F [ ψ ( y ) ¯ ψ ( x )] F (cid:27)+ = iTr[ T a T b T c ] f AS;1 ( x , y ) + iTr[ T c T b T a ] f AS;2 ( x , y )= 12 (cid:20) − d abc Re f AS;1 ( x , y ) + f abc Im f AS;1 ( x , y ) (cid:21) . (B.2) The second in the above string of equations implicitly defines the two traces of quarkpropagators (devoid of flavour structure) as f AS;1 ( x , y ) and f AS;2 ( x , y ). In the lastequation we have made use of the fact that the two traces of propagators are complexconjugates of each other which is a consequence of the γ -Hermiticity property of Wilsonfermion propagators. Finally, the fact that the above correlation function is invariantunder charge conjugation leads to the vanishing of the term proportional to f abc in thelast expression. Hence, we obtain a X x , y h A a ( x ) S b ( y ) O c i = − d abc Re f AS;1 ( x , y ) = − d abc f AS ( x , y ) , (B.3)which implicitly defines f AS ( x , y ). The correlation functions f P and f AS are schematicallydrawn in Fig. 1. Analogously, from the mass dependent term of the Ward identity we alsodefine ˜ f PS ( x , y ); the summation over all times from t up to t (see eq. (3.4)) is includedin its definition. It is important to note that d abc appears in both eqs. (B.1) and (B.3).Therefore, it cancels out in the Ward identity, which becomes an expression between tracesof propagators, without any flavour indices. Putting everything together, we eventuallyobtain eq. (3.5). C Comparison of Z determinations and scaling tests In this appendix we present more details on our Z results, listed in Table 2. In Fig. 8 andTable 6 our preferred determination for Z , namely Z { T/ } is compared to Ref. [14] (deDivitiis et al.), Z determined at two values of the gauge coupling in Ref. [32] (Bali et al.)and to a Z estimate that we work out from the results of Refs. [48] and [15] (Heitger etal.). In particular, we extract the axial current normalisation Z A at our couplings from theinterpolation formula of Ref. [48] and combine it with the ratio Z S /Z P of the pseudoscalarand scalar renormalisation constants from Ref. [15].Our result agrees with the other determinations at weaker bare couplings, whiledisagreements are seen at stronger couplings. These are attributed to lattice artefactsassociated with intrinsic ambiguities of O( a ) or higher between different determinations.Agreement is generally better between our results and those of Ref. [14] (de Divitiis et al.).In order to confirm this claim of consistency (leaving aside higher cut-off effects) weconstruct ratios of different determinations and investigate their behaviour as a function ofthe lattice spacing. Interestingly, rather than O( a ), leading cut-off effects of O( a ) canbe identified in the ratio Z { T/ } /Z { T/ } , as seen in Fig. 9 (left). The scaling behaviour of23 .
55 1 .
60 1 .
65 1 .
70 1 .
75 1 . g . . . . . . . . Z = Z { T/ } , this work Z , LCP-0, de Divitiis et al. (2019)1 /Z A · Z P /Z S , Heitger et al. (2020) Z , Bali et al. (2016) Figure 8: Z results, obtained with different methods, as a function of the squared bare coupling g . The preferred determination of this work is Z = Z { T/ } (squares). The pentagons are obtainedby combining results from Refs. [48] and [24] (Heitger et al.). The two Z estimates determinedin Ref. [32] (Bali et al.) are depicted by triangles. The circles correspond to the Z results fromRef. [14] (de Divitiis et al.). One-loop perturbation theory is illustrated by the dotted line, Ref. [14]. β Z = Z { T/ } this work Z , LCP-0de Divitiis et al. Z , LCP-1de Divitiis et al. 1 /Z A · Z P /Z S Heitger et al.
Table 6:
Comparison of our preferred Z determination with results from Ref. [14] (de Divitiis etal.) and the combination of results from Refs. [48] and [24] (Heitger et al.). .
000 0 .
002 0 .
004 0 .
006 0 .
008 0 . a in fm . . . . . f ( a ) = 1 + ca χ / d . o . f . = 0.87 Z { T/ } /Z { T/ } .
000 0 .
002 0 .
004 0 .
006 0 .
008 0 .
010 0 . a in fm . . . . . f ( a ) = 1 + ca , χ / d . o . f . = 0.35 f ( a ) = 1 + ca , χ / d . o . f . = 0.37 Z { T/ } /Z , LCP-0,de Divitiis et al. (2019)(1 /Z A · Z P /Z S ) /Z { T/ } ,Heitger et al. (2020) Figure 9:
Left:
Ratio of our results Z { T/ } /Z { T/ } , fitted as a function of the lattice spacing. Right:
Ratios of Z { T/ } and Z computed in previous works, fitted as a function of the latticespacing. The coarsest lattice spacing is excluded from the fits. The dashed lines are the fits whilethe horizontal dotted lines indicate the expected continuum results. our results compared to those of previous works is shown in Fig. 9 (right). All ratios arefitted with an ansatz 1 + ca , excluding the coarsest lattice spacing. When adding a termlinear in the lattice spacing, its fit parameter vanishes within its uncertainty in all cases.In conclusion, these scaling tests indicate that our results for Z are in accordance with thetheoretical expectation of O( a ) ambiguities or higher which by virtue of the imposed lineof constant physics decrease monotically towards the continuum limit. References [1] M. Bochicchio, L. Maiani, G. Martinelli, G. C. Rossi and M. Testa,
Chiral Symmetry on theLattice with Wilson Fermions , Nucl. Phys. B (1985) 331.[2] L. Maiani, G. Martinelli, M. L. Paciello and B. Taglienti,
Scalar Densities and Baryon MassDifferences in Lattice QCD With Wilson Fermions , Nucl. Phys. B (1987) 420.[3] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas,
A General method fornonperturbative renormalization of lattice operators , Nucl. Phys. B (1995) 81,[ hep-lat/9411010 ].[4] C. Sturm, Y. Aoki, N. H. Christ, T. Izubuchi, C. T. C. Sachrajda and A. Soni,
Renormalization of quark bilinear operators in a momentum-subtraction scheme with anonexceptional subtraction point , Phys. Rev. D (2009) 014501, [ ].[5] S. Capitani, M. Lüscher, R. Sommer and H. Wittig, Non-perturbative quark massrenormalization in quenched lattice QCD , Nucl. Phys. B (1999) 669, [ hep-lat/9810063 ].[Erratum: Nucl. Phys. B 582, 762 (2000)].[6] M. Dalla Brida, S. Sint and P. Vilaseca,
The chirally rotated Schrödinger functional:theoretical expectations and perturbative tests , JHEP (2016) 102, [ ].[7] G. M. de Divitiis and R. Petronzio, Nonperturbative renormalization constants on the latticefrom flavor nonsinglet Ward identities , Phys. Lett. B (1998) 311, [ hep-lat/9710071 ].
8] T. Bhattacharya, S. Chandrasekharan, R. Gupta, W.-J. Lee and S. R. Sharpe,
Nonperturbative renormalization constants using Ward identities , Phys. Lett. B (1999) 79,[ hep-lat/9904011 ].[9] T. Bhattacharya, R. Gupta, W.-J. Lee and S. R. Sharpe,
Order a improved renormalizationconstants , Phys. Rev. D (2001) 074505, [ hep-lat/0009038 ].[10] M. Guagnelli, R. Petronzio, J. Rolf, S. Sint, R. Sommer and U. Wolff, Non-perturbative resultsfor the coefficients b m and b a − b P in O( a ) improved lattice QCD , Nucl. Phys. B (2001)44, [ hep-lat/0009021 ].[11] T. Bhattacharya, R. Gupta, W. Lee and S. R. Sharpe,
Scaling behavior of discretization errorsin renormalization and improvement constants , Phys. Rev. D (2006) 114507,[ hep-lat/0509160 ].[12] P. Fritzsch, J. Heitger and N. Tantalo, Non-perturbative improvement of quark massrenormalization in two-flavour lattice QCD , JHEP (2010) 074, [ ].[13] G. S. Bali, E. E. Scholz, J. Simeth and W. Söldner, Lattice simulations with N f = 2 + 1 improved Wilson fermions at a fixed strange quark mass , Phys. Rev. D (2016) 074501,[ ].[14] G. M. de Divitiis, P. Fritzsch, J. Heitger, C. C. Köster, S. Kuberski and A. Vladikas, Non-perturbative determination of improvement coefficients b m and b A − b P and normalisationfactor Z m Z P /Z A with N f = 3 Wilson fermions , Eur. Phys. J. C (2019) 797, [ ].[15] J. Heitger, F. Joswig and A. Vladikas, Ward identity determination of Z S /Z P for N f = 3 lattice QCD in a Schrödinger functional setup , Eur. Phys. J. C (2020) 765, [ ].[16] G. Bali, S. Bürger, S. Collins, M. Göckeler, M. Gruber, S. Piemonte et al., NonperturbativeRenormalization in Lattice QCD with three Flavors of Clover Fermions: Using Periodic andOpen Boundary Conditions , .[17] T. Bhattacharya, R. Gupta, W. Lee, S. R. Sharpe and J. M. Wu, Improved bilinears in latticeQCD with non-degenerate quarks , Phys. Rev. D (2006) 034504, [ hep-lat/0511014 ].[18] M. Constantinou, M. Hadjiantonis, H. Panagopoulos and G. Spanoudes, Singlet versusnonsinglet perturbative renormalization of fermion bilinears , Phys. Rev. D (2016) 114513,[ ].[19] M. Lüscher and P. Weisz, On-Shell Improved Lattice Gauge Theories , Commun. Math. Phys. (1985) 59. [Erratum: Commun. Math. Phys. 98, 433 (1985)].[20] B. Sheikholeslami and R. Wohlert, Improved Continuum Limit Lattice Action for QCD withWilson Fermions , Nucl. Phys. B (1985) 572.[21] M. Bruno et al.,
Simulation of QCD with N f = 2 + 1 flavors of non-perturbatively improvedWilson fermions , JHEP (2015) 043, [ ].[22] M. Bruno, T. Korzec and S. Schaefer, Setting the scale for the CLS flavor ensembles , Phys. Rev. D (2017) 074504, [ ].[23] D. Mohler, S. Schaefer and J. Simeth, CLS flavor simulations at physical light- andstrange-quark masses , EPJ Web Conf. (2018) 02010, [ ].[24] J. Heitger, F. Joswig, A. Vladikas and C. Wittemeier,
Non-perturbative determination of c V , Z V and Z S /Z P in N f = 3 lattice QCD , EPJ Web Conf. (2018) 10004, [ ].[25] J. Heitger and F. Joswig,
The renormalised O( a ) improved vector current in three-flavourlattice QCD with Wilson quarks , .
26] J. Bulava, M. Della Morte, J. Heitger and C. Wittemeier,
Nonperturbative renormalization ofthe axial current in N f = 3 lattice QCD with Wilson fermions and a tree-level improved gaugeaction , Phys. Rev. D (2016) 114513, [ ].[27] J. Bulava, M. Della Morte, J. Heitger and C. Wittemeier, Non-perturbative improvement ofthe axial current in N f = 3 lattice QCD with Wilson fermions and tree-level improved gaugeaction , Nucl. Phys. B (2015) 555, [ ].[28] L. Chimirri, P. Fritzsch, J. Heitger, F. Joswig, M. Panero, C. Pena et al.,
Non-perturbativerenormalization of O ( a ) improved tensor currents , PoS
LATTICE2019 (2020) 212,[ ].[29] M. Bruno, I. Campos, P. Fritzsch, J. Koponen, C. Pena, D. Preti et al.,
Light quark masses in N f = 2 + 1 lattice QCD with Wilson fermions , Eur. Phys. J. C (2020) 169, [ ].[30] J. Heitger, F. Joswig and S. Kuberski, Determination of the charm quark mass in lattice QCDwith flavours on fine lattices , .[31] G. S. Bali et al., The strange and light quark contributions to the nucleon mass from LatticeQCD , Phys. Rev. D (2012) 054502, [ ].[32] G. S. Bali, S. Collins, D. Richtmann, A. Schäfer, W. Söldner and A. Sternbeck, Directdeterminations of the nucleon and pion σ terms at nearly physical quark masses , Phys. Rev. D (2016) 094504, [ ].[33] S. Aoki et al., FLAG Review 2019: Flavour Lattice Averaging Group (FLAG) , Eur. Phys. J. C (2020) 113, [ ].[34] K. Ottnad, Excited states in nucleon structure calculations , in , 11, 2020. .[35] J. Green,
Systematics in nucleon matrix element calculations , PoS
LATTICE2018 (2018)016, [ ].[36] M. Lüscher, S. Sint, R. Sommer and P. Weisz,
Chiral symmetry and O( a ) improvement inlattice QCD , Nucl. Phys. B (1996) 365, [ hep-lat/9605038 ].[37] J. Bulava and S. Schaefer,
Improvement of N f = 3 lattice QCD with Wilson fermions andtree-level improved gauge action , Nucl. Phys. B (2013) 188, [ ].[38] M. Lüscher and S. Schaefer. http://luscher.web.cern.ch/luscher/openQCD .[39] A. D. Kennedy, I. Horvath and S. Sint,
A New exact method for dynamical fermioncomputations with nonlocal actions , Nucl. Phys. Proc. Suppl. (1999) 834,[ hep-lat/9809092 ].[40] M. A. Clark and A. D. Kennedy, Accelerating dynamical fermion computations using therational hybrid Monte Carlo (RHMC) algorithm with multiple pseudofermion fields , Phys. Rev.Lett. (2007) 051601, [ hep-lat/0608015 ].[41] M. Lüscher, S. Sint, R. Sommer, P. Weisz and U. Wolff, Non-perturbative O(a) improvementof lattice QCD , Nucl. Phys. B (1997) 323, [ hep-lat/9609035 ].[42] L. Del Debbio, H. Panagopoulos and E. Vicari, θ dependence of SU(N) gauge theories , JHEP (2002) 044, [ hep-th/0204125 ].[43] P. Fritzsch, A. Ramos and F. Stollenwerk, Critical slowing down and the gradient flowcoupling in the Schrödinger functional , PoS
Lattice2013 (2014) 461, [ ].[44] U. Wolff,
Monte Carlo errors with less errors , Comput. Phys. Commun. (2004) 143,[ hep-lat/0306017 ]. [Erratum: Comput. Phys. Commun. 176, 383 (2007)].
45] S. Schaefer, R. Sommer and F. Virotta,
Critical slowing down and error analysis in latticeQCD simulations , Nucl. Phys. B (2011) 93, [ ].[46] A. Ramos,
Automatic differentiation for error analysis of Monte Carlo data , Comput. Phys.Commun. (2019) 19, [ ].[47] P. T. Boggs and J. E. Rogers,
Orthogonal distance regression , tech. rep., National Institute ofStandards and Technology, Gaithersburg, MD, 1989. 10.6028/NIST.IR.89-4197.[48] M. Dalla Brida, T. Korzec, S. Sint and P. Vilaseca,
High precision renormalization of theflavour non-singlet Noether currents in lattice QCD with Wilson quarks , Eur. Phys. J. C (2019) 23, [ ].].