Direct determination of the strange and light quark condensates from full lattice QCD

# Direct determination of the strange and light quark condensates from full lattice QCD

C. McNeile Bergische Universität Wuppertal, Gaussstr. 20, D-42119 Wuppertal, Germany    A. Bazavov Physics Department, Brookhaven National Laboratory, Upton NY 11973, USA    C. T. H. Davies SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK    R. J. Dowdall SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow, G12 8QQ, UK    K. Hornbostel Southern Methodist University, Dallas, Texas 75275, USA    G. P. Lepage Laboratory of Elementary-Particle Physics, Cornell University, Ithaca, New York 14853, USA    H. D. Trottier Physics Department, Simon Fraser University, 8888 University Drive, Burnaby, BC, V5A 1S6, Canada TRIUMF, 4004 Westbrook Mall, Vancouver, BC, V6T 2A3, Canada
July 21, 2019
###### Abstract

We determine the strange quark condensate from lattice QCD for the first time and compare its value to that of the light quark and chiral condensates. The results come from a direct calculation of the expectation value of the trace of the quark propagator followed by subtraction of the appropriate perturbative contribution, derived here, to convert the non-normal-ordered to the scheme at a fixed scale. This is then a well-defined physical ‘nonperturbative’ condensate that can be used in the Operator Product Expansion of current-current correlators. The perturbative subtraction is calculated through and estimates of higher order terms are included through fitting results at multiple lattice spacing values. The gluon field configurations used are ‘second generation’ ensembles from the MILC collaboration that include flavors of sea quarks implemented with the Highly Improved Staggered Quark action and including sea quarks down to physical masses. Our results are : , , where is a light quark with mass equal to the average of the and quarks. The strange to light quark condensate ratio is 1.08(16). The light quark condensate is significantly larger than the chiral condensate in line with expectations from chiral analyses. We discuss the implications of these results for other calculations.

## I Introduction

A critical feature of the nonperturbative dynamics of QCD at zero temperature is the condensation of quark-antiquark pairs in the vacuum, spontaneously breaking the chiral symmetry of the action. The value of the chiral condensate (the quark condensate at zero quark mass) is then an important parameter for low energy QCD Reinders et al. (1985). The well-known Gell-Mann, Oakes, Renner (GMOR) relation Gell-Mann et al. (1968):

 f2πM2π4=−mu+md2⟨0|¯¯¯uu+¯¯¯dd|0⟩2 (1)

connects the quark masses times condensate to the square of the mass times decay constant for the Goldstone boson of the spontaneously broken symmetry. Eq.( 1) has normalisation such that MeV. The GMOR relation holds in the limit of . A value for this chiral condensate can be derived from the chiral extrapolation of lattice QCD results for light meson masses and decay constants. See, for example, the recent result of for the chiral condensate in the scheme at 2 GeV using SU(2) chiral perturbation theory Borsanyi et al. (2012).

The determination of the quark condensate for non-zero quark masses is more problematic because, depending on the method used, there are various sources of unphysical quark mass dependence and a careful definition of the condensate is required. This definition must be phrased in terms of the Operator Product Expansion (OPE) since this is the context in which the condensate appears Shifman et al. (1979); Reinders et al. (1985); Shifman (1998). The OPE allows us to separate short and long-distance contributions in, for example, a short-distance current-current correlator. The expansion is in terms of a set of matrix elements of local operators multiplied by coefficient functions. The aim is for all the long-distance contributions (with scale ) to be contained in the matrix elements and the short distance contributions (with scale ) in the coefficient functions. A key matrix element, since it corresponds to a relatively low-dimensional () operator, is that of the quark condensate. The clean separation of scales in the OPE only works if the local operators are not normal ordered Chetyrkin et al. (1985); Jamin and Munz (1993). Then the coefficient functions are analytic in the quark masses and therefore free of infrared sensitivity. This means, however, that the quark mass dependent mixing of with the unit operator must be taken into account and that the vacuum matrix element of is not cut-off independent. The quantity that appears in the OPE is the vacuum matrix element in, for example, the scheme at the scale . We can derive this matrix element from lattice QCD and we give results here for = 2 GeV. The results can easily be run to other scales, as appropriate.

The value of the condensate for quarks of non-zero mass up to that of the strange quark is needed in a number of calculations involving light quark correlators. In lattice QCD it is frequently easier and statistically more precise to use strange quarks than very light quarks in contexts where the quark mass is not expected to be important. Then the strange quark condensate appears in the calculation, however. Examples include the matching to continuum QCD perturbation theory of lattice QCD calculations of moments of heavy-light meson correlators Koponen et al. (2010) and of light meson correlators at large space-like  Shintani et al. (2010). Such calculations are used to extract quark masses and the strong coupling constant, . A continuum example where the strange quark condensate is needed is in the determination of the strange quark mass, , from hadronic decays Gamiz et al. (2003).

Current estimates of the value of the strange quark condensate vary by almost a factor of two Albuquerque et al. (2010); Maltman (2009). It is not even clear whether the strange condensate is larger or smaller than the light quark condensate. For very large quark masses, , say, so that the quark mass dominates the propagator, it seems clear that the condensate should fall to zero, but this does not help in determining the slope of the condensate with for small quark masses.

Here we address the determination of the strange (or other non-zero mass) quark condensate by direct calculation in full lattice QCD. By direct we mean that we determine the vacuum expectation value of the strange quark propagator as well as the light quark propagator on a range of gluon field configurations at different values of the lattice spacing and sea quark masses. To isolate the low-energy nonperturbative value of the condensate from these results requires the subtraction of a perturbative contribution. The perturbative contribution in lattice QCD has two pieces. One diverges as and dominates the vacuum expectation value of the strange quark propagator, particularly on our finer lattices. The second piece contains infrared sensitive logarithms of the quark mass which cancel against similar terms in continuum perturbation theory allowing an infrared safe definition of the condensate for use in the OPE, as discussed above.

The error in the final result then depends on how well this subtraction can be done. Here we use an explicit calculation of the perturbative pieces through and fit for unknown higher order terms. The known quark mass and dependence of these unknown terms helps in constraining them along with the very small statistical errors in our lattice results. We also use a particularly good discretisation of the Dirac action known as the Highly Improved Staggered Quark (HISQ) formalism Follana et al. (2007) on ‘second generation’ gluon field configurations so that discretisation errors in the physical nonperturbative results are small.

The paper is laid out as follows. In Section II we describe the theoretical background to direct calculations of the quark condensate in lattice QCD. Section III gives our lattice QCD results on gluon configurations with 2+1+1 flavors of sea quarks, describing the calculation of the perturbative contribution that is subtracted and then the procedure for fitting the remaining nonperturbative condensates as a function of quark mass and lattice spacing. We also give results from configurations including 2+1 flavors of sea quarks over a wider range of lattice spacing values but studying only the strange quark condensate in Appendix D. In Section IV we compare to previous values and discuss the implications of our results for both zero and finite temperature QCD calculations. Section V gives our conclusions.

## Ii Theoretical Background

The direct determination of the chiral condensate in lattice QCD requires the calculation of the expectation value over an ensemble of gluon fields, , of where is the lattice discretisation of the Dirac matrix. The quark action for a given quark flavor,

 Sf=¯¯¯¯ψMfψ (2)

and

 ⟨¯¯¯¯ψψ⟩=⟨0|¯¯¯¯ψfψf|0⟩=−1V⟨TrMf(U)−1⟩U, (3)

where the trace is over spin, color and space-time point and the gluon fields in the ensemble used for the average include the effect of sea quarks (of all flavors, not just ) in their probability distribution. is the lattice volume, . For a naive discretisation of the Dirac action takes the form:

 M=γμΔμ+m (4)

where is a covariant finite difference on the lattice:

 Δμψx=12a(Uμ(x)ψ(x+^μ)−U†μ(x−^μ)ψ(x−^μ)) (5)

and is the quark mass for that flavor. Because of fermion doubling this formalism describes 16 ‘tastes’ of quarks in 4-dimensions rather than just 1 and we must divide the right-hand side of Eq. (3) by the number of tastes, . The staggered formalism is derived from this naive formalism by a rotation which allows the spin degree of freedom to be dropped. In that case the quark field becomes a 1-component spinor, which is numerically very efficient, and .

For the eigenvalues of for either naive or staggered quarks, are purely imaginary and come in pairs. Therefore, in the absence of exact zero modes,

 −⟨¯¯¯¯ψψ⟩ = 1Nt∑λ(1m+iλ+1m−iλ) (6) = 1Nt∑λ2mm2+λ2.

A calculation at on a finite volume lattice would then give an answer for the quark condensate of zero. This does not mean that chiral symmetry is unbroken, however. The problem arises because the broad distribution of non-zero eigenvalues (ignoring topological near-zero modes) drops to zero near the origin in a way that depends strongly on the volume. Once the quark mass is below the minimum of the non-zero eigenvalues the result for will be distorted. A more careful consideration of limits must be made. If is taken to infinity before is taken to zero then the sum over eigenvalues can be replaced by an integral and the Banks-Casher relation Banks and Casher (1980) is obtained. This connects the zero-mass condensate to the spectral density at the origin:

 Σ=−⟨¯¯¯¯ψψ(m=0)⟩=πρ(0)Nt. (7)

Thus can be obtained from studies of the spectral density and, more recently, has also been obtained from matching the distribution of low eigenmodes to random matrix theory in the regime Leutwyler and Smilga (1992); Shuryak and Verbaarschot (1993); Damgaard (2002). Here we are more concerned with extracting a condensate at non-zero quark mass, for example at the strange quark mass, and so the issue above is not relevant. We will work on large volume lattices (over ten times larger than the study in Follana et al. (2005) that looked at staggered eigenvalues in the regime) at quark mass values that are well within the distribution of non-zero eigenvalues. Our results for then include both the effects of a non-zero value for and a non-zero quark mass.

A second issue arises, however, in the extraction of a physical nonperturbative value for the condensate at non-zero quark mass. A perturbative contribution appears from mixing between the scalar operator and the identity since the identity operator has a vacuum expectation value in the trivial perturbative vacuum. This perturbative contribution vanishes at zero quark mass since chiral symmetry is not broken in perturbation theory (for the same reason as that given on the lattice above in Eq. (6)). At non-zero quark mass it contains odd powers of starting with a quadratically ultra-violet divergent term linear in . This can be illustrated with a tree-level calculation in the continuum to give:

 −⟨¯¯¯¯ψψ⟩ = ∫Λ0d4k(2π)412mk2+m2 = 34π2(mΛ2+m3logm2Λ2+m2).

The quadratic ultraviolet divergence depends on the scheme used but the term is universal since it arises from the infrared part of the integral. The coefficient above agrees with that obtained for the scheme in Braaten et al. (1992) and on the lattice for highly improved staggered quarks to be described in section III.1. We stress that a perturbative contribution of this kind is present for all lattice regularisations of QCD, whatever their chiral symmetry properties, and so must be calculated and subtracted to give a physical result. The quadratic divergence present for Ginsparg-Wilson fermions is demonstrated in the quenched approximation in, for example, Hernandez et al. (1999) and the additional divergences for the Wilson formalism with broken chiral symmetry in David and Hamber (1984).

This subtraction is somewhat analogous to subtracting perturbative contributions to the mean plaquette to obtain the nonperturbative gluon condensate. That, however, is extremely difficult to do because the nonperturbative condensate contribution to the plaquette is so small. This contribution is given at leading order by:

 δPcond=−π236a4⟨αsG2/π⟩. (9)

If we take the value of the gluon condensate as then . On very coarse lattices for which , this contributes less than 1% to the value of the plaquette. On finer lattices the nonperturbative condensate contribution is even smaller because it falls as while the perturbative contribution falls only as for some scale . This means that the plaquette is in fact a very good variable to use for the determination of from lattice QCD calculations but not for the determination of the gluon condensate Davies et al. (2008). For larger Wilson loops the gluon condensate contribution is larger, being proportional to the square of the area of the loop, but the coefficients in the perturbative series also become larger.

The determination of the nonperturbative quark condensate from is in much better shape than this for several reasons. The main one is that the nonperturbative condensate contribution to in lattice units is only a factor of smaller than the leading perturbative contribution rather than . In addition the perturbative contribution is suppressed by the quark mass, which is small for the and quarks we will consider here. The perturbative contribution is a well-defined function of the quark mass at every order in perturbation theory and so results at several values of the quark mass, and the lattice spacing, can be used to constrain unknown higher orders, beyond the that we have explicitly calculated, and we will make use of that here.

## Iii Lattice QCD calculation on nf=2+1+1 gluon configurations

The gluon field configurations used here are listed in Table 1. They were generated by the MILC collaboration Bazavov et al. (2010a) using a tadpole-improved Lüscher-Weisz gauge action with coefficients corrected perturbatively through including pieces proportional to , the number of quark flavors in the sea Hart et al. (2009). The gauge action is then improved completely through . Sea quarks are included with the highly improved staggered quark (HISQ) action Follana et al. (2007) which has been designed to have very small discretisation errors. Discretisation errors are formally removed through but higher order errors, particularly staggered taste-changing errors, are seen to be smaller with HISQ than with the earlier asqtad staggered quark action Follana et al. (2007); Bazavov et al. (2010a). The HISQ action used here has two smearing steps for the gluon field appearing in the quark action with a U(3) projection of the smeared links between the two smearing steps. The configurations include a sea charm quark in addition to up, down and strange. These configurations are then said to have 2+1+1 flavors in the sea, since the and quarks are taken to have the same mass (denoted here). This is heavier than the average mass in the real world on most of the configuration sets but there are three for which the mass has its physical value (3, 6 and 8). The and masses are tuned as closely as possible to their correct values on each set. The tuning of the sea quark mass is accurately done – typically to better than 5% – so the quark mass can be accurately calibrated in terms of the quark mass for chiral extrapolations. The values of the lattice spacing for most ensembles were determined in Dowdall et al. (2011) using the decay constant of the meson. The values vary from 0.15 fm to 0.09 fm as we go from the very coarse to the fine lattices. The spatial volumes are large, from when to when .

On each of these ensembles we determine for HISQ valence quarks for various quark masses. To do this we use an identity that relates the quark propagator for staggered quarks to a product of quark propagators:

 1amqTrM−100 = ∑nTr[M−10nM−1n0](−1)n (10) = ∑nTr|M−10n|2

Here 0 and are arbitrary lattice sites and is the quark mass in lattice units used for the quark propagator. The righthand side of Eq. (10) is simply the correlator between 0 and for the Goldstone pseudoscalar meson made of a quark and antiquark of mass . Summing over projects on to zero spatial momentum and sums over timeslices. Thus, dividing both sides by 4, the number of tastes for staggered quarks, we obtain:

 −a3⟨¯¯¯¯ψψ⟩0=(amq)∑tCπ(t). (11)

The raw condensate value on the left-hand side of this expression is normalised to the single flavor case and the pion correlator on the right-hand side is the usual zero-momentum Goldstone meson correlator. This allows us to determine by summing over the Goldstone pseudoscalar correlators calculated in Dowdall et al. (2011). Eq. (10) is derived in Kilcup and Sharpe (1987) for a single propagator origin, 0, but the derivation can trivially be extended to hold for the random wall source that we use for our correlators in Dowdall et al. (2011). The identity holds configuration by configuration for lattice QCD quark formalisms with sufficient chiral symmetry and in the continuum for a specific gauge field background. We give an explicit proof of this in Appendix A. Since our Goldstone pseudoscalar correlators are sums of positive numbers they are particularly precise and this precision then carries over to our condensate results. For the light condensate we use the pion correlators made of light quarks and for the strange condensate we use the correlators made of strange quarks. We stress that what we calculate using Eq. (10) is the vacuum expectation value of the condensate (and not a specific ‘in-meson’ value) despite the fact that we determine it for convenience from a meson correlator.

Table 2 gives the valence quark masses used in our calculation and the raw results for the condensate obtained from Eq. (10). The correlator calculations used 16 ‘random wall’ time sources on approximately one thousand configurations in each ensemble (somewhat fewer on sets 6 and 8) and so the results have very small statistical errors. The valence light quark masses are equal to those in sea (except for a very small change on set 3) but we have shifted the valence strange quark masses slightly to be closer to the physical strange quark mass, following Dowdall et al. (2011). On sets 1 and 2 we give results for two different values of the strange quark mass, to help in constraining the valence mass dependence of the condensate.

Errors on the condensate values are determined after binning over adjacent sets of at least five configurations, following analysis of the autocorrelation function. An example plot is shown, for coarse set 6, in Fig. 1. The autocorrelation function is defined as:

 CΔT=⟨xixi+ΔT⟩−⟨xi⟩⟨xi+ΔT⟩⟨x2i⟩−⟨xi⟩2. (12)

Here represents a condensate value on configuration and that on a configuration a further time units along in the ordered ensemble. thus corresponds to adjacent configurations. Fig. 1 shows that nearby configurations in the ensembles are correlated and thus binning is necessary to obtain a reliable statistical error. A similar analysis applies to masses and decay constants as discussed in Dowdall et al. (2011).

In the next section we describe the perturbative calculation of the condensate which we will then subtract from the raw results of Table 2 to enable the nonperturbative condensate to be determined.

### iii.1 Perturbative calculation of ⟨TrM−1⟩

We computed the perturbative contribution to the chiral condensate for the HISQ action through first-order in :

 −a3⟨¯¯¯¯ψψ⟩PT,HISQ = am0× [c0(am0)+c1(am0)αs+O(α2s)],

where is the bare quark mass parameter that appears in the HISQ action. The Feynman diagrams required to this order are shown in Fig. 2. The perturbative quadratic ultraviolet divergence discussed in eq.(II) shows up as finite values for the perturbative coefficients, as defined above, in the limit .

We computed the coefficients from numerical evaluation of the lattice loop integrals over a range of masses that includes the light and strange quark masses that we have used. A representative sample of our results is given in Table 3, and is illustrated in Fig. 3.

An excellent fit to the perturbative coefficients in this range of small quark masses, , can be obtained using the following parameterizations

 c0(am0)=c00+(am0)2[c01log(am0)+c02], (14)

and

 c1(am0) =c10+(am0)2[c11log2(am0) +c12log(am0)+c13]. (15)

Higher order terms in appear as discretisation errors in the comparison to to be done below and so can be ignored - they are negligible for the masses we are using in any case. The leading logarithm of at each order originates entirely from the infrared region of the loop momenta, and the respective coefficients and can easily be computed analytically. The values of these coefficients must and do agree with the values in the scheme Braaten et al. (1992). At one-loop we also have a constraint on the sub-leading (single) logarithm of since, as discussed in Appendix C, all terms must vanish in the difference between the vacuum expectation values of in perturbation theory in the continuum and on the lattice. Allowing for the renormalisation between the mass and the HISQ bare mass:

 ¯¯¯¯¯m(μ)=m0(1+αs[−2πlogaμ+0.1143(3)]+…), (16)

we find that should have the value 0.2307(2). With the logarithmic terms fixed to their known values we can obtain the other coefficients in eqs. (14) and (15) from a fit to the values for and as a function of . We find:

 c00 =0.38366(1), c01 =3/(2π2), c02 =−0.153(1), (17)

and

 c10 =0.03657(7), c11 =−6/π3, c12 =0.2307(2), c13 =0.308(15). (18)

These fits are illustrated in Fig. 3 and reproduce our results for the coefficients to within their numerical integration errors, which are smaller than about for , and for .

The perturbative determination of the vacuum expectation value of has also been done in the scheme, in Braaten et al. (1992). The power divergence is missing in this case but, as discussed above, there are terms proportional to .  Braaten et al. (1992) finds:

 − ⟨¯¯¯¯ψψ⟩(μ)PT,¯¯¯¯¯¯¯MS=¯¯¯¯¯m3(μ)× [d01lm+d02+αs(d11l2m+d12lm+d13)+…]

where and

 d01 = c01=32π2 d02 = −34π2 d11 = c11=−6π3 d12 = 5π3 d13 = −52π3 (20)

As discussed in Appendix B we must subtract the difference between the lattice QCD and perturbative calculations from our lattice QCD results to obtain the nonperturbative condensate in the scheme at the scale . We work with the combination which would be RG-invariant in the absence of this perturbative contribution and it is convenient to derive the subtraction needed in lattice units and as a function of the bare lattice quark mass. Using eq. (16) we obtain

 ΔPT = −a4(⟨m0¯¯¯¯ψψ⟩PT,HISQ−⟨¯¯¯¯¯m(μ)¯¯¯¯ψψ⟩PT,¯¯¯¯¯¯¯MS) = c00(am0)2+αsc10(am0)2 + (am0)4[c01lμ−0.077(1)] + αs(am0)4[c11l2μ+0.1340(2)lμ+0.406(15)]+…,

where . This difference of perturbative expansions is now free of all logarithms of and therefore well-defined and infrared safe.

### iii.2 Determining the nonperturbative strange and light quark condensates

#### iii.2.1 A first look at the results

The physical condensate in the scheme at the scale is then defined by:

 ⟨m¯¯¯¯ψψ⟩NP,¯¯¯¯¯¯¯¯MS(μ)=a−4(a4⟨m¯¯¯¯ψψ⟩0−ΔPT), (22)

where is the numerical result from lattice QCD and is given through as a function of the quark mass in Eq. (III.1). will also contain unknown higher order pieces in that we can try to determine from a fit to the lattice QCD results. First we look at the effect of the calculated tree-level and one-loop contributions.

is a strong function of the quark mass in lattice units, dominated by the terms that give rise to the quadrative divergence with inverse lattice spacing. This means that the relative size of the subtraction compared to the raw results varies strongly with quark mass and with lattice spacing, and this is reflected in the raw results before the subtraction is made. In Fig. 4 the open squares show the unsubtracted results (i.e. setting to zero in Eq. (22)) as a function of the square of the inverse lattice spacing for quarks at the four different masses that we have results for in Table 2: strange quarks and light quarks of masses (sets 1, 4 and 7), (sets 2 and 5) and the physical value, (sets 3, 6 and 8).

Instead of plotting the condensate results directly, the -axis in Fig. 4 is:

 Rl=−4ml⟨¯¯¯¯ψψl⟩(f2πM2π) (23)

for light quarks and

 Rs=−4ms⟨¯¯¯¯ψψs⟩(f2ηsM2ηs) (24)

for strange quarks. The values of the raw unsubtracted are determined directly from Table 2 using the , , and values given there.

The ratio is a good quantity to plot (and later to use in our fits) for a number of reasons:

• is a physical renormalisation-group invariant quantity as , up to discretisation errors, as is clear from the GMOR relation. The division by the square of the meson decay constant times its mass makes a dimensionless ratio which is convenient but it is also one that (from the GMOR relation) we expect to be close to 1.

• Using the ratio also reduces the effect of any slight mistuning of quark masses since the quark mass multiplied in the numerator cancels against the square of the meson mass in the denominator. The tuning of uses the decay constant, as described in Dowdall et al. (2011). This means that, by definition, does not contain discretisation errors that would mask the identification of the pieces that diverge as .

• Finally the ratio has reduced finite-volume effects over that in either the numerator or denominator. This is expected from the fact that chiral loop effects, which are sensitive to the volume, cancel in  Gasser and Leutwyler (1985); Jamin (2002). This is illustrated in Fig. 5 where we show results Toussaint (private communication) for pion mass, decay constant and (unsubtracted) light quark condensate as well as the ratio of Eq. (23) for ensembles with and and three different spatial volumes. The spatial volumes correspond to a spatial length in lattice units of 24, 32 and 40. The set with is our set 4 (see Table 1). For each quantity we plot the ratio of the value at to that at . It is clear from the plot that the finite volume dependence in each of , and is cancelled to a very high level of accuracy (0.1(1)% for set 4) in .

Figure 4 shows clearly the presence of a quadratic divergence with in the raw results. This is very ‘clean’ in our calculations because the form of the divergence is very constrained. Only a term of the form is allowed in for staggered quarks, i.e. no term of the form can appear. In the ratio this term takes the form where depends on the meson mass and decay constant. The HISQ formalism has very small discretisation errors, as is clear from the decay constant and meson mass results in Dowdall et al. (2011), and so there is little additional -dependence to confuse the analysis of the divergent pieces.

Because the power divergence is so dominant it is tempting to try to fit the unsubtracted results for to a very simple form: . This is in fact possible (it is important to include the error in the inverse lattice spacing when doing this since this is larger than the error in ) and we obtain which is the dashed line in lefthand plot of Fig. 4. We also obtain for with and for with , shown in the next two plots in Fig. 4. These fits are too naive to be useful, as we shall see below, because they miss out many important terms. Consequently the value and error of the intercept, , is unreliable for extracting a nonperturbative result for , especially in the quark case. However, the fits do illustrate that the ratio of slopes is that expected for a term that behaves as (although the simple fit does not allow for the running of the lattice bare quark mass with scale). The ratio of slopes between that for and for with is 3.2 which corresponds approximately to 5 (for the ratio of one power of the quark masses when the other power is cancelled by the square of the meson mass) divided by the ratio of the square of the decay constants from Table 2.

Figure 4 compares results for in which the tree-level piece of has been subtracted from the raw values of following Eq. 22. We take to be 2 GeV. These results are indicated by pluses. Now the slope in is much smaller since the most of the divergence has been removed. This makes the results more sensitive to the form of the remaining pieces of the divergence and the simple linear fits that were made to the unsubtracted data are no longer possible.

The perturbative contribution is very small for HISQ quarks and makes very little difference to the perturbative subtraction. The crosses show the results taking to be the full calculated perturbative subtraction through given in Eq. III.1. We have used  McNeile et al. (2010) for the value multiplying the one-loop coefficient, but the coefficient itself is so small that variations in scale for make no difference. The crosses are barely distinguishable from the pluses giving the tree-level subtracted numbers.

It is clear from Fig. 4 that there is still some divergence in left in after subtraction of the perturbative contribution through one-loop. This is not surprising since we know that will have higher order terms in . The challenge now is to fit the one-loop subtracted allowing for these higher order terms and thereby obtain the physical, nonperturbative, results for the strange and light quark condensates. We will fit both and simultaneously and use the known mass dependence of the unknown higher order terms to constrain them. At the same time we will allow for higher-order non-divergent mass-dependent terms from perturbation theory as well as physical, non-perturbative, dependence on the quark mass. Possible dependence on positive powers of , i.e. discretisation errors, must also be included.

#### iii.2.2 Determining a physical result from fitting

We now describe the full fit to the results that we use to determine the final physical values for () and () at the physical strange and light quark masses (where ). We take the following form for the ratio :

 Rq,0(a,amq)=R(q)NP,phys+δRPT+δRa2+δRχ+δRvol. (25)

are the raw results obtained from Table 2, is the final physical result in the scheme at 2 GeV. The terms represent fitted or known dependence on and . We use Bayesian techniques Lepage et al. (2002) to perform the fits so that we can add many higher order terms as part of each with constrained coefficients. This makes sure that the final error on is not underestimated by ignoring the existence of higher order corrections.

contains the known tree-level and one-loop perturbative results given in section III.1. In addition we include unknown higher-order terms. For the divergence these take the form:

 δRPT,div=an4αns(amq)2(afπ)2(aMπ)2 (26)

with the analogous term for the strange quark case, with the same . is a coefficient whose prior we take to be and we allow for , 3 and 4. Note that a prior width of 4.0 is conservative given the size of the corresponding coefficient at tree-level and one-loop. For the non-divergent pieces we take:

 δRPT,non−div=cn4αns(amq)4(afπ)2(aMπ)2 (27)

again with the analogous term for the strange quark case, with the same coefficient . We take , 3 and 4. in principle contains a sum of powers of up to . However, since is small for and our range of lattice spacings, these pieces are negligible and do not affect the fit and we simply take a prior on of 0.0(4.0). Again this is conservative given the results at tree-level and one-loop. For we use  McNeile et al. (2010) and discuss below the dependence of the results on changing to a different scale. takes values from 0.35 on the very coarse lattices to 0.26 on the fine ensembles.

allows for discretisation errors. We take the form

 δRa2=2∑i=1di(Λaπ)2i. (28)

Only even powers of appear in discretisation errors for staggered quarks and we take their scale to be set by 1 GeV. Since all tree level errors at are removed in the HISQ formalism we take the prior for to be i.e. 0.0(0.3). Higher order are given the prior 0.0(1.0). We include and terms, but have checked that higher order terms have very little effect. In addition we include a mass-dependent discretisation error in the form , giving a prior of 0.0(1.0). This allows for a number of effects, one of which could be mixing with a gluon condensate. This has negligible impact.

includes the valence and sea quark mass dependence that allows us to extrapolate to physical light quark masses and interpolate between the strange quark masses that we have to the physical strange quark mass. The chiral corrections to the GMOR relation were analysed in Gasser and Leutwyler (1985) (see also Jamin (2002)). The leading corrections are particularly simple because the chiral logarithms cancel to leave a correction proportional to . We allow for both and terms in our light quark mass fits by defining a chiral expansion parameter

 xl=M2π2(Λχ)2, (29)

with 1.0 GeV, and taking

 δRχ,val=2∑i=1g(l)ixil. (30)

We fit the , and results with this form taking the prior on the coefficients to be 0.0(2.0). This allows for a linear term of approximately the size expected in Jamin (2002). Higher order terms than have no effect.

The chiral expansion of Eq. (30) combines with the parameter in Eq. 25 to define the physical non-perturbative light condensate with its mass dependence. Since the GMOR relation is exact as we enforce this by taking the prior on to be 1.0000(5). Differences from 1 because of residual lattice finite volume effects up to 0.5% are allowed for as described in the paragraph above.

The data for is fit simultaneously with that for the light quarks, because they share parameters for the perturbative subtraction. However, we largely decouple the physical parameters because the strange quark is relatively far from the chiral limit and we would need a lot of parameters for a chiral expansion to connect light and strange quarks. Instead we allow a separate parameter for with the broad prior of 1.0(5). We take the same form for the valence mass dependence as in Eq. (30) where now

 xs=(M2ηs−(0.6893(12))2)2(Λχ)2, (31)

but this now simply allows for slight mistuning of the strange quark on some ensembles, and the fact that we have two values for the strange quark mass on sets 1 and 2. 0.6893 is the physical value for the mass determined in Dowdall et al. (2011). The priors on are the same as those on . Finite volume effects are expected to be completely negligible for because they are negligible for the components of .

The strange quark results in Fig. 4 show some very small sensitivity to the sea light quark masses if we compare results on sets 1 and 2 and sets 3 and 4. We therefore allow an additional linear dependence on the light quark mass in the sea in of the form

 δRχ,sea=2∑i=1ki(δmsea10ms,phys)i (32)

where

 δmsea=(2ml,sea+ms,sea)−(2ml,phys+ms,phys). (33)

We take  Bazavov et al. (2010b) and values determined from as in Dowdall et al. (2011). The prior for coefficient is taken as 0.0(1.0) which is conservative given the small effects observed in the results. We take to be common to both and for the light quarks. We note here also that the absence of chiral loop effects at this order means that staggered quark taste-changing effects are also absent. They can be handled, if necessary, with a sea-quark mass dependent term Dowdall et al. (2011). Including such a term here makes no difference to the physical result.

allows for remaining finite volume effects. These are small, as demonstrated in Fig. 5. They do produce a small systematic effect, however, because the lattice size in units of the pion mass, , is somewhat smaller on the lattices with smallest . We take:

 δRvol=ve−ML (34)

where is the pseudoscalar meson mass made of that quark ( for and for ) and is the linear extent of the lattice from Table 1. Coefficient is taken to have prior 0.0(0.2), consistent with Fig. 5.

Fitting and , , and , results simultaneously to the form in Eq. (25) readily produces good fits with for 18 degrees of freedom. The final fitted result for and (evaluated from and taken at the appropriate physical masses) is robust to the addition of higher order terms in the various corrections.

We take our final results from using for the scale of . The results do not change significantly as this is varied (although the fitted coefficients do change). Our fits return a substantial value for the coefficient of the power divergent term at , , of around 2.0, for . This is substantially larger than that seen at one-loop but not a particularly large value for a perturbative coefficient in general. It would simply imply that the small coefficient at one-loop for the HISQ action is not repeated at higher orders. We also find that the chiral correction to the GMOR for light quarks is substantial and negative (). This will be discussed further below.

The fit results are shown in Fig. 6. The data points (crosses) correspond to the lattice QCD results after subtraction of the perturbative contribution through (as in Fig. 4). The filled bands show the fitted curves when the full fitted perturbative contribution () is subtracted and masses and decay constants are set to the physical values corresponding to the quark and the light quark (for this we use ). These bands include the full error from the fit.

Our final physical results for are the key results from this paper.

 Rl,phys = −4ml⟨¯¯¯¯ψψl⟩¯¯¯¯¯¯¯¯MS(2GeV)(f2πM2π) Rs,phys = −4ms⟨¯¯¯¯ψψs⟩¯¯¯¯¯¯¯¯MS(2GeV)(f2ηsM2ηs) (35)

We find:

 Rs,phys = 0.574(86) Rl,phys = 0.985(5) Rs,physRl,phys = 0.583(84). (36)

The complete error budgets for , and their ratio are given in Table 4. The substantial 15% error that we have in reflects the difficulty of extracting a physical result from a power divergent quantity. For the error is 17 times better largely because the slope of the divergent piece is 15 times smaller. Errors in are dominated by errors from the lattice spacing and from fitting the remaining power divergent subtraction terms. There are also substantial errors from statistics and from tuning to the light and strange physical mass points. This is done by tuning the appropriate meson masses through the term in Eq. 30. This term depends on the lattice spacing through the definition of (Eq. 29) and (Eq. 31), because the meson masses appear in GeV units in these terms. The uncertainties in these terms then becomes correlated with the fit to the power divergence, increasing the uncertainty. For the power divergence is much less of an issue, but these same terms dominate the final error there as well.

The results for and can be converted to values of the condensate using the lattice result for = 0.1819(5) GeV and = 0.6893(12) GeV Dowdall et al. (2011), and experimental values for (0.1304(2) GeV) and (0.13498 GeV). We obtain:

 ms⟨¯¯¯¯ψψ⟩¯¯¯¯¯¯¯¯MSs(2GeV) = −2.26(34)×10−3GeV4 ml⟨¯¯¯¯ψψ⟩¯¯¯¯¯¯¯¯MSl(2GeV) = −7.63(4)×10−5GeV4. (37)

The ratio of the two values above is slightly more accurate than a naive combination, giving 29.6(4.3).

Using the precise determinations for light quark masses now available from lattice QCD we can finally obtain condensate values. We take = 92.2(1.3) MeV Davies et al. (2010a); McNeile et al. (2010) and  Bazavov et al. (2009); Bazavov et al. (2010b). These give:

 ⟨¯¯¯ss⟩¯¯¯¯¯¯¯¯MS(2GeV) = −0.0245(37)(3)GeV3 = −(290(15)MeV)3 ⟨¯ll⟩¯¯¯¯¯¯¯¯MS(2GeV) = −0.0227(1)(4)GeV3 (38) = −(283(2)MeV)3,

where the second error for each condensate in comes from the error in the quark masses.

For the ratio of strange to light condensate we have:

 ⟨¯¯¯ss⟩¯¯¯¯¯¯¯¯MS(2GeV)⟨¯ll⟩¯¯¯¯¯¯¯¯MS(2GeV)=1.08(16)(1), (39)

where the first error comes from and has the error budget given in Table 4 and the second error comes from the strange to light quark mass ratio.

#### iii.2.3 Approach of R to the chiral limit

The relationship of the light quark condensate to the chiral condensate is also important. is defined to have the value 1 from the GMOR relation in the chiral limit but the results of Eq. 36 indicate that it approaches this limit from below as the light quark mass is reduced. Supporting evidence for this is found by studying the quantity derived from the combination of condensates used by the HOTQCD collaboration in their study of finite temperature QCD Bazavov et al. (2012). We define by:

 Rδ=4mlf2πM2π⟨¯¯¯¯ψψl⟩−mlms⟨¯¯¯¯ψψs⟩1−mlms. (40)

The quadratic divergence with lattice spacing cancels between the two condensates because it is linear in the quark mass to all orders in perturbation theory. The non-divergent perturbative contributions proportional to the cube of the quark mass are completely negligible here, from the perturbative analysis in section III.1, and so we do not need to include them in making a physical quantity. can then simply be calculated from the raw data in Table 2.

A plot of against is shown in Fig. 7. In the limit on an infinite volume as does. can be determined more precisely than , however, because of the nonperturbative cancellation of the power divergence and it clearly approaches 1 from below. differs from by a term which is proportional to and to the difference between and . Both the dependence of on and the difference between and then contribute to the slope with seen in Fig. 7. We cannot separate them and therefore unambiguously identify the slope of with . We can however use this for a consistency check.

We fit to the simple form:

 Rδ = 1.000(1)+c1(1+c2a2+c3a4)mlms (41) + c4(mlms)2+c5e−mπL.

This allows for linear and quadratic terms in with discretisation errors. We take priors on , and to be 0.0(1.0) and prior on to be 0.0(0.3) consistent with behaviours. The final term allows for finite volume effects dependent on the combination . As for our fit to , we take the prior on to be 0.0(0.2) for consistency with Fig. 5.

The fit gives of 0.75 for 10 degrees of freedom and a physical slope, , of -0.51(4). This value is consistent with the difference between and 1 in Eq. 36, and indeed with the difference between and 1. This consistency between the results from and indicates that the difference between and (which would upset this consistency) cannot be large. This is indeed what we also find in Eq. 39.

#### iii.2.4 The chiral susceptibility

A further quantity that is of particular interest in studies of QCD at finite temperature is the chiral susceptibility for a quark of flavor :

 χf=∂∂mf(−⟨¯¯¯¯ψψf⟩). (42)

We give results here for the chiral susceptibility for zero temperature QCD to fill out the physical picture of the condensate. From differentiation of the path integral for the condensate it is clear that the chiral susceptibility is given by the flavor-singlet scalar correlator. It is convenient to split this into two contributions which we call and 111In Bazavov et al. (2012) these are called and .. comes from two scalar operators connected by quark lines, which is the flavor-nonsinglet scalar meson correlator. comes from two scalar operators connected only by gluons, in which the disconnected contribution is cancelled.

 χ = χq+χg (43) χq = 1Nt∑nTr[M−10nM−1n0] χg = −1N2tV(⟨(TrM−1)2⟩−⟨TrM−1⟩2).

The factors of number of tastes, , above are specific to naive/staggered quarks.