Determination of the quark condensate from heavy-light current-current correlators in full lattice QCD

We derive the Operator Product Expansion whose vacuum expectation value gives the time-moments of the pseudoscalar heavy-light current-current correlator up to and including terms in $\alpha_s^2$ multiplying $\langle\overline{\psi}\psi\rangle/M^3$ and terms in $\alpha_s$ multiplying $\langle \alpha_s G^2 \rangle/M^4$, where $M$ is the heavy-quark mass. Using lattice QCD results for heavy-strange correlators obtained for a variety of heavy quark masses on gluon field configurations including $u$, $d$ and $s$ quarks in the sea at three values of the lattice spacing, we are able to show that the contribution of the strange-quark condensate to the time-moments is very substantial. We use our lattice QCD time-moments and the OPE to determine a value for the condensate, fitting the 4th, 6th, 8th and 10th time-moments simultaneously. Our result, $\langle \overline{s}s \rangle^{\overline{\text{MS}}}(2 \text{GeV}) = -(296(11) \,\mathrm{MeV})^3$, agrees well with HPQCD's earlier, more direct, lattice QCD determination~\cite{McNeile:2012xh}. As well as confirming that the $s$ quark condensate is close in value to the light quark condensate, this demonstrates clearly the consistency of the Operator Product Expansion for fully nonperturbative calculations of matrix elements of short-distance operators in lattice QCD.

We derive the Operator Product Expansion whose vacuum expectation value gives the timemoments of the pseudoscalar heavy-light current-current correlator up to and including terms in α 2 s multiplying ψψ /M 3 and terms in αs multiplying αsG 2 /M 4 , where M is the heavy-quark mass. Using lattice QCD results for heavy-strange correlators obtained for a variety of heavy quark masses on gluon field configurations including u, d and s quarks in the sea at three values of the lattice spacing, we are able to show that the contribution of the strange-quark condensate to the time-moments is very substantial. We use our lattice QCD time-moments and the OPE to determine a value for the condensate, fitting the 4th, 6th, 8th and 10th time-moments simultaneously. Our result, ss MS (2GeV) = −(296(11) MeV) 3 , agrees well with HPQCD's earlier, more direct, lattice QCD determination [1]. As well as confirming that the s quark condensate is close in value to the light quark condensate, this demonstrates clearly the consistency of the Operator Product Expansion for fully nonperturbative calculations of matrix elements of short-distance operators in lattice QCD.

I. INTRODUCTION
The physics of the strong interaction clearly exhibits the features of spontaneously broken chiral symmetry induced by the condensation of quark-antiquark pairs in the vacuum. The fact that the vacuum expectation value of ψψ is non-zero is readily demonstrated in the fully nonperturbative approach to QCD provided by lattice QCD [2]. Calculating an accurate value for 0|ψψ|0 is, however, far from simple as we will discuss further below.
Each of the light quark flavours, u, d and s, gives a condensate, which could differ in value because of the different quark masses. The value of the condensate at zero quark mass (the chiral condensate) is an important parameter of low energy QCD [3]. Coefficients corresponding to derivatives of the chiral condensate with respect to quark mass then appear in higher order terms in the chiral expansion. Alternatively, and more simply, one can determine values for the condensates of specific quarks at their physical quark masses and these are then important input for analyses using QCD sum rules. For example, the values of both the s quark condensate and the light (u/d) quark condensate are key ingredients in the determination of |V us | from hadronic τ decay [4]. These condensates can be determined most accurately by lattice QCD calculations [1] and we give a new method for doing this here. We will apply the method to determine the strange quark condensate, although the approach could also be used for the light quark case.
One way in which the nonzero value of the quark condensate feeds into strong interaction physics is through its appearance in sub-leading ('nonperturbative') terms in the Operator Product Expansion (OPE) of shortdistance quantities [3,[5][6][7]. The OPE means that care must be taken to understand condensate contributions when determining such quantities in a fully nonperturbative calculation in lattice QCD. It also means that values for condensates can be determined from such calculations provided that sufficient accuracy is available to separate the different orders in the OPE. The determination of the same condensates from multiple short-distance quantities provides a strong test of the OPE approach.
The quark condensate is the lowest dimensional gaugeinvariant one, so might be expected to dominate the subleading terms of calculations of gauge-invariant quantities 1 . In fact quark condensate contributions are often rather small because ψψ will typically appear in a chirally-symmetric combination with the quark mass, m [3]. This then requires an additional inverse power of the high momentum parameter, Q, which defines the short-distance. For light quarks with m << Λ QCD the m/Q factor provides significant suppression. When a valence light quark is combined with a heavy quark, however, the spin of the heavy quark can be flipped with very little penalty, to compensate the light quark spin flip in ψψ and no factor of m appears. This implies that shortdistance quantities derived from heavy-light correlation functions should provide a good signal for the light quark condensate, if an accurate OPE can be constructed.
Here we demonstrate such an analysis using timemoments of pseudoscalar heavy-strange current-current correlators and an α 2 s -accurate OPE in the MS scheme that extends through NNLO in inverse powers of the heavy-quark mass (the short-distance scale here). The moments can be calculated very precisely in lattice QCD for a range of heavy quark masses up to the b. We use moments whose values are independent of the lattice spacing, up to discretisation effects that can be removed. The observed dependence on the heavy quark mass, matched to that expected from the OPE then allows us to obtain the strange quark condensate in the MS scheme. We fit the four lowest moments simultaneously for our final answer.
We can compare this method to earlier results from a more direct lattice QCD calculation [1]. The earlier calculation determined the expectation value of the trace of the quark propagator in lattice QCD. Naively this quantity would appear to be the quark condensate, but the mixing between the ψψ operator and the identity means that an OPE had also to be developed for this case. This was used to relate the lattice result to the strange quark condensate in the continuum in the MS scheme. The agreement between our new result and the earlier one is then a test of the OPE approach.
The paper is laid out as follows: Section II gives the theoretical background and the construction of the OPE, Section III describes the lattice calculation and Section IV the analysis that combines the two. Finally Section V compares our result for the condensate to earlier values and gives our conclusions.

II. THEORETICAL BACKGROUND
The time-moments of current-current correlators at zero spatial momentum can be related to q 2 -derivative moments of the corresponding polarisation functions. They are physical quantities, i.e. they do not depend on the ultraviolet regulator, and their values can be determined from the continuum and physical u/d quark mass limit of lattice QCD calculations. The HPQCD collaboration has pioneered the lattice calculation of such moments and their comparison with values, in the vector current case, derived from experimental results for the cross-section for e + e − → hadrons via a virtual photon as a function of centre-of-mass energy [12][13][14]. See also Refs. [15,16].
When the currents contain heavy quark fields, the low moments are perturbative, with the scale of α s being set by the heavy quark mass. Each moment is given by a power series in α s , multiplying an inverse power (depending on the moment number) of the quark mass.
For heavyonium currents (made of a heavy quark and antiquark) the perturbative series is known to high order (through α 3 s ) [17][18][19][20][21]. This then provides an accurate method for the determination of heavy quark masses by comparing the perturbation theory to nonperturbative results for the moments either from experiment for e + e − → hadrons (see, for example, Ref. [22]) or from lattice QCD [15,[23][24][25][26]. The results using experimental information are necessarily restricted to the vector (electromagnetic) current but lattice QCD results are available for currents with other spin-parity. The pseudoscalar current-current correlator is particularly useful in that case since, in lattice QCD formalisms with sufficient chiral symmetry, the quark mass times current is absolutely normalised as in continuum QCD, removing systematic uncertainties from current normalisation.
For the pseudoscalar case the moments are: Here M is the heavy quark mass, J 5 = ψ h γ 5 ψ h with ψ h the heavy quark field and the integral over x projects onto zero spatial momentum. The first four q 2 -derivative moments, 1, 2, 3 and 4 correspond to n = 4, 6, 8 and 10.
In perturbation theory M n is given by where g n is a power series in α s . An important advantage of the heavyonium currentcurrent correlator technique for determination of quark masses is the insensitivity to nonperturbative contributions coming from condensates. The radiative corrections to the heavy-quark vacuum polarisation and its derivatives close to q 2 = 0 [3,6] are sensitive to phase-space regions of quark momentum, p, and gluon momentum k where p 2 → 0 and k 2 → 0. Since heavy quarks are still highly virtual as p 2 → 0, it is only the k 2 → 0 region that is sensitive to long-distance physics. This generates a contribution to heavyonium current-current correlator moments given by the gluon condensate divided by four powers of the heavy quark mass. For the low moments used for the determination of the quark mass, the effect of this term on the moments is below 0.05% even when the heavy quark is a charm quark, and it has negligible impact [25].
In this work we study time-moments of heavy-light current-current correlators. We will see that nonperturbative condensate contributions, in this case coming from the light quark condensate, have about one hundred times more impact. This is largely because the leading quark condensate term appears with only three inverse powers of the heavy quark mass and no powers of α s /π (and, as discussed earlier, no powers of the light quark mass). Heavy-light moments then in fact present an opportunity for a determination of the light quark condensate. To do this we first need to determine an accurate OPE to compare to our lattice QCD results. This we do in the next section.

A. The OPE for heavy-light current-current corrrelators
Appendix B of Ref.
[1] discusses what is needed to derive the OPE for heavy-light current-current correlators. For ease of reference, we summarise the pieces that we will use here.
Consider moments of two pseudoscalar densities, J 5 = ψ h γ 5 ψ, composed of a heavy quark (mass M ) and a light quark (mass m M ) where the heavy quark fields are contracted with each other: As discussed above in the heavyonium case, the M 2 factor makes O (n) independent of the ultraviolet regulator provided n ≥ 4. This means that lattice and continuum calculations should agree in the limit of zero lattice spacing. Continuum results derived from lattice calculations can then be compared to continuum expressions derived from continuum QCD perturbation theory.
Operator O (n) is short-distance, dominated by length scales of order 1/M , provided the heavy-quark is sufficiently heavy and the light quarks have momenta small compared with M . Consequently the OPE implies that O (n) can be expressed in terms of a set of local operators in an effective theory, with cutoff scale Λ < M , and coefficient functions that depend only upon physics between scales Λ and M : 1 (Λ) is the unit operator and note that ψψ is not normalordered [27,28]. Working in the continuum we take the effective theory on the right-hand side of Eq. (2.5) to be QCD with an MS regulator and Λ = µ. Then we can express the right-hand side in terms of masses and couplings at scale µ, often chosen so that µ = M (µ). The coefficient functions c n , d n and e n are perturbative when M is large, and analytic in α s and m/M . They can be determined by perturbative matching of the matrix elements of the left-and right-hand sides of Eq. (2.5) between a given set of states. What we expect to see is that the matrix element of the left-hand side exhibits some infrared sensitivity that can be recognised as being part of the perturbative expression for one of the condensates on the right-hand side, requiring a particular value for the coefficient multiplying that condensate.
Here we can use vacuum matrix elements for the perturbative matching because the perturbation theory that we need in that case has already been done. We take the perturbative expressions for the vacuum matrix element of the left-hand-side of Eq. (2.5) as a function of x ≡ m/M from Refs. [29,30]. We can identify the condensate contributions from the non-analytic pieces that contain ln x multiplied by powers of x [31]. This can be done straightforwardly and unambiguously through multiple orders in α s and x because sufficiently accurate perturbative expressions for the light quark condensate are also known [32,33]. The gluon condensate does not appear until α 2 s x 4 ln x and so there are no issues with the separation of terms that can be identified with each condensate.
The starting point for the OPE derivation is (from Eq. (2.5)) where we take out a factor of n! for future convenience. The perturbative expansion of the left-hand side is given by the expressions in Refs. [29,30] for 3/(4π) 2 ×C p k for 2k + 2 = n.
The perturbative expansion through α 2 s (µ) for the vacuum matrix element of ψψ for a quark of mass m is given in Ref. [33]: where L m = ln(m/µ) and n f = n l +n m , with n l massless quarks and n m quarks with mass m ≡ m(µ) in the sea. This allows us to demonstrate simply how the OPE is derived using part of the perturbative expansion for one of the moments. We consider the zeroth and first order in α s expansions for the first (n=4) moment for the heavylight current-current correlator. The expressions given in Ref. [29] as a function of x = m/M can be expanded out to fourth order in x to give, through O(α s ) for the fourth The terms with ln x indicate infrared sensitivity of the series, which will be made explicit when these terms are replaced by the quark condensate. In Eq. (2.7) we can replace L m with ln x + ln(M/µ) to expose terms in ln x. This allows us to substitute ψψ for the leading x 3 ln x term on the second line of Eq. (2.8). In doing that we find that the x 3 ln 2 x at O(α s ) is automatically absorbed into ψψ . Likewise the x 4 ln x term at zeroth order can be replaced by a term proportional to x ψψ and this automatically absorbs the term at O(α s ) of the form x 4 ln 2 x.
Notice that none of this could be done if the perturbative calculation had set the light quark mass (and hence x) to zero. Our series for the fourth moment through x 4 then becomes The expansion clearly has the form expected 2 in Eq. (2.6) where the short-distance coefficients c n , d n , e n and their higher order counterparts are power series in α s with only polynomial dependence on x. We identify the series for d 4 in the square brackets of the third line. Notice that the leading term of d 4 is independent of x, i.e. ψψ appears without powers of the quark mass. It is also clear that the size of the coefficients appearing in d 4 , compared to those in c 4 , emphasise the condensate contribution to M 4 . This is what allows us to 'see' the effect of the condensate in the heavy-light moments so clearly.
The results in Refs. [29] and [33] in fact allow us to derive the c n and d n series through α 2 s and, respectively, x 4 and x. These series are given in Tables V and VI in Appendix A. We give results for n = 4, 6, 8 and 10. The relative size of the d n coefficients compared to c n grows with n. This is consistent with the expectation that longer-distance nonperturbative effects play an increasing rôle in the higher moments where there is a bigger contribution from larger t values.
The perturbation theory for heavy-light moments has been extended through α 3 s in [35] for massless quarks. Although this does not allow us to determine condensate contributions it does allow the leading x-independent term to be added to c n and we will test the impact of that in Section IV.
At O(α 2 s ) quarks in the sea start to contribute. It is then important to distinguish between the number of massless quarks, n l , the number of quarks of mass m, n m , and the number of heavy quarks, n h , active in the sea. Ref. [29] provides a Mathematica script giving their calculation at α 2 s with these pieces explicit. Although we derive here the OPE for heavy-light correlator moments for a light quark of mass m in terms of the condensate for that quark, the massless quarks in the sea will also condense. Because their masses have been set to zero in the perturbative expansion it is not possible to explicitly extract the contribution from their condensate. This is, however, a very small effect because these quark condensate contributions must appear multiplied by the light quark mass, as a consequence of chiral symmetry. Hence they can be neglected. This is discussed further in Section IV.
The perturbative expansion at α 2 s allows the term multiplying the gluon condensate to be determined at leading order, using the perturbative expansion for the gluon condensate from Ref. [32] We find e n = 1/(12π) in agreement with [3,34].
The OPE tells us that the behaviour of the heavy-light moments is very different from expectations from a naive application of the x-dependent perturbative expansion. We will be able to demonstrate that clearly in Section III from our lattice QCD results. For light quarks of small mass x will be small; it is less than 0.05 for the range of quark masses we use in our calculation. Powers of x and powers of x multiplying ln x then have relatively small impact. The OPE tells us, however, that in a fully nonperturbative scenario such as the real world or a lattice QCD calculation, these terms will be replaced by much larger terms coming from the light quark condensate. The coefficients with which the condensate terms appear, especially the fact that they enter at O(α 0 s ), means that the terms are clearly visible in the results for the heavylight moments. The OPE requires that the light quark condensate that enters here is the same matrix element as appears in the OPE for other operators and this consistency check of the OPE framework is an important one. Accurate nonperturbative results and an OPE that goes well beyond leading order are both required for this. The OPE derived here will be used in the next section along with our lattice QCD calculation to determine a value for the s quark condensate. I. Simulation parameters for the MILC n f = 2 + 1 gluon field ensembles that we use, labelled by set number in the first column. The lattice spacing values are given in units of a parameter known as r1, derived from the static quark potential, in the second column [36]. The physical value for r1 is 0.3133(23) fm [37]. Ls/a and Lt/s give the lattice dimensions. u0am sea l and u0am sea s give the sea u/d and s quark masses respectively in lattice units in the MILC convention, where u0 is the plaquette tadpole parameter. Set 1 will be referred to as 'fine', set 2 as 'superfine' and set 3 as 'ultrafine'. We use around 200 configurations from each ensemble and increase statistics by using 2 or 4 'random wall' time sources for propagators on each configuration.

III. THE LATTICE QCD CALCULATION
In the Section II we discussed how the OPE leads to time-moments of heavy-light current-current correlators containing important terms proportional to the light quark condensate. Here we describe how we calculate these time moments fully nonperturbatively in a lattice QCD calculation. We will focus on heavy-strange correlators and the determination of the s quark condensate, because that is computationally simpler than for u/d quarks. In fact, as we shall see, it is convenient to take a ratio of the heavy-strange moments to those of heavy-charm current-current correlators. Both are calculated in the same way in lattice QCD, simply with different quark masses being chosen in the solution of the Dirac equation on a given background gluon field to give the quark propagator, and different quark propagators combined to give the current-current correlator.
We work on ensembles (sets) of gluon field configurations that include the effect of u, d and s quarks in the sea with u and d quarks taken to have the same mass (denoted m l ). The configurations were generated by the MILC collaboration who used an O(a 2 ) improved discretisation of the gluon action and the improved staggered (asqtad) action for the quarks in the sea [36]. The main parameters of the ensembles used in this work are listed in Table I. The ensembles include three values of the lattice spacing, a, ranging from a ≈ 0.09 fm to a ≈ 0.045 fm. The sea s quark mass is close to the physical s quark mass but the u/d quark masses are not at their physical values, but instead at the heavier m l /m s = 0.2 point. This corresponds to a value for the pion mass, M π , around 300 MeV.
On each ensemble, we calculate quark propagators for the s quark, the c quark and then for a set of heavier quark masses heading towards that of the b quark. For these valence quarks we use the Highly Improved Staggered Quark (HISQ) action [39], and the valence quark  Table I) and then Columns 2 and 3 give the valence s and c quark masses in lattice units respectively. Column 4 gives the list of heavy quark masses used and Column 5 the corresponding values (taken from Ref. [38]) for the mass of the pseudoscalar heavyonium meson, η h , made from those heavy quark propagators. On Set 2, we use amc = 0.273, which is the tuned one. The data for amc = 0.28 was used to study the dependence on the charm quark mass, which was found to be negligible (see Section IV B).  (7) masses that we use are given in Table II. The tuning of the valence s and c quark masses to their physical values is discussed in Ref. [38]. The HISQ action was developed by HPQCD [39] to have very small discretisation errors and this has been demonstrated to be the case in many calculations (for example [24,40,41]). This reduction in discretisation errors means more accurate results across the board for a given value of the lattice spacing but this is particularly important in the heavy quark regime. There discretisation errors are controlled by the quark mass in lattice units (am) and the HISQ action has a higher reach in the quark mass for a given lattice spacing than actions that are less improved. This has allowed the HPQCD collaboration to initiate a programme of heavy quark physics using the HISQ action that stretches from c to b physics [24,41,42]. Values of ma that correspond to b quarks can only be reached on very fine lattices. With results for multiple masses at multiple lattice spacings, however, it is possible to fit a functional form to the results that also includes a functional form for discretisation effects and this can then be used to determine a value at the b quark mass in the continuum limit [24,42].
Here we are interested in correlation functions for pseudoscalar mesons, made by multiplying together the quark propagators discussed above. We calculate correlators made from a variety of combinations of quark masses, including one heavy quark whose mass we denote am h . We define the pseudoscalar current-current correlator at zero spatial momentum by where j 5 = ψ 1 γ 5 ψ 2 for the case here of two different quark flavours. We are working in the staggered formalism and so the γ 5 matrix is implemented through phasefactors [39]. Here we use the local pseudoscalar operator, i.e. the one with 'Goldstone' taste or, using staggered spin-taste notation, γ 5 ⊗ γ 5 . This means that, with the mass factors above, the current-current correlator is absolutely normalised, and no lattice current renormalisation factor is needed. The pseudoscalar correlator defined in (3.1) creates a meson at time 0 and destroys it at time t. In the large t limit the correlator is dominated by the ground-state meson with that valence quark content and spin-parity quantum numbers. Fitting the correlators and extracting the parameters of the large t behaviour enables the masses and decay constants of the ground-state pseudoscalar mesons to be determined. In Refs. [38,42] HPQCD showed that results at multiple heavy quark masses for multiple lattice spacings could be combined to determine the physical dependence on the heavy quark mass of the heavy-strange and heavy-charm meson masses and decay constants. This dependence could then be evaluated at the b quark mass to yield M Bs and M Bc for a comparison at the few MeV level with experiment, as well as f Bs and f Bc as inputs to flavour physics 3 .
Here we study the time-moments of the same pseudoscalar current-current correlators. These are defined by the lattice version of Eq. (2.3) where t respects the lattice periodicity so that [23] t In contrast to the results described in the paragraph above, the time-moments emphasise the small-t region of the correlator where the product of the two current operators is short-distance and an OPE is required, as described in Section II. Following Ref. [23], to reduce the discretisation errors in the time-moments we divide each moment by its tree level value, G n , calculated with the gluon fields set to the unit matrix. For the heavy-charm case G (0) n can be simply calculated on the same lattices as those used for the interacting gluon field ensembles. For the heavystrange case we must use space-time volumes that are a factor of 3 larger in each spatial direction in order to eliminate finite-volume effects in the free case [45]. We define the reduced moment R n as where m 0h can either be the bare mass or the tree-level pole mass of the heavy quark, which differ by discretisation effects [39]. Which is used here is irrelevant since this factor will be cancelled below. The continuum limit of the lattice QCD results gives the value for the moments derived in continuum QCD (3.5) The dimensionless function g n is the vacuum expectation value for the appropriate OPE for that moment and g n its tree-level component (see Eq. (A4)). If g n is derived in the MS scheme, then m h (µ) is the heavy quark mass in that scheme. For the heavy-strange case, as discussed in Section II A, the vacuum expectation value for the OPE has sizeable contributions from the s quark condensate that we want to determine. In order to emphasise these, and reduce systematic effects from the heavy quark masses, it is useful to take a ratio of the reduced moments R n for heavy-strange to those of heavy-charm. The heavy-charm case has only a small gluon condensate contribution, similar to that of heavyonium discussed in Section II and so makes a useful denominator. Another advantage of taking this ratio is that some other systematic errors are partly cancelled, for example those from discretisation effects and missing higher orders in perturbation theory.
For the ratio of heavy-strange moments to heavycharm we have, in the continuum limit, As explained in Secion II A, continuum QCD perturbation theory calculations through O(α 2 s ) as a function of light to heavy quark mass ratio (x) enable us to determine an accurate expansion of g n,hs given in Appendix A. g n,hs consists a leading-order (in powers of the heavy-quark mass) perturbative series c n with a subsidiary perturbative series d n multiplying the s quark condensate divided by the cube of the heavy quark mass and e n multiplying the gluon condensate divided by the fourth power of the heavy quark mass (Eq. (2.6)). g n,hs is the leading tree-level term in g n,hs , noting that no condensates appear when the gluon field is set to the unit matrix. In the heavy-charm case the same O(α 2 s ) heavy-light perturbation theory can be used [29,30] but now evaluated numerically for the specific charm to heavy quark mass ratios corresponding to each set of heavy-charm moments calculated. The values of the coefficients in g n,hc for each case used here are given in Table VII in Appendix A. The gluon condensate contribution to the heavy-charm moments will be discussed further in Section IV.
We will study the 4th, 6th, 8th and 10th timemoments. Table III gives our lattice results for the ratios for these moment numbers on each gluon field ensemble and for each heavy quark mass. Statistical uncertainties FIG. 1. Ratios of reduced moments for heavy-strange to heavy-charm calculated on the lattice (symbols) and the ratio of the appropriate perturbative series with no condensate contribution. The dotted and dashed lines show the ratio of the cn series for heavy-strange and the gn series for heavy-charm through order αs and α 2 s , respectively. The circles, squares, and triangles show the lattice results on ultrafine, superfine, and fine lattice ensembles respectively. are shown-they are very small. In addition our fits include the statistical correlations between the results on a given ensemble; these are not included in the Table. Figure 1 shows the lattice results for R n (hs)/R n (hc) at our three values of the lattice spacing. The x-axis is the cube of the inverse of m h ≡ m h (m h ), i.e. the heavy quark mass in the MS scheme at its own scale. The values for this are obtained from the values for M η h using the function derived in Ref. [24]. Figure 6 of that reference shows the relationship graphically. m h is proportional to M η h (up to very small corrections that we include) over the range of masses we use here from 2m c upwards.
Also plotted on Figure 1 is the appropriate ratio of perturbative series for the ratio of moments in which the condensate contributions are ignored. The dotted and dashed lines show increasing orders in the perturbative expansion. It is clear from this that the lattice QCD results do not agree with the perturbation theory when the condensate contributions are ignored. At small 1/m h the discrepancy is linear in 1/m 3 h also indicating that the discrepancy can be traced to a quark condensate contribution.
In the next section we discuss how we can fit the lattice QCD results to the perturbative expansion including condensate terms derived from the OPE, allowing for systematic uncertainties from higher order perturbative and nonperturbative effects. This allows us to determine a value for the s quark condensate quite accurately. It is clear from Figure 1 that this should be possible given the clear and significant contribution coming from that condensate.

A. Fit function
Our goal is to extract the strange-quark condensate from ratios of reduced moments R n (hs)/R n (hc) by comparing lattice-QCD results to theoretical expectations from the OPE and perturbation theory. To this end, we need to take into account discretization errors, scale uncertainties and mistunings of quark masses from lattice QCD and truncation of the OPE and truncation of the perturbation theory in the continuum.
Combining Eqs. (3.6) and (2.6) the continuum theoretical expectation for R n (hs)/R n (hc) is R n (hs) R n (hc) = g The output from the fit will be the fit parameter m s ss . We have modified g, c, d and e to absorb the appropriate g (0) factor and to allow for discretisation effects in the lattice determination of R n (hs) and R n (hc) and missing higher order terms in the perturbation theory, to be discussed below. We have expressed the condensate contribution in terms of m s ss for later convenience, using the identity m s (m h ) = x m h (since ratios of lattice quark masses in a particular quark formalism give the ratios of quark masses for a given continuum formalism at the same scale, up to discretisation effects [46]). As discussed above in Section III, m h is fixed from the η h mass. Uncertainties in the value of the lattice spacing (here set by r 1 , see Table I ) are fed in at this point. In the third line R n,hoOPE models our ignorance of higher dimension condensates and R n,mist , along with d mist , contain the effects of mistunings in sea and valence quark masses. Both of these will be described further below but we begin with terms on the second line. In Eq. (4.2)c andd are treated in a very similar way. We takec where the first term is the perturbative series for that piece of R n (hs) and the second term in square brackets allows for discretisation effects.
For the perturbative series we use the expressions tabulated in Tables V and VI, for c n,hs is the tree-level term in the expansion of M n , including the non-analytic terms in x that become the condensate in the interacting theory; these terms have very little effect since x is so small. α s is evaluated in the MS scheme at the scale µ = m h using n f = 3 evolution from the starting point α s (MS, n f = 3, µ = 5 GeV) = 0.2034(21) [24]. x is the ratio of s to heavy quark masses, both at the scale µ. Up to discretisation effects this is equal to the ratio of lattice bare quark masses x = am s /am h [46] obtained using the values from Table II. To include higher-order terms by treating the coefficients of these terms as fit parameters, we define, for i ≥ 3 with a prior of 0 ± 1 for eachc (i,j) n with i = 3 and 4, and 0±n 3 for eachd (i,j) n with i = 3. These priors are justified by comparing similar expressions for i < 3. Note that the expansion in Eq. (4.4a) is organized in powers of nx rather than x because a factor of n typically accompanies each power of x in the known expansions of c (i) n (x) for i < 3. The form of expansion, however, does not have any noticeable effects in our final result for the strangequark condensate.
The second term, in square brackets, of each of Eqs (4.3a) and (4.3b) allows for discretisation effects in the lattice QCD calculation of R n (hs). We allow separately for discretisation effects coming from the matrix element of the unit operator and the matrix element of ψψ since they have slightly different form. We view the discretisation effects in the d n term as being discretisation effects in the condensate and hence we do not allow the d lat coefficients to depend on n. We allow for discretisation effects that depend on aΛ and on am h , taking Λ = 0.3 GeV. Note that c We now discussg n,hc which allows for systematic uncertainties in R n (hc). The continuum perturbation theory for R n (hc) that we use is missing terms at α 3 s and above. However in the ratio R n (hs)/R n (hc) those terms can be reabsorbed into a single perturbative series with the missing α 3 s terms from R n (hs). Hence such terms are already taken into account by allowing for missing higher-order perturbative terms in Eqs (4.4a) and (4.4b) and we do not include any additional terms for g n,hc . This would also be true for discretisation effects except that the discretisation effects in R n (hc) can take a different form from those of R n (hs) because of the presence of the charm quark mass to set a higher scale for them. We takẽ g n,hc = g n,hc (4.5) g n,hc and g (0) n,hc are evaluated from the coefficients given in Table VII and the appropriate value of α s (µ). Since we expect the leading heavy quark discretization errors to cancel in the ratio of R n (hs) and R n (hc) we take g (0,l) n,lat = c 0,l n,lat . As above, since all tree-level errors are removed in the reduced moments, we take the priors of the remaining parameters to be O(α s ), namely 0.0(0.3).
In the ratio of heavy-strange moments to heavy-charm, both numerator and denominator have contributions from the gluon condensate. For the heavy-strange case, we have determined the leading order coefficient, e n , which is 1/(12π), independent of n. The heavy-charm case is similar to that of heavyonium, discussed in Section II in that the contribution comes from gluon propagators where k 2 → 0. Taking the results of Ref. [34] for the vacuum polarisation function for a valence quark and antiquark of different masses and expanding in powers of z = q 2 /m 2 h allows us to determine the gluon condensate contribution to the heavy-charm moments. For simplicity we also make an expansion in powers of 1/m c and collect terms in x −1 c and x 0 c where x c = m c /m h . We collect all contributions to the gluon condensate in one term with coefficientẽ n (see Eq. (4.2)). Theñ In order to treat m s ss and α s G 2 as fit parameters in Eq. (4.2) we must fix their scales, rather than allowing them to run with µ. To this end, we first rewrite them as Exploiting the perturbative expansion in Eqs. (2.7) and (2.10), we obtain (1 − 3L), (4.8b) where L = ln(µ/ν). The effect of the mixing of m s ss and α s G 2 with the identity or with each other, which generates the running in Eqs. (4.7a) and (4.7b), is very small since powers of the s quark mass appear. Any mixing of α s G 2 with the sea u/d quark condensate is completely negligible as a result. Then, setting µ = m h and ν = 2 GeV, we plug Eqs. (4.7a) and (4.7b) in Eq. (4.2). Finally, we treat m s ss (2 GeV) as a fit parameter with a prior of 0 ± 1 GeV 4 . For the gluon condensate we set the prior of αs π G 2 (2 GeV) to 0 ± 0.012 GeV 4 [47]. To allow for the effect of condensates of higher dimension than the ones covered by our OPE, we use (4.9) In the fit we set Λ = 0.3 GeV and take 0 ± 1 for the prior values of the dimensionless coefficients r (k) n,OP E . Note that we include a factor of n for each power of Λ because we see such a pattern for some fermionic condensates (see, for example, [28]). Our fit is insensitive to how many additional condensates we include beyond k = 8. Our results do not give a signal for any condensate other than the quark condensate.
Note that, as discussed in Section II A we can neglect the impact of a light quark condensate from the sea u/d quarks. This can appear at α 2 s but only in combination with the u/d quark mass and then divided by four powers of m h . Any such term would have negligible impact here.
(4.11) The values of δm sea /m s,phys for the ensembles that we use are tabulated in Table VI of Ref. [41]. Since we are using unphysial sea u/d quark masses the values of δm sea are non-zero and vary from 0.31 to 0.59. We also use, following Ref.
[1], a prior of 0.00(1) for r (k) n,mist . Eq. (4.10) incorporates the leading mistuning effects into our analysis. Since the strange-quark condensate itself can be affected by mistuning in the sea quark masses as well the valence strange-quark mass, we allow for this through the parameter d mist in Eq.  In this analysis we use M η s,phys = 0.6858 GeV [37], and following Ref.
[1] we use a prior of 0.0(1) for d (1) n,mist and 0.0(5) for d (2) mist . Using the fit function of Eq. (4.2) we proceed to fit the lattice QCD results of Table III to obtain a fitted value for m s ss (2 GeV).

B. Results
As described on Sec. III we have three sets of lattice gluon field configurations for which we have calculated R n (hs)/R n (hc) for n =4, 6, 8 and 10; results are tabulated on Table III. To fit these we use the fit function described in Section IV A and perform a combined fit as well as separate fits to the data for each n value. Figure 2 shows the lattice values and the band gives the fit results, extrapolated to zero lattice spacing and physical sea masses, from the combined fit. The combined fit has χ 2 /dof = 7/32. If we drop the quark condensate from the fit function the combined fit gives χ 2 /dof = 102/32 showing our inability to fit the lattice results without including this term. The fit yields a value for the fit parameter m s ss (2 GeV). To obtain the strange-quark condensate we divide m s ss (2 GeV) by m s (2 GeV; n f = 3) = 92.2(1.3) MeV [24], which is determined using the lattice-QCD ensembles used in this analysis and additional ensembles.
We ignore the correlation between the value of the strange-quark mass and our final result for m s ss (2 GeV) because the uncertainty in m s is considerably smaller than that in m s ss (2 GeV). Our value for the strange-quark condensate from a combined fit to (4.13) Figure 3 shows the stability of ss (2GeV) as we change the lattice data used in our work. It shows the results when a simultaneous fit or separate fits to the 4th, 6th, 8th and 10th time-moments are performed. The uncertainties in the separate fits are larger for n = 8 and 10 than for n = 4 and 6.
The dominant sources of error for ss (2GeV), obtained from the simultaneous fits to the 4th to 10th timemoments, are tabulated in Table. IV. The contributions from most sources are of similar size, about 1 to 2 percent. The dominant four uncertainties come, not surprisingly, from the gluon condensate, higher order condensates, higher order terms in perturbation theory and discretisation effects (a → 0 extrapolation).
To investigate the stability of our analysis we performed the following tests: • to investigate the impact of α 3 s corrections, we fixed c (3,0) n in Eq. (4.4a) to the known α 3 s correction to R n (hs) for massless light quarks [35], ignoring the corresponding (unknown) correction in R n (hc). We find a 1σ shift in ss under this test, but it should be emphasized that the α 3 s correction being added here is incomplete, and therefore not correct. Indeed we would expect some cancellation of α 3 s terms between R n (hs) and R n (hc) as we see at O(α s ) and O(α 2 s ) (see Tables V and VII). • to test the impact of varying the charm quark mass, we repeated the analysis using the alternative 'ultrafine' data with the mistuned mass value am c = 0.28 instead of those with the tuned am c = 0.273 (see Table II). The value of ss decreases by 0.2σ.
• to reduce the effect of lattice artifacts, we dropped the 'fine' data. The value of ss decreased by 0.6σ. • to test sensitivity to the presence of the gluon condensate, we repeated the analysis doubling the prior width of the gluon condensate. The value of ss then increases by 0.4σ.
These tests show that our result for the strange-quark condensate is stable under variations in the fit function and lattice data. It is also noteworthy that the posterior values of almost all fit parameters are within the 1σ width of their prior distributions. This implies that there is no tension between our choice of priors and the lattice data.

V. CONCLUSIONS
We have developed a new method for determining the light quark condensate using time-moments of heavylight current-current correlators calculated in lattice QCD. We fit lattice QCD results for the heavy-strange case as a function of heavy quark mass and lattice spacing to the expectation based on an OPE accurate through α 2 s and 1/m 4 h . We are able to extract a result for ss with a total uncertainty of 11%, or 3.7% in | ss | 1/3 . Our result, in the MS scheme at 2 GeV and with n f = 3 is The error budget is given in Table IV. The reasons that this accuracy is possible in this calculation are: • the light quark condensate is the leading (in powers of 1/m h ) condensate in the heavy-light moments. It appears at tree-level, suppressed by three powers of the heavy quark mass but with no suppression by the light quark mass.
• we have an α 2 s -accurate OPE so that the coefficients of both the leading 'perturbative' term (from the unit operator) and the coefficient of the quark condensate are under good perturbative control and the coefficient of the gluon condensate, at the next order in 1/m h , is known to leading order.
• we have accurate lattice QCD results for the timemoments at multiple values of the lattice spacing and multiple heavy quark masses. The HISQ action that we use has small discretisation errors (although they are visible here) and allows us to push to high quark masses. Using multiple heavy quark masses allows the identification of the quark condensate term from its functional dependence. Indeed we have shown that it is not possible to fit the results without including a quark condensate (Section IV B).
• we have used a ratio of heavy-strange to heavycharm correlator moments that removes systematic uncertainties from overall powers of the heavy quark mass that would otherwise appear. We also fit multiple moments simultaneously to improve the accuracy.
We can compare our result for the strange quark condensate to that obtained earlier by HPQCD from a completely independent method [1]. That method used lattice QCD results for the vacuum expectation value of the trace of the quark propagator (Tr(M −1 )) along with an O(α s )-accurate OPE for that case. In the OPE for the Tr(M −1 ) the relative behaviours of the unit operator and ψψ terms are rather different from the method introduced here. Here the short-distance scale is physical because it is set by m h and we can use lattice results at a variety of lattice spacing values and masses to pin down the condensate term. In the Tr(M −1 ) method the shortdistance scale is set by the lattice spacing. Then the term coming from the unit operator is suppressed, relative to ψψ , by a power of the light quark mass but it diverges, relatively, by two powers of the inverse lattice spacing as a → 0. This means that the most useful results are those on coarse lattices using improved actions for small discretisation effects and this is the approach taken in [1]. Results with multiple light quark masses were used to fit/constrain higher order terms but the strange quark condensate was less accurately determined than that of the light u/d quarks, because of the size of the contribution from the unit operator.
The value for the strange quark condensate obtained from the Tr(M −1 ) method is -(290(15) MeV) 3 [1]. Our new result in Eq. (5.1) agrees well with this and is more accurate. Note that the earlier result includes 4 flavours of sea quarks and here we include 3 flavours. We expect the impact of that change to be negligible, however, given that it produces a 0.2% change in m s from perturbation theory. This agreement between two very different calculations is strong validation of the OPE approach for short-distance quantities in fully nonperturbative QCD. The analysis here also 'closes the loop' on understanding both the small-t and large-t behaviour of the heavylight current-current correlators from lattice QCD, as discussed in Section III and as has already been achieved for heavyonium correlators [23]. Figure 4 compares our new result for | ss | 1/3 to the earlier result from Ref.
[1] for both | ss | 1/3 and | ll | 1/3 , where l has the average of the u and d quark masses.
Our new results confirm that ss and ll are close in value, with a 1σ preference for ss to be slightly larger in magnitude. Given that the chiral condensate (at zero quark mass) is smaller in magnitude than that of the light quark [1], it is clear that the magnitude of the condensate increases with quark mass for small quark mass. At heavy quark mass, however, the condensate magnitude must fall with quark mass [48]. Where, in the middle of this picture, the strange quark condensate sits is not yet completely clear. We believe that our new method can be used in future to explore/constrain further the dependence of the condensate on quark mass between light and strange.

ACKNOWLEDGMENTS
We are grateful to the MILC collaboration for the use of their configurations and to Matthias Steinhauser for useful discussions. The computing for this project used the Argonne Leadership Computing Facility at the Argonne National Laboratory and the Data Analytic system at the University of Cambridge, operated by the High Performance Computing Service on behalf of the STFC DiRAC HPC Facility. It was funded by a BIS national E-infrastructure capital grant and STFC capital and operations grants to DiRAC. We are grateful to the HPCS support staff for assistance. Funding for this work came from the Science and Technology Facilities Council.
Appendix A: Coefficients of the OPE for heavy-light current-current correlators The OPE leads to an expansion for the time-moments of heavy-light current-current correlators that can be organized as (repeating Eq. (2.6)) The short-distance coefficients c n , d n , e n and their higher order counterparts are power series in α s with only polynomial dependence on x. So, for example, with c   We use n l , n m , and n h to denote the number of massless, mass m and heavy mass M quarks in the sea, respectively. In the case of heavy-strange correlators calculated on gauge configurations with (2 + 1)-flavors that we use here (see Section III), we set n l = 2, n m = 1, and n h = 0.  [29,30,33]. With an O(α 2 s ) analysis we can only access e (0) n and we find e (0) n = 1/(12π), for all n. Note that in Eqs. (2.7) and (2.10) the gauge coupling α s corresponds to a theory with n l + n m active quarks, i.e., n h = 0, while the gauge coupling in Ref. [30] corresponds to a theory with n l + n m + n h active quarks. Therefore, when the matching is performed at scale µ,  to match the two couplings.
where we ignore, for now, condensate contributions. That issue is considered as part of the systematic error analysis in Section IV. We use the perturbative expansion given in Refs. [29,30] for unequal mass quarks, evaluating each of the coefficients at the mass ratio corresponding to that of the lattice charm and heavy quark masses. The values of the g (i) n at each of these ratios are tabulated in Table VII. Note that here, for 2+1 light flavours in the sea, n l = 3, n m = 0 and n h = 0.