Obsidian/Zotero/001_artiklid/takauranewapproachquarkmass2025.md

---

category: literaturenote citekey: takauranewapproachquarkmass2025 title: A new approach to quark mass determination using the gradient flow authors: "Takaura, Hiromasa; Harlander, Robert V.; Lange, Fabian" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.09537v1" zotero_key: WH8F5QGM zotero_storage: 5BZ7VJW8 collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Takaura, Hiromasa"

status: converted

A new approach to quark mass determination using the gradient flow

Hiromasa Takaura,^a Robert V. Harlander,^b Fabian Langec,d

E-mail: hiromasa.takaura@yukawa.kyoto-u.ac.jp,

robert.harlander@rwth-aachen.de, fabian.lange@physik.uzh.ch

Abstract: We propose a new method to determine quark masses using ratios of the vacuum-expectation values (VEVs) of flowed quark bilinear operators. They can be expressed as functions of the flow time t and the MS quark mass m, which can then be determined by matching with the corresponding lattice results. Motivated by this, we evaluate these VEVs perturbatively through next-to-leading order in the strong coupling. We provide the results as expansions in the limits of small and large m² t. To this end, we develop a new expansion technique based on the Laplace transform. Additionally, we present numerical results with the exact mass dependence over a wide range of m² t. We discuss the expected perturbative precision for the mass determination based on our next-to-leading order perturbative calculations, and possible non-perturbative corrections.

^aCenter for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

b Institute for Theoretical Particle Physics and Cosmology, RWTH Aachen University, Sommerfeldstraße 16, D-52056 Aachen, Germany

^cPhysik-Institut, Universität Zürich, Winterthurerstrasse 190, 8057 Zürich, Switzerland

^dPSI Center for Neutron and Muon Sciences, Forschungsstrasse 111, 5232 Villigen PSI, Switzerland

Contents

1		Quark mass determination using the gradient flow	1
2		Calculational method	3
	2.1	Theoretical framework	3
	2.2	Expressing the VEVs in terms of scalar integrals	4
	2.3	Expansion by Laplace transform	5
		2.3.1 Laplace transform	5
		2.3.2 One-loop example	6
		2.3.3 Two-loop examples	7
	2.4	Differential equations	9
	2.5	Numerical evaluation	11
	2.6	Small-flow-time expansion	11
3	Results		12
	3.1	Perturbative results	12
	3.2	Lattice observables and uncertainty estimates for quark mass extractions	14
4		Discussion of non-perturbative effects	19
5		Conclusions	21
A		List of loop integrals	21
B		Choosing v0 for the integral J17	24
C		Formulae	24
	C.1	Laplace transform	24
	C.2	Loop integrals	25
D		Perturbative series	26

1 Quark mass determination using the gradient flow

The gradient-flow formalism (GFF) [1–3] has been playing a crucial role in recent developments in lattice QCD. It allows for efficient scale setting, provides a gauge-invariant renormalization scheme which can be implemented on the lattice as well as in dimensional regularization, and enables one to assume the symmetries of the continuum even in lattice calculations.

In this paper, we propose a way to determine the quark masses with the GFF. Quark masses are fundamental parameters of QCD, and their precise determination is essential in contexts such as Higgs and flavor physics. Our proposal is to consider the ratios of the vacuum expectation values (VEVs) of the bilinear operators

$$S(t) \equiv \langle \bar{\chi}(t,x)\chi(t,x) \rangle$$ and $R(t) \equiv \langle \bar{\chi}(t,x) \overleftrightarrow{\mathcal{D}} \chi(t,x) \rangle$ (1.1)

where χ(t, x) and χ¯(t, x) denote the flowed quark fields. Matrix elements of S(t) and R(t) are individually ultraviolet (UV) divergent in general, even after renormalizing the QCD parameters. However, the remaining UV divergences can be fully absorbed into a flowed-quark wave function renormalization factor [3]. Thus, the divergences cancel when considering ratios of such matrix elements [4]. Therefore, the continuum limit in a lattice determination of these ratios exists. On the other hand, one can calculate these ratios in perturbation theory, where the result is expressed as a function of the flow time t, the strong coupling αs, and the quark mass m, renormalized in the modified minimal-subtraction scheme (MS). By matching the perturbative and lattice results, one should thus be able to determine the quark mass m.

The determination of the quark masses has been the subject of lattice calculations for many years. We refer to ref. [5] for an overview of available results. The main methodologies behind these studies are the regularization-independent momentum subtraction (RI-MOM) scheme [6] or its variant, the regularization-independent symmetric momentum subtraction (RI-SMOM) scheme [7, 8], the quark-current (or current-like) correlator method [9], or the minimal renormalon-subtracted (MRS) scheme [10]. Similar to the RI-(S)MOM schemes, our method is also based on quark bilinear operators. However, while in the RI-(S)MOM schemes three-point functions such as ⟨ψ¯(z1)ψψ¯ (x)ψ(z2)⟩ are calculated to study the renormalization of the bilinear operators, our proposal is based on the flowed one-point functions S(t) and R(t) of eq. (1.1). Moreover, the quantities of interest in our approach are manifestly gauge invariant and thus no gauge-variant operators are required for the renormalization or the short-distance (small-flow-time) expansion, which allows for a systematic discussion of nonperturbative corrections.

For the quark mass determination in this framework, one has to match lattice results with perturbative results. While S(t) including its mass dependence is known exactly at leading order in perturbation theory [3] and R(t) can easily be obtained in the same manner, higher-order results are only known in the limit of a small quark mass m. S(t) vanishes in the strictly massless limit; its linear mass term, together with the massless limit of R(t), was calculated at O(αs) in ref. [4] 1 and at O(α 2 s ) in ref. [11]. In addition, the leading mass correction of O(m² t) is available through O(α 2 s ) from ref. [12].

However, the size of the dimensionless parameter m² t can span a wide range, particularly when one is interested in heavy (i.e., charm and bottom) quarks. This can be seen as follows. Assuming a lattice spacing a, let us impose the conditions

$$a^2 \ll 8t \ll \Lambda_{\rm QCD}^{-2}$$ (1.2)

on the range of the flow time t, where ΛQCD ∼ 0.3 GeV is the mass scale of QCD derived from dimensional transmutation. The first condition is necessary for reasonable lattice

¹The result for S(t) was first presented in the arXiv versions v1 and v2 of ref. [4], whereas R(t) was presented in later arXiv versions and the published version.

simulations of eq. (1.1), the second one is imposed by perturbativity. Then, if a ⁻¹ = 4 GeV is assumed as a reference,2 0.1 ≪ 8m² c t ≪ 20 in the charm-quark case and 1.0 ≪ 8m² b t ≪ 200 in the bottom-quark case. Therefore, the known results valid for m² t ≪ 1 are not always sufficient.

In this paper, we evaluate the mass dependence of S(t) and R(t) at two-loop level, i.e., at O(αs) in perturbation theory. We provide the first several terms of the small- and largem² t expansions for these quantities.3 Additionally, we present the full mass dependence over a wide range of m² t, based on a numerical computation using ftint [14].

To expand loop integrals in either m² t or 1/(m² t), we develop a new method based on the Laplace transform. The standard technique to expand loop integrals under a certain hierarchy is known as "strategy of regions" [15]. Its application to the GFF in the smallflow-time limit has been described in ref. [16]. The opposite limit should be accessible as well with this method, following analogous considerations as described in ref. [17], albeit in a completely different context. The strategy employed here, on the other hand, is very symmetric in its application to small and large flow times.

The paper is organized as follows. In sec. 2 we explain our methods for the perturbative calculations. Sec. 3 presents our main results and a study of the perturbative behavior. We also provide an estimate of the achievable precision for the mass extraction using our nextto-leading order perturbative results. In sec. 4, we discuss the size of non-perturbative corrections and a way to properly suppress them for light quarks. Sec. 5 is devoted to our conclusions. Appendices A–C collect some detailed formulæ and discussion. Our full perturbative results are provided in app. D and also in ancillary files.

2 Calculational method

In sec. 2.1 we briefly summarize the perturbative approach to the GFF and then obtain the expressions for S(t) and R(t) in terms of scalar two-loop integrals in sec. 2.2. These loop integrals are evaluated using two approaches. On the one hand, we expand them for m² t ≪ 1 and m² t ≫ 1, leading to semi-analytic expressions in these two limits. The methodologies adopted here are explained in sec. 2.3. We also solve differential equations to evaluate some integrals, as explained in sec. 2.4. On the other hand, we evaluate the integrals numerically while keeping their full mass dependence, as explained in sec. 2.5. In sec. 2.6, we discuss how the first few orders in the expansions for small m² t can also be extracted from known results based on the small-flow-time expansion of the relevant bilinear operators.

2.1 Theoretical framework

The QCD action, before gauge fixing, is defined by

$$S_{\text{QCD}} = \int d^d x \left( \frac{1}{4g0^2} F{\mu\nu}^a(x) F_{\mu\nu}^a(x) + \bar{\psi}(x) (\not!!D + m) \psi \right)$$ (2.1)

²When one is interested in the bottom quark mass, even finer lattice spacings seem to be required. We consider a typical lattice spacing available today.

³These are updated results of our proceedings article [13], where the first few terms of the expansions for S(t) were presented.

where g⁰ denotes the bare coupling and

$$F{\mu\nu}(x) = \partial{\mu}A{\nu}(x) - \partial{\nu}A{\mu}(x) + [A{\mu}(x), A{\nu}(x)]$$ and $D{\mu} = \partial{\mu} + A{\mu}$ . (2.2)

We write the gauge field as A^µ = A^a µT a , where the normalization of the color generators T a in the fundamental representation is given by

$$\operatorname{Tr}\left(T^{a}T^{b}\right) = -\frac{1}{2}\delta^{ab}.\tag{2.3}$$

and

$$(T^a T^a)_{ij} = -C_F \delta^{ij}. (2.4)$$

with C^F = (N^{2 c} − 1)/(2Nc).

In the GFF, one defines generalizations of the gauge and the fermion fields in terms of the flow equations [1–3]

$$\partial{t}\chi(t,x) = (\mathcal{D}{\mu}\mathcal{D}{\mu} - \alpha{0}\partial{\mu}B{\mu}(t,x))\chi(t,x),$$

$$\partial{t}\bar{\chi}(t,x) = \bar{\chi}(t,x)\left(\overleftarrow{\mathcal{D}}{\mu}\overleftarrow{\mathcal{D}}{\mu} + \alpha{0}\partial{\mu}B{\mu}(t,x)\right),$$

$$\partial{t}B{\mu}(t,x) = \mathcal{D}{\nu}G{\nu\mu}(t,x) + \alpha{0}\mathcal{D}{\mu}\partial{\nu}B{\nu}(t,x),$$ (2.5)

and the boundary conditions

$$\chi(t=0,x) = \psi(x)$$ , $\bar{\chi}(t=0,x) = \bar{\psi}(x)$ , $B{\mu}(t=0,x) = A{\mu}(x)$ . (2.6)

The flowed covariant derivative is defined by

$$\mathcal{D}{\mu} = \partial{\mu} + B{\mu}, \qquad \overleftarrow{\mathcal{D}}{\mu} = \overleftarrow{\partial}{\mu} - B{\mu}, \qquad \mathcal{D}{\mu} = \partial{\mu} + [B_{\mu}, \cdot]$$ (2.7)

for the flowed quark, anti-quark, and gauge field, respectively. Also, for flowed quark bilinear operators, we use ←→D ^{µ = Dµ −}

$$\overleftrightarrow{\mathcal{D}}{\mu} = \mathcal{D}{\mu} - \overleftarrow{\mathcal{D}}_{\mu}. \tag{2.8}$$

α⁰ is a gauge fixing parameter which cancels in physical quantities. Throughout our calculation, we fix α⁰ = 1 and employ Feynman gauge for the gluon propagator.

2.2 Expressing the VEVs in terms of scalar integrals

Straightforward algebraic operations based on the perturbative method for solving the flow equations (2.5) [18] allow us to express S(t) and R(t) at NLO as

$$\begin{split} S|_{2-\text{loop}} &= g0^2 N\text{c} C_\text{F} \Big[ 8(d-2)I_1 - 32I_2 - 16I_3 - 16I_4 \ &- 8I_5 + 32I7 - 4(d-2)I{10} + 16I{11} \Big] \;, \ R|{2-\text{loop}} &= g0^2 N\text{c} C_\text{F} \Big[ 16J_1 + (4d-8)J_2 + 2(8-4d)J_3 + (16d-32)J_5 \ &- 64J_6 + 32J_7 - 32J_8 + 16J9 - 32J{10} + 64J{11} + 32J{12} \ &- 64J{13} + 16J{14} + (8-4d)J{15} - (8-4d)J{16} - (8-4d)J{17} \ &- (16-8d)J{18} + (16-8d)J{19} - 32J{20} \Big] \;. \end{split}$$

The two-loop scalar integrals Iⁿ and Jⁿ are listed in app. A, where we use the short-hand notation

$$\intp \equiv \int \frac{\mathrm{d}^d p}{(2\pi)^d} \quad \text{and} \quad \int{p,k} \equiv \int \frac{\mathrm{d}^d p}{(2\pi)^d (2\pi)^d} \,, \tag{2.10}$$

with d = 4−2ϵ the number of spacetime dimensions. We selectively made use of integrationby-parts identities to express some integrals in terms of others [11, 19, 20].

2.3 Expansion by Laplace transform

In this section, we explain how to expand the scalar integrals listed in app. A in the limits of small and large m² t. Our main technique, described in sec. 2.3.1, is applied to the integrals I1, I2, I3, I4, I5, I6, I8, I10, J1, J2, J3, and J17. The remaining integrals can either be expressed in terms of the already evaluated ones, up to the addition of m-independent integrals that can be directly evaluated, or are linked through differential equations (DEs), whose solutions are discussed in Sec. 2.4.

2.3.1 Laplace transform

It is convenient to write a general integral depending on m² and t as

$$I(m^2, t) = t^{-\alpha/2} \hat{I}(m^2 t),$$ (2.11)

where α is the mass dimension of the original integral. We now perform a Laplace transformation of the dimensionless integral ˆI(z) with respect to z ≡ m² t:

$$\tilde{I}(v) \equiv \int_0^\infty dz \, z^{-v-1} \hat{I}(z) \,. \tag{2.12}$$

The inverse transformation is given by

$$\hat{I}(z) = \frac{1}{2\pi i} \int_{-i\infty + v_0}^{i\infty + v_0} dv \, \tilde{I}(v) \, z^v, \qquad (2.13)$$

where v⁰ is a real number chosen such that the integral of eq. (2.12) converges for v = v0. If ˆI = O(z ^as ) for small z and I = O(z ⁻a^l ) for large z, one may choose −a^l < v⁰ < as. This also implies that ˜I(v) develops singularities in the v-plane which encode the structure of the expansions of I(z) in z and in 1/z, as we will now explain.

For z = m² t ≪ 1, we consider closing the integration contour of the inverse transform (2.13) in the right half of the v-plane, since z ^v → 0 as v → ∞ in this case. With the help of Cauchy's residue theorem we then arrive at

$$\hat{I}(z \ll 1) = -\sum{v{\text{sing}} > v0} \text{Res}[\tilde{I}(v)z^v]|{v=v_{\text{sing}}},$$ (2.14)

where we sum over the singularities of ˜I(v) located at vsing > v0. Res represents the residue, and the singularities are assumed to be poles. In turn, for z = m² t ≫ 1, we close the contour on the opposite side and obtain

$$\hat{I}(z \gg 1) = \sum{v{\text{sing}} < v0} \text{Res}[\tilde{I}(v)z^v]|{v = v_{\text{sing}}},$$ (2.15)

where we sum over the singularities of ˜I(v) located at vsing < v0. The factor z v in eqs. (2.14) and (2.15) allows us to distinguish individual orders of the expansions of ˆI(z) in z and 1/z, as we will see explicitly in the following examples. To obtain these expansions, we need to calculate the singularities of the integrand in eq. (2.13). Note that both the small- and large-m² t expansion can be obtained from the same quantity ˜I(v).

Our approach is analogous to the ideas found in refs. [21, 22], for instance. However, this is the first application of the inverse Laplace transform to higher loop integrals, employed for expanding them under certain hierarchy. Eqs. (2.12) and (2.13) are familiar from the standard Mellin-Barnes method for Feynman integrals [23, 24]. However, in this method transformations like eq. (2.12) are applied to rewrite every propagator by introducing auxiliary integration variables, whereas we introduce the single variable v by integrating over the kinematic scale z. Expansions in the context of standard Mellin-Barnes integrals have been discussed in refs. [25, 26].

2.3.2 One-loop example

Let us illustrate this method with a simple example. At the leading order in perturbation theory, S is given by

$$S|_{1-\text{loop}} = -4N_c \int_p \frac{m}{m^2 + p^2} e^{-2tp^2} \equiv -4N_c t^{1/2 - d/2} \hat{I}(m^2 t).$$ (2.16)

Applying the Laplace transform to the loop integrand yields a massless propagator with a power depending on v and we can trivially perform the loop integral to obtain

$$\tilde{I}(v) = t^{d/2} \int{p} e^{-2tp^{2}} \int{0}^{\infty} dz \frac{z^{-v-1/2}}{z + tp^{2}} = t^{d/2} \frac{\pi}{\cos(\pi v)} \int_{p} \left(\frac{1}{tp^{2}}\right)^{v + \frac{1}{2}} e^{-2tp^{2}}$$

$$= \frac{\pi}{\cos(\pi v)} \frac{1}{(4\pi)^{2-\epsilon}} \frac{\Gamma(3/2 - \epsilon - v)}{\Gamma(2 - \epsilon)} 2^{v - 3/2 + \epsilon}.$$ (2.17)

In this example, one may take v⁰ = 0, because each step of the calculation is well defined.

Now we study the singularities of ˜I. There are two origins of singularities. One is π/ cos(πv), a factor of the Laplace transform of the loop integrand. The other is Γ(3/2 − ϵ − v), which emerges after the loop integration. The singularities of the latter reflect infrared or UV divergences of the integral over p, given in the second line of the equation. In this case, we do not have UV divergences due to the factor e −2tp²

The positive (v > v⁰ = 0) singularities of ˜I(v, t) contribute to the small-(m² t) expansion. They are located at

$$v = n - \frac{1}{2}$$ and $v = n + \frac{1}{2} - \epsilon$ , $n \in \mathbb{N} = {1, 2, 3, ...}$ , (2.18)

which give

$$-\operatorname{Res}[\tilde{I}(v)z^{v}] = \begin{cases} \frac{z^{1/2}}{32\pi^{2}} + \mathcal{O}(\epsilon), & v = 1/2\ \frac{z^{3/2}}{16\pi^{2}} \left(\frac{1}{\epsilon} + 1 + 3\log 2 + \log \pi\right) + \mathcal{O}(\epsilon), & v = 3/2\ -\frac{z^{3/2 - \epsilon}}{16\pi^{2}} \left(\frac{1}{\epsilon} + 1 - \gamma_{E} + \log(4\pi)\right) + \mathcal{O}(\epsilon), & v = 3/2 - \epsilon\ \vdots \end{cases}$$ (2.19)

The 1/ϵ terms cancel between the v = k and v = k −ϵ contributions. As a result, we obtain the small-m² t expansion of eq. (2.16) as

$$S|_{1-\text{loop}} = -\frac{Nc}{8\pi^2} \frac{m}{t} \left[ 1 + 2m^2 t L{mt} + 4(m^2 t)^2 (-1 + L_{mt}) + \mathcal{O}((m^2 t)^3) \right] + \mathcal{O}(\epsilon), \quad (2.20)$$

where

$$L{mt} = \log(2m^2t) + \gamma{\rm E}. \tag{2.21}$$

The large-m² t expansion arises from the negative singularities located at

$$v = -n + \frac{1}{2}, \qquad n \in \mathbb{N}. \tag{2.22}$$

In this case there are no ϵ-shifted singularities. We obtain the large-m² t expansion

$$S|{1-\text{loop}} = -\frac{N{\text{c}}}{16\pi^2} \frac{1}{mt^2} \left[ 1 - \frac{1}{m^2 t} + \frac{3}{2} \frac{1}{(m^2 t)^2} + \mathcal{O}((m^2 t)^{-3}) \right] + \mathcal{O}(\epsilon). \tag{2.23}$$

These expansions agree with the result of ref. [16], obtained with the strategy-of-regions method [15]. In particular, the cancellation of the 1/ϵ terms has been discussed within that context. Since the exact result of eq. (2.16) can easily be obtained and reads

$$S|{1-\text{loop}} = -\frac{1}{(4\pi)^{d/2}} 4N{\text{c}} m^{d-1} e^{2m^2 t} \Gamma(1 - d/2, 2m^2 t), \qquad (2.24)$$

expanding it also allows us to verify the correctness of the new method.

2.3.3 Two-loop examples

Let us consider a few two-loop examples for the application of our method. Like in the one-loop example of sec. 2.3.2, we perform the Laplace transform on the integrand of the two-loop integrals, and subsequently evaluate the loop integrations. It turns out that we can set v⁰ = 0⁻ (slightly smaller than 0) in the inverse Laplace transform of eq. (2.13) for all integrals. We discuss a subtle point on this for J¹⁷ at the end of this section.

The integral I2. Our first example is the integral

$$I_2 = \int_0^t ds \int0^s ds' \int{p,k} \frac{mp^2}{k^2(p^2 + m^2)} e^{-(2t - s + s')p^2 - (s + s')k^2 - (s - s')(k - p)^2}.$$ (2.25)

Using the formulæ for the Laplace transform and the loop integrals of app. C, we find

$$\tilde{I}{2}(v) = -\frac{2^{-17/2 - v + 4\epsilon} \pi^{-3 + 2\epsilon}}{(-1 + \epsilon) \Gamma(2 - \epsilon)} \frac{\Gamma(\frac{5}{2} - v - \epsilon)}{\cos(\pi v)} \times \int{0}^{1} \int{0}^{1} du \, du' \, u^{\epsilon} (4 - u(1 - u')^{2})^{-\frac{3}{2} + v + \epsilon} {}{2}F_{1}(1, v - \frac{1}{2}; 2 - \epsilon, \frac{1}{4}u(1 - u')^{2}), \quad (2.26)$$

where we introduced the dimensionless flow-time integration variables s = tu and s ′ = su′ , and ²F¹ is the hypergeometric function. Since the integral over the flow-time parameters is finite, all singularities of ˜I2(v) in v are given by the factor Γ( ^{5 2} − v − ϵ)/ cos(πv), which has only simple poles at

$$v = \pm \left(n - \frac{1}{2}\right)$$ and $v = n + \frac{3}{2} - \epsilon$ , $n \in \mathbb{N}$ . (2.27)

The residues of ˜I² can thus be obtained by evaluating the u, u′ -integrals at these singular points. We use HypExp [27, 28] to expand the hypergeometric functions in ϵ, which allows us to obtain the coefficients for the 1/ϵ terms analytically. It turns out that they cancel in the final result of I2. Concerning the constant term of O(ϵ 0 ), we managed to obtain an analytical (numerical) result for the residues of the negative (positive) singularities.

The integral I5. Our second example is

$$I5 = \int{p,k} \frac{m}{k^2 (p^2 + m^2)((p-k)^2 + m^2)} e^{-tp^2 - tk^2 - t(k-p)^2}.$$ (2.28)

As opposed to the integral I² discussed above, the integral I⁵ contains two massive propagators. We thus first combine them using Feynman parametrization:

$$\frac{1}{(p^2+m^2)((p-k)^2+m^2)} = \int_0^1 \mathrm{d}x \frac{1}{[xp^2+(1-x)(k-p)^2+m^2]^2}. \tag{2.29}$$

After the Laplace transform and the loop integrations, the Feynman parameter integral can be evaluated analytically using Mathematica [29]. The only remaining integration is then over the parameter β (see eq. (C.5)) which we map to the interval y ∈ [0, 1] via β = y/(1−y). The y-integrand f(y) is rather complicated and not shown here. For positive v, it may develop singularities at y = 1 which one can write as

$$f(y) = \sum_{a=1}^{2} \frac{(1-y)^{-v-a\epsilon+1/2} g_a(y)}{\Gamma(\frac{1}{2}+v)\cos(\pi v)},$$ (2.30)

where the ga(y) are regular at y = 1. To identify the singularities of the Laplace transform, we perform the y-integral as

$$\int{0}^{1} dy f(y) = \frac{1}{\Gamma(\frac{1}{2} + v)\cos(\pi v)} \sum{a=1}^{2} \left[ \int{0}^{1} dy (1 - y)^{-v - a\epsilon + 1/2} \left[ g{a}(y) - \mathcal{T}{1,N} g{a}(y) \right] + \sum{k=0}^{N} \frac{(-1)^{k}}{k!} \frac{g{a}^{(k)}(1)}{v - k - 3/2 + a\epsilon} \right],$$ (2.31)

where T1,N ga(y) denotes the Taylor expansion of ga(y) around y = 1 up to power (1 − y) N , g (k) ^a (y) is the k th derivative of ga(y), and N an integer part of v − 3/2 (or v − 3/2 ≤ N < v − 1/2).

The remaining y integral is finite and can be calculated numerically. It contributes to the residues of the simple poles arising from the pre-factor 1/(Γ(1/2+v) cos (πv)) which are located at v = n/2 with n ∈ N. On the other hand, one can directly read off the residues (−1)^k g (k) ^a (1)/k! of the poles at v = k + 3/2 − aϵ arising from the second term in square brackets.

Similarly, for negative v, the integrand f(y) may develop singularities at y = 0, which can be treated in an analogous manner. Specifically, we find

$$\int{0}^{1} dy f(y) = \frac{1}{\Gamma(\frac{1}{2} + v) \cos(\pi v)} \left[ \int{0}^{1} dy \, y^{v+1/2} \left[ h(y) - \mathcal{T}{0,N} h(y) \right] + \sum{k=0}^{N} \frac{1}{k!} \frac{h^{(k)}(0)}{v+3/2+k} \right],$$ (2.32)

with the integer N = −v − 3/2. Note that, since the pre-factor only develops poles at positive v, the residues at negative v are given entirely by the Taylor expansion coefficients h (k) (0) of the function h(y) which can be easily evaluated analytically.

We managed to solve all the integrals Iⁿ and Jⁿ targeted by this method. For I7, I11, J13, and J15, we employ differential equations, as will be explained in sec. 2.4.

The integral J17. This integral involves an additional complication, in that its Laplace transform contains singularities at v = 0 and v = −ϵ, meaning that v⁰ in eq. (2.13) has to be chosen carefully. As argued in app. B, the choice −ϵ < v⁰ < 0 properly separates the small- and large-(m² t) expansions, assuming ϵ > 0.

Spurious 1/t² terms. Some of the Jⁿ develop terms ∼ 1/t² in the large-m² t expansion, even though the one-loop result of R(t) starts at order 1/t³ only:

$$R|_{1-\text{loop}} = -\frac{N_c}{8\pi^2} \frac{1}{m^2 t^3} \left( 1 - \frac{3}{2m^2 t} + \frac{3}{(m^2 t)^2} + \mathcal{O}((m^2 t)^{-3}) \right) + \mathcal{O}(\epsilon), \tag{2.33}$$

It turns out though that all 1/t² -terms indeed cancel among the Jⁿ when combining them to R(t) through eq. (2.9).

2.4 Differential equations

To calculate the remaining integrals I7, I11, and J13, it is more convenient to make use of differential equations (DEs) rather than applying the Laplace transform. In particular, I⁷ has a flow-time integral in addition to two massive propagators, which makes the calculation using the Laplace transform more complicated. It is easy to see that

$$(\partial_t - 2m^2)I_7 = m^2(I_5 - 2I_6),$$

$(\partialt - 2m^2)I{11} = -2m^2I_8,$ (2.34)

I⁶ is introduced for this purpose, although it does not directly contribute to S(t).

Let us consider a general DE of the form

$$(\partialt - 2m^2)I(t) = I{\text{ext}}(t), \tag{2.35}$$

where the expansions of Iext(t) are assumed to be known. By introducing i(t) as

$$I(t) \equiv e^{2m^2t}i(t), \tag{2.36}$$

we can derive the DE

$$i'(t) = e^{-2m^2t}I_{\text{ext}}(t)$$ (2.37)

for i(t).

We obtain the small-m² t expansion of I as follows. Using the known result of the small-m² t expansion of Iext(t), eq. (2.37) is given in the form

$$i'(t) = \sum_{k=0}^{\infty} (a_k m^{2k+3} t^{-1+k+2\epsilon} + b_k m^{2(k-\epsilon)+3} t^{-1+k+\epsilon} + c_k m^{2(k-2\epsilon)+3} t^{-1+k}).$$ (2.38)

Note that we keep the powers in t exact without expanding them in ϵ at this stage, while the coefficients ak, b^k and c^k are expanded in ϵ. Solving this DE, we obtain

$$i(t) = \sum_{k=0}^{\infty} \left( \frac{a_k}{k+2\epsilon} m^{2k+3} t^{k+2\epsilon} + \frac{b_k}{k+\epsilon} m^{2(k-\epsilon)+3} t^{k+\epsilon} + \frac{c_k}{k} m^{2(k-2\epsilon)+3} t^k \right) + \text{const.} \quad (2.39)$$

We note that the coefficient a⁰ is required up to O(ϵ) in order to obtain the result of I(t) = e2m² t i(t) up to O(ϵ 0 ). This means that the integrals on the right-hand side of eqs. (2.34) need to be evaluated through O(ϵ) using the Laplace transform method described in sec. 2.3.1. We also note that our calculation gives b⁰ = c⁰ = 0 for I⁷ and I11.

The constant term in eq. (2.39) is determined by

const. = $$i(0) = I(t = 0)$$ , (2.40)

with an explicit calculation of I(t = 0), which corresponds to a regular two-loop massive tadpole integral which can be calculated by standard methods.

Concerning the large-m² t expansion, we can write

$$i(t) = -\int{t}^{\infty} dt' i'(t') = -\int{t}^{\infty} dt' I_{\text{ext}}(t') e^{-2m^{2}t'}$$

$$= -e^{-2m^{2}t} \frac{I{\text{ext}}(t)}{2m^{2}} - \frac{1}{2m^{2}} \int{t}^{\infty} dt' I'_{\text{ext}}(t') e^{-2m^{2}t'}$$ (2.41)

where we used that t n e −2m^{2 t} → 0 for t → ∞ and arbitrary n. Furthermore, in the first step we assumed that i(t) = e−2m² t I(t) → 0 for t → ∞, which is certainly true as this limit does not introduce any singularities in a flow-time integral I(t). Repeatedly iterating the last step in eq. (2.41), we arrive at

$$I(t) = -\left(\frac{I{\text{ext}}(t)}{2m^2} + \frac{I'{\text{ext}}(t)}{(2m^2)^2} + \frac{I''_{\text{ext}}(t)}{(2m^2)^3} + \cdots\right). \tag{2.42}$$

Following this procedure, we obtain the small- and large-m² t expansions of I⁷ and I11. In an analogous manner, we use the DE

$$(\partialt - 2m^2)J{13} = m^2(J4 - 2J{12}), (2.43)$$

to determine the expansions of J13.

We remark that similar DEs can be used to obtain or check some of the results. For example, J¹⁵ and J¹⁸ are related via

$$(\partialt - 2m^2)J{18} = -2m^2J_{15}. (2.44)$$

2.5 Numerical evaluation

Since we expect the asymptotic expansions of the integrals to fail in the intermediate region, we also evaluate the integrals numerically on a one-dimensional grid in m² t. This is achieved with the help of the program ftint [14]. It uses Schwinger parametrization of the momentum integrations in order to map the integrals onto a multi-dimensional hypercube. The integrand becomes the product of polynomials in the Schwinger and flow-time parameters, raised to d-dependent powers. Such integrals can be evaluated algorithmically by sector decomposition [30, 31], which expresses them in terms of a Laurent series in d − 4 whose coefficients are convergent integrals which can be evaluated numerically. For the sector decomposition and the numerical integration, ftint relies on pySecDec [32–35]. Evaluation of all integrals to a relative precision of 10−⁸ or better on a grid of 200 points takes only a few hours on a currently average computer.

2.6 Small-flow-time expansion

In addition to a direct calculation of the mass expansion outlined above, one can also determine the expansion order-by-order by employing the small-flow-time expansion on the operator level [18] as was briefly mentioned in ref. [12]. This effectively amounts to an operator product expansion with t → 0. For S(t) it reads

$$\bar{\chi}(t,x)\chi(t,x) = \zeta{\rm S}^{(1)}(t)\frac{m}{t}\mathbb{1} + \zeta{\rm S}^{(3)}(t)m^3\mathbb{1} + \zeta_{\rm S}(t)\bar{\psi}(x)\psi(x) + \mathcal{O}(t)$$ (2.45)

up to mass dimension three [36] where the small-flow-time expansion separates local, unflowed operators from flow-time dependent matching coefficients. Eq. (2.45) holds both on the bare and the renormalized level. By taking the VEV we can obtain S(t) in terms of the matching coefficients ζ(t) and VEVs of unflowed operators.

Similarly, the small-flow-time expansion for R(t) reads

$$\bar{\chi}(t,x) \overleftrightarrow{\mathcal{D}} \chi(t,x) = \zeta_2^{(0)}(t) \frac{1}{t^2} \mathbb{1} + \zeta_2^{(2)}(t) \frac{m^2}{t} \mathbb{1}$$

$$+ \zeta{21}(t) \frac{1}{g^2} F{\mu\nu}^a(x) F{\mu\nu}^a(x) + \zeta{22}(t) \bar{\psi}(x) \overleftrightarrow{\mathcal{D}} \psi(x) + \zeta_{23}(t) m^4 \mathbb{1} + \mathcal{O}(t) \tag{2.46}$$

up to mass dimension four [37]. Using the equations of motion we can replace ψ¯(x) ←→ D ψ/ (x) = −2mψ¯(x)ψ(x) so that the same operator contributes to both expansions.

ζ S (t) in eq. (2.45) corresponds to the leading contribution to S(t) and was computed through NNLO in α^s in refs. [4, 11].4 ζ (3) S (t) is available through NNLO from ref. [38], and

⁴We again refer to the arXiv version v2 of ref. [4].

ζS(t) through NNLO from refs. [38–41]. The VEV of ψ¯(x)ψ(x) was computed through NNLO already a long time ago [42–44], but only for a theory with a single, massive quark flavor, i.e. without additional massless quarks.

Similarly, ζ (0) 2 (t) of eq. (2.46) corresponds to the leading contribution to R(t) and is available through NNLO from refs. [4, 11]. ζ (2) 2 (t) was computed through NNLO in ref. [37]. In the same reference, ζ23(t) was computed to the same order. ζ21(t) and ζ22(t) through NNLO are available from refs. [4, 37, 45]. Again, the VEV of g ⁻2F a µν(x)F a µν(x) was calculated through NNLO a long time ago [44, 46, 47], but also in a theory without additional massless quarks. Since the expansions of both ζ21(t) and the VEV of g ⁻2F a µν(x)F a µν(x) start at O(αs), they only contribute from NNLO onwards and are required to lower orders than the other contributions.

We use this strategy to extract analytic expressions up to O(αsm³ ) for S(t) and up to O(αsm⁴ ) for R(t), confirming our numerical expansions to this order. C S,≪1 1,0 , C S,≪1 1,1 , C R,≪1 1,0 , C R,≪1 1,1 , and C R,≪1 1,2 in app. D are obtained in this manner.

Extending the extraction to NNLO in α^s is more subtle: Eqs. (2.45) and (2.46) in this form implicitly assume a sum over the quark flavors for all quark operators. While up to NLO there is no mixing between heavy- and light-quark flavors and this subtlety is irrelevant, it becomes important to distinguish the flavors at NNLO. However, the published results of refs. [37, 38] do not contain sufficient information to fully disentangle the contributions. In addition, the full NNLO result for ⟨ψ¯(x)ψ(x)⟩ is missing as mentioned above.5

3 Results

3.1 Perturbative results

The renormalized mass and coupling in the MS scheme are given through

$$g_0 = \left(\frac{\mu e^{\gamma_E/2}}{\sqrt{4\pi}}\right)^{\epsilon} Z_g g(\mu^2),$$

$$m = Z_m \overline{m}(\mu), \tag{3.1}$$

where

$$Z_g = 1 + \mathcal{O}(\alpha_s),$$

$$Z_m = 1 - \frac{\alpha_s(\mu^2)}{4\pi} \frac{6C_F}{2\epsilon} + \mathcal{O}(\alpha_s^2),$$ (3.2)

with αs(µ 2 ) = g 2 (µ 2 )/(4π). In addition, the renormalized flowed quark fields in the MS scheme are defined by [3]

$$\chi{\rm R}(t,x) = Z{\chi}^{1/2} \chi(t,x)$$ and $\bar{\chi}{\rm R}(t,x) = Z{\chi}^{1/2} \bar{\chi}(t,x),$ (3.3)

where

$$Z{\chi} = 1 + \frac{\alpha{\rm s}}{4\pi} \frac{6C{\rm F}}{2\epsilon} + \mathcal{O}(\alpha{\rm s}^2). \tag{3.4}$$

⁵The required expressions could be easily obtained with current methods, but their calculation is beyond the scope of this paper.

Unless otherwise stated, m denotes m(µ) in what follows.

In this section, we analyze our two-loop results of the renormalized VEVs

$$\overline{S} := \langle \overline{\chi}{R}(t, x) \chi{R}(t, x) \rangle, \overline{R} := \langle \overline{\chi}{R}(t, x) \overleftrightarrow{\mathcal{D}} \chi{R}(t, x) \rangle$$ (3.5)

in the MS scheme. The explicit expansions in the limits of large and small m² t are presented in app. D and available in Mathematica readable format in ancillary files accompanying this manuscript. We checked that the leading terms in the small-m² t limit are correctly reproduced [4, 11, 12].6 In addition, we verified that they satisfy the renormalization group equations

$$\mu^{2} \frac{\mathrm{d}}{\mathrm{d}\mu^{2}} \langle \bar{\chi}{\mathrm{R}}(t,x) \chi{\mathrm{R}}(t,x) \rangle = \gamma{\bar{\chi}\chi} \langle \bar{\chi}{\mathrm{R}}(t,x) \chi_{\mathrm{R}}(t,x) \rangle,$$

$$\mu^{2} \frac{\mathrm{d}}{\mathrm{d}\mu^{2}} \left\langle \bar{\chi}{\mathrm{R}}(t,x) \overleftrightarrow{\not{D}} \chi{\mathrm{R}}(t,x) \right\rangle = \gamma{\bar{\chi}\chi} \left\langle \bar{\chi}{\mathrm{R}}(t,x) \overleftrightarrow{\not{D}} \chi_{\mathrm{R}}(t,x) \right\rangle, \tag{3.6}$$

where γχχ¯ = −3CFαs/(4π) + O(α 2 s ), which follow from eqs. (3.3) and (3.4).7

Another strong check of our results is the consistency of the expansions with the numerical evaluation of the integrals described in sec. 2.5, of course. We provide these numerical results for m² t ∈ [2−¹¹ , 2 ¹¹] in ancillary files as well. The explanation of the files is given in app. D. The left part of fig. 1 shows the perturbative coefficients of Sexp/Snum, where Sexp represents the small- and large-m² t expansion, whereas Snum is the numerical result retaining the full mass dependence. Schematically, the ratio is given by

$$\frac{\bar{S}{\text{exp}}}{\bar{S}{\text{num}}} = \frac{\bar{S}{\text{exp,LO}} + \bar{S}{\text{exp,NLO}}\alpha_s + \mathcal{O}(\alphas^2)}{\bar{S}{\text{num,LO}} + \bar{S}_{\text{num,NLO}}\alpha_s + \mathcal{O}(\alpha_s^2)}$$

$$= \frac{\bar{S}{\text{exp,LO}}}{\bar{S}{\text{num,LO}}} + \frac{\bar{S}{\text{exp,NLO}}\bar{S}{\text{num,LO}} - \bar{S}{\text{exp,LO}}\bar{S}{\text{num,NLO}}}{\bar{S}_{\text{num,LO}}^2}\alpha_s + \mathcal{O}(\alpha_s^2). \tag{3.7}$$

We observe the expected behavior, where the O(α s ) coefficient approaches one and the O(αs) coefficient approaches zero in the small- and large-m² t limits. Similarly, the right figure shows the analogous result for R. The small-m² t expansion given here is valid for m² t ≲ 0.2, and the large-m² t expansion is valid for m² t ≳ 6.

Let us make a remark on the large-m² t expansion. At one loop, the all-order series in 1/(m² t) reads [16]

$$S|{1-\text{loop}} = -\frac{N{\text{c}}}{16\pi^2} \frac{1}{\overline{m}t^2} \sum_{k=0}^{\infty} (-1)^k \frac{(k+1)!}{(2\overline{m}^2 t)^k}.$$ (3.8)

This is a sign-alternating divergent series, whose convergence radius is zero. Therefore, extending the large-m² t expansion to higher orders does not necessarily improve its approximation. The series gives the optimal approximation when it is truncated at an optimal order k∗, which is estimated as

$$k_* \sim 2\overline{m}^2 t - 2 \tag{3.9}$$

depending on the size of m² t one is interested in.

⁶We again refer to the arXiv version v2 of ref. [4] for S(t).

⁷Note also that µ ²dαs/dµ ² = −ϵα^s + O(α 2 ^s ), following from eq. (3.1).

Figure 1. The perturbative coefficients of the ratio of S (left) and those of R (right), obtained using two approaches, i.e., expansions and numerical evaluation, as a function of m² t, cf. eq. (3.7). The lines in the region m² t < 1 correspond to the small-m² t expansion (blue and purple lines) and those in m² t > 1 to the large-m² t expansion (red and orange lines). Since the O(αs) coefficient also depends on µ/m in addition to m² t, we show results with a few choices: µ/m = 100 (purple solid), 5 (purple dashed), 1 (purple dotted), and µ/m = 1 (orange solid), 2 (orange dashed), 1/2 (orange dotted).

3.2 Lattice observables and uncertainty estimates for quark mass extractions

From the discussion at the beginning of sec. 3 it is clear that ratios of VEVs of bilinear operators are finite and renormalization group invariant. Specifically, we consider three observables

$$r_a(m) = \frac{S(t,m)}{R(t,m)}, \qquad r_b(m) = \frac{R(t,m)}{R(t,m=0)}, \qquad r_c(m) = m \frac{\mathrm{d}}{\mathrm{d}m} \left(\frac{S(t,m)}{R(t,m)}\right)$$ (3.10)

as examples. All of them should be accessible on the lattice (the m = 0 limit may have to be approximated by a very small quark mass though). Their perturbative expressions can be obtained from our expansions of S¯ and R¯. They can also be obtained from the purely numerical results of sec. 2.5 (r^c requires a numerical derivative which is unproblematic though).

To see how these observables behave in the small- and large-m² t limits, we show the first few terms of the expansions in these limits. We express the perturbative coefficients in terms of

$$Lt \equiv \log(2\mu^2 t) + \gamma{\rm E}$$ and $L_m \equiv \log(\mu^2/\overline{m}^2)$ . (3.11)

For r^a the perturbative series reads

$$r{a}(\overline{m}) = \overline{m}t \left{ 1 + C{F}\alpha{s}(\mu^{2})(0.276671 + 0.238732L{t}) \right.$$

$$+ \overline{m}^{2}t \left[ 2 - 2L{m} + 2L{t} + C{F}\alpha{s}(\mu^{2})(0.645693 - 0.30742L{m} - 0.95493L{m}^{2} + 2.69474L{t} + 0.477465L{m}L{t} + 0.477465L{t}^{2}) \right]$$

$$+ \mathcal{O}((\overline{m}^{2}t)^{2}) \left. \right} + \mathcal{O}(\alpha_{s}^{2})$$

$$(3.12)$$

for m² t ≪ 1, and

$$r{a}(\overline{m}) = \overline{m}t \left{ 0.5 + C{F}\alpha{s}(\mu^{2})(-0.0554033 + 0.0596831L{m} + 0.0596831L{t}) + \frac{1}{\overline{m}^{2}t} \left[ 0.25 + C{F}\alpha{s}(\mu^{2})(-0.0465134 - 0.0596831L{m}) \right] + \mathcal{O}((\overline{m}^{2}t)^{-2}) \right} + \mathcal{O}(\alpha_{s}^{2})$$ (3.13)

for m² t ≫ 1. The perturbative series for rc(m) can be straightforwardly obtained from this expression. For r^b the perturbative series reads

$$\begin{split} rb(\overline{m}) &= 1 - \overline{m}^2 t \bigg{ 2 + C{\rm F} \alpha_{\rm s}(\mu^2) (0.589536 + 0.95493 L_t) \ &+ \overline{m}^2 t \Big[ - 4 L_m + 4 Lt \ &+ C{\rm F} \alpha_{\rm s}(\mu^2) (-1.60014 - 1.41061 L_m - 1.90986 L_m^2 \ &+ 3.32047 L_t + 1.90986 Lt^2) \Big] \ &+ \mathcal{O}((\overline{m}^2 t)^2) \bigg} + \mathcal{O}(\alpha{\rm s}^2) \end{split} \tag{3.14}$$

for m² t ≪ 1, and

$$r_b(\overline{m}) = \frac{1}{\overline{m}^2 t} \left{ 1 - C_F \alpha_s(\mu^2) (0.0753887 + 0.358099 L_m + 0.119366 L_t) + \frac{1}{\overline{m}^2 t} \left[ -1.5 + C_F \alpha_s(\mu^2) (0.3496 + 1.01461 L_m + 0.417782 L_t) \right] + \mathcal{O}((\overline{m}^2 t)^{-2}) \right} + \mathcal{O}(\alpha_s^2)$$ (3.15)

for m² t ≫ 1. Higher orders in m² t and in 1/(m² t) can be easily obtained from the results presented in app. D. In addition, the analytical expressions for these coefficients are available. The above equations are presented only to illustrate the structure.

The dependence of these results on the renormalization scale µ is formally of higher orders in perturbation theory. Nevertheless, in order to study their numerical behavior, we need to choose a value for µ. In the small-t region, it is natural to set µ ² ∼ 1/t, since this drives αs(µ 2 ) into the perturbative region and also eliminates potentially large logarithms. In the large-t region, however, we need to prevent µ from entering the non-perturbative regime. Unless stated otherwise, our central choice for µ which interpolates between these two requirements is

$$\mu_{\rm int}(t) \equiv \sqrt{\mu_t^2 + \overline{m}^2(\overline{m})},$$ (3.16)

with

$$\mu_t \equiv \frac{e^{-\gamma_E/2}}{\sqrt{2t}} \,. \tag{3.17}$$

We do not claim that eq. (3.16) is an optimal choice for general 1/ √ t and m. For the scope of this manuscript we prefer a simple prescription over a dedicated study of the optimal renormalization scale.

As input parameters, we adopt [48]

$$\begin{split} &\alpha{\mathrm{s}}(M{\mathrm{Z}}^2) = 0.1180 \,, \ &\overline{m}{\mathrm{s}}(2 \; \mathrm{GeV}) = 93.5 \; \mathrm{MeV} \,, \ &\overline{m}{\mathrm{c}}(\overline{m}{\mathrm{c}}) = 1.2730 \; \mathrm{GeV} \,, \ &\overline{m}{\mathrm{b}}(\overline{m}_{\mathrm{b}}) = 4.183 \; \mathrm{GeV} \,. \end{split} \tag{3.18}$$

Using the five-loop beta function [49–51] and the four-loop decoupling relations for the coupling [52], we obtain the running coupling at our central scale choice. Similarly, the running mass at the central scale is obtained based on the four-loop anomalous dimension of the MS mass [53].

The perturbative results of the observables defined in eq. (3.10) for the strange quark are shown in fig. 2 where we subtract one from r^b to suppress non-perturbative effects as discussed in sec. 4. Here we use the small-m² t expansion, as m² s t ≲ 0.005 in the range of these figures. Therefore, rather than using µint(t) of eq. (3.16), we simply set the central renormalization scale to µ = µ^t , defined by eq. (3.17). We confirmed that the results agree with the numerical results retaining full mass dependence at larger t in the figures. The black dashed lines show the leading term in the small-m² t expansion at O(αs). The perturbative uncertainties are indicated by the bands whose width is determined by

$$\Delta = \max_{\mu = \mu_t \text{ or } 2\mut} r - \min{\mu = \mu_t \text{ or } 2\mu_t} r.$$ (3.19)

We do not consider µ = µt/2 here, as it would extend into the non-perturbative region. For this variation of µ, we use the one-loop mass anomalous dimension at O(α 0 s ) and the one-loop beta-function and the two-loop mass anomalous dimension at O(αs). The band is then symmetrized around the central value for ra, obtained by setting µ = µ^t

The analogous results for the charm and bottom quark are shown in fig. 3 and 4, respectively. Here we use the numerical results retaining full mass dependence rather than the expansions, since the validity of the latter does not extend to the full range in t. Also, we only consider r^a and r^b in this case for reasons that will be discussed in sec. 4. 8 As the central renormalization scale, we choose µ = µint(t) in this case, defined by eq. (3.16). For reference, we remark that the two terms µ^t and m contributing to µint are equal at t ≈ 0.173 GeV−² for the charm quark, and at t ≈ 0.016 GeV−² for the bottom quark. The perturbative uncertainties are estimated in analogy to eq. (3.19), but using µ = µint(t)/2 and 2µint(t).

Let us now study the expected accuracy of a mass determination using the three observables ra, rb, and rc. To estimate this accuracy, let us for the moment assume that an infinitely precise lattice measurement of a quantity r is available, denoted by rexp, and that the theoretical prediction r(m) of this quantity depends on the mass m. Neglecting also theoretical uncertainties, we could determine the physical value m as the solution of

$$r(\overline{m}) = r_{\exp}. \tag{3.20}$$

⁸For the charm we also subtract the trivial one from r^b for better visibility.

Figure 2. Perturbative results for the strange quark, with error bands estimated by varying the renormalization scale by a factor of two around its central value µt, see eq. (3.17). The number of flavors is set to n^f = 3. The O(αs) result with the leading term in the m² t-expansion is shown by the black dashed line for reference.

Figure 3. Perturbative results for the charm quark, with error bands estimated by varying the renormalization scale by a factor of two around its central value µint(t), see eq. (3.16). The number of flavors is set to n^f = 4. The small-m² t expansion at O(αs) is shown by the black dashed line for reference. The large-m² t expansions have no relevance in this t-range.

However, if the theoretical prediction carries an uncertainty δr, then all one can say is that, to some level of confidence,

$$r(\overline{m}) \in [r{\exp} - \delta r, r{\exp} + \delta r].$$ (3.21)

Figure 4. Perturbative results for the bottom quark, with error bands estimated by varying the renormalization scale by a factor of two around its central value µint(t), see eq. (3.16). The number of flavors is set to n^f = 5. The small- and large-m² t expansions at O(αs) are shown by the black dashed lines for reference. The small-m² t expansion for the left figure is not shown as it is valid only in a tiny region.

Thus, we can determine m only up to ±∆m, where

$$r(\overline{m} \pm \Delta \overline{m}) \approx r(\overline{m}) \pm r'(\overline{m}) \Delta \overline{m} \in [r{\exp} - \delta r, r{\exp} + \delta r],$$ (3.22)

and therefore

$$\Delta \overline{m} = \left| \frac{\delta r}{r'(\overline{m})} \right| . \tag{3.23}$$

In fig. 5, we study this quantity for the observables ra, rb, and rc. We only take into account the perturbative uncertainty in δr though; eventually, one would have to take into account the uncertainty from the lattice input, of course. The figures show the expected precision ∆m/m (solid black), the relative (perturbative) uncertainty δr/r (PT uncertainty), and the dimensionless derivative mr′ (m)/r × 0.01 (where 0.01 is multiplied for the sake of visibility in the figure) (Sensitivity to mass), as a function of t, fixing m at the values given in eq. (3.18).

In the case of the strange and charm quarks, where the small-m² t is more relevant, we see the tendency that ∆m/m decreases towards smaller values of t, because the perturbative uncertainty decreases in this case. For the strange quark, the effect of the larger perturbative uncertainty for r^b compared to r^c is partly compensated by the larger derivative, leading to a similar expected accuracy on m. In the case of the bottom quark, the figure on r^b indicates the potential to extract its mass value at 1−2 % already at this perturbative order. We note that these are crude estimates, based on the scale variation around our central scale choices, eq. (3.16) or (3.17). Therefore they are understood as order-of-magnitude estimates.

Besides the perturbative uncertainty due to missing higher orders, we also have to take into account the uncertainties from the lattice determinations for a full uncertainty estimate. In addition, there could also be non-perturbative corrections which we discuss in the next section.

Figure 5. Estimate of the expected accuracy on mass determination at the current perturbative order [O(αs)]. See the text for details.

4 Discussion of non-perturbative effects

So far we have neglected non-perturbative effects which will be particularly relevant for light quarks with m ≪ ΛQCD. For m² t ≪ 1 we can estimate them based on the small-flow-time expansion discussed in sec. 2.6. Taking the expressions from eqs. (2.45) and (2.46), we can split the observables into perturbative and non-perturbative contributions and find

$$\bar{S} = \bar{S}{PT} + \zeta{S}(t) \left[ \langle \bar{\psi}\psi \rangle{R} - \langle \bar{\psi}\psi \rangle{R,PT} \right] + \mathcal{O}(t)$$ (4.1)

and

$$\bar{R} = \bar{R}{PT} + \zeta{21}(t) \frac{1}{g^2} \left[ \langle F{\mu\nu}^a F{\mu\nu}^a \rangle{R} - \langle F{\mu\nu}^a F{\mu\nu}^a \rangle{R,PT} \right]

• \zeta{22}(t) \left[ \langle \bar{\psi} \overleftrightarrow{p} \psi \rangle{R} - \langle \bar{\psi} \overleftrightarrow{p} \psi \rangle_{R,PT} \right] + \mathcal{O}(t),$$ (4.2)

where S¯ PT and R¯ PT represent the perturbative results calculated in this paper up to O(αs). The differences [⟨O⟩^R − ⟨O⟩R,PT] between the full, exact evaluation of the VEV and its perturbative contributions, which are non-zero for m ̸= 0, single out the non-perturbative effects.

For S¯ in eq. (4.1), the non-perturbative contribution, ⟨ψψ¯ ⟩^R −⟨ψψ¯ ⟩R,PT, is of O(Λ³ QCD) and can practically be more significant than the leading perturbative contribution S¯|PT = O(m/t). One can see this as follows. For the strange quark we are in the regime where m² ≪ Λ 2 QCD ≪ (8t) ⁻¹ ≪ a ⁻² and the ratio of the non-perturbative contribution to the perturbative contribution, ^Λ3 QCD m/t , can exceed one due to

$$\frac{\Lambda{\rm QCD}}{\overline{m}} \frac{\Lambda{\rm QCD}^2 a^2}{8} \ll \frac{\Lambda{\rm QCD}^3}{\overline{m}/t} \ll \frac{\Lambda{\rm QCD}}{\overline{m}} \frac{1}{8}.$$ (4.3)

Similarly, the leading finite mass effect in R¯|PT, that is O(m2/t), can be subdominant compared to the non-perturbative effects of O(Λ⁴ QCD). This is not optimal for a mass determination.

Compared to the simple ratios r^a and r^b in eq. (3.10), observables like

$$r_c(m) = m \frac{\mathrm{d}}{\mathrm{d}m} \left( \frac{S(t,m)}{R(t,m)} \right)$$ and $1 - r_b(m) = 1 - \frac{R(t,m)}{R(t,m=0)}$ (4.4)

could be more suitable, because they suppress non-perturbative effects. Assuming that the non-perturbative effects consist of possible power terms of m and ΛQCD, and for instance, as for ψψ¯ ,

$$\langle \bar{\psi}\psi\rangle{\rm R} - \langle \bar{\psi}\psi\rangle{\rm R,PT} = c0\Lambda{\rm QCD}^3 + c1m\Lambda{\rm QCD}^2 + c2m^2\Lambda{\rm QCD},$$ (4.5)

where c0,1,² are constants, one can see that the largest non-perturbative contribution, c0Λ 3 QCD, which is independent of m, disappears due to the derivative in rc. Hence, the leading non-perturbative effect reduces to O(mΛ 2 QCD), whose magnitude relative to the perturbative effect m/t reduces to O(tΛ 2 QCD), which is much smaller than one and can be well controlled. As for 1 − R/(R|m=0), the leading non-perturbative effect is O(mΛ 3 QCD), which is smaller than that found in R/(R|m=0), but still can be large compared with the perturbative effect m2/t. Based on these observations, we conclude that r^c is the most suitable quantity for light quark mass determination among ra,b,c and 1 − rb.

For heavy quarks, the non-perturbative effects are properly suppressed already in r^a and rb, as long as m² t ≪ 1, by applying a similar argument as above for Λ 2 QCD ≪ m² ≪ (8t) ⁻¹ ≪ a −2 . This is why r^c is studied only for the strange quark.

For m² t ≫ 1 we know little about non-perturbative effects. This study would be worth pursuing in the future.

5 Conclusions

We proposed a method to determine quark masses using the gradient flow by considering ratios of VEVs of flowed quark-bilinear operators. They are scheme and gauge invariant and can be calculated both on the lattice and in perturbation theory. We evaluated these VEVs perturbatively up to two loops, i.e. O(αs), including finite mass effects. Compared to the already known perturbative results in the small-m² t limit, we found that the higherorder terms of the small-m² t expansion, the large-m² t expansion, and full mass dependence are crucial for studying particularly charm and bottom quarks. The perturbative behavior of the VEV ratios exhibit good perturbative behavior. We also estimated the expected accuracy achieved in the mass determination from the perturbative viewpoint.

We presented the small- and large-m² t expansions of these VEVs, as well as fully numerical perturbative series retaining the full mass dependence. To calculate the expansions, we developed a method based on the Laplace transform rather than the often employed technique called "strategy of regions." The Laplace transform with respect to a combination of the external scales m and t is considered, whose singularities and residues give expansions of loop integrals. Although the Laplace transforms have been used in somewhat analogous manners in the literature, to our understanding, this is the first application to expanding integrals beyond one loop in terms of external scales.

In order to improve the accuracy of the quark mass determination and to be competitive with existing methods, it may be required to compute higher orders in the perturbative series. Even though the technology for such calculations is available (see, e.g. refs. [11, 14]), they are beyond the scope of the current paper.

Finally, the method proposed in this paper requires lattice results of the flowed quark bilinear operators for determining the quark mass. These quantities have been simulated in refs. [54, 55] for light quarks (u, d, s). It will be interesting to study which accuracy can be achieved already with these simulations following the proposed method and what is required for a competitive extraction.

Acknowledgments

We thank Robert Mason for checking parts of our numerical results. H.T. thanks Yukinari Sumino, Hiroshi Suzuki, and especially Masakiyo Kitazawa for helpful discussions. The work of H.T. was supported by JSPS KAKENHI Grant Numbers JP19K14711 and JP23K13110. H.T. is the Yukawa Research Fellow supported by Yukawa Memorial Foundation. The work of R.V.H. was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 - TRR 257. The work of F.L. was supported by the Swiss National Science Foundation (SNSF) under contract TMSGI2\_211209.

A List of loop integrals

In this appendix we list all of the integrals appearing in our calculation of S or R described in sec. 2. Some of the following integrals do not contribute directly, but are introduced to facilitate the evaluation of other integrals.

The scalar integrals relevant for S are

$$I{1} = \int{0}^{t} ds \int_{p,k} \frac{m}{k^{2}(p^{2} + m^{2})} e^{-2tp^{2} - 2sk^{2}},$$

$$I{2} = \int{0}^{t} ds \int{0}^{s} ds' \int{p,k} \frac{mp^{2}}{k^{2}(p^{2} + m^{2})} e^{-(2t - s + s')p^{2} - (s + s')k^{2} - (s - s')(k - p)^{2}},$$

$$I{3} = \int{0}^{t} ds \int_{p,k} \frac{m}{k^{2}(p^{2} + m^{2})} e^{-(2t - s)p^{2} - sk^{2} - s(k - p)^{2}},$$

$$I{4} = \int{0}^{t} ds \int{0}^{t} ds' \int{p,k} \frac{m(p - k)^{2}}{k^{2}((p - k)^{2} + m^{2})} e^{-(2t - s - s')p^{2} - (s + s')k^{2} - (s + s')(k - p)^{2}},$$

$$I{5} = \int{p,k} \frac{m}{k^{2}(p^{2} + m^{2})((p - k)^{2} + m^{2})} e^{-tp^{2} - tk^{2} - t(k - p)^{2}},$$

$$I{6} = \int{0}^{t} ds \int_{p,k} \frac{m}{k^{2}((p - k)^{2} + m^{2})} e^{-(2t - s)p^{2} - sk^{2} - s(k - p)^{2}},$$

$$I{7} = \int{0}^{t} ds \int_{p,k} \frac{m^{3}}{k^{2}(p^{2} + m^{2})((p - k)^{2} + m^{2})} e^{-(2t - s)p^{2} - sk^{2} - s(k - p)^{2}},$$

$$I{8} = \int{p,k} \frac{m}{k^{2}(p^{2} + m^{2})^{2}} e^{-2tp^{2}},$$

$$I{9} = \int{p,k} \frac{m}{k^{2}(p^{2} + m^{2})^{2}} e^{-2tp^{2}},$$

$$I{10} = \int{p,k} \frac{m}{k^{2}(p^{2} + m^{2})^{2}((p - k)^{2} + m^{2})} e^{-2tp^{2}},$$

$$I{11} = \int{p,k} \frac{m^{3}}{k^{2}(p^{2} + m^{2})^{2}((p - k)^{2} + m^{2})} e^{-2tp^{2}}$$

$$(A.1)$$

and the scalar integrals relevant for R are

$$\begin{split} J_1 &= \int0^t \mathrm{d}s \int{p_{,k}} \frac{(p-k)^2}{k^2((p-k)^2+m^2)} \mathrm{e}^{-(t-s)p^2-(t+s)k^2-(t+s)(p-k)^2} \ J2 &= \int{p_{,k}} \frac{1}{(p^2+m^2)((p-k)^2+m^2)} \mathrm{e}^{-tp^2-tk^2-t(p-k)^2} \ J3 &= \int{p_{,k}} \frac{1}{k^2(p^2+m^2)} \mathrm{e}^{-tp^2-tk^2-t(p-k)^2} \ J4 &= \int{p_{,k}} \frac{m^2}{k^2(p^2+m^2)((p-k)^2+m^2)} \mathrm{e}^{-tp^2-tk^2-t(p-k)^2} \ J_5 &= \int0^t \mathrm{d}s \int{p_{,k}} \frac{p^2}{k^2(p^2+m^2)} \mathrm{e}^{-2tp^2-2sk^2} \ J_6 &= \int0^t \mathrm{d}s \int{p_{,k}} \frac{p^2}{k^2(p^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J_7 &= \int0^t \mathrm{d}s \int{p_{,k}} \frac{p^2}{k^2(p^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J_8 &= \int0^t \mathrm{d}s \int{p_{,k}} \frac{p^2}{k^2(p^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J_9 &= \int0^t \mathrm{d}s \int{p{,k}} \frac{(p-k)^2}{k^2((p-k)^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J{10} &= \int0^t \mathrm{d}s \int{p{,k}} \frac{p^2}{k^2(p^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J{11} &= \int0^t \mathrm{d}s \int{p{,k}} \frac{m^2}{k^2(p^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J{12} &= \int0^t \mathrm{d}s \int{p{,k}} \frac{m^2}{k^2(p^2+m^2)} \mathrm{e}^{-(2t-s)p^2-sk^2-s(k-p)^2} \ J{13} &= \int0^t \mathrm{d}s \int{p{,k}} \frac{m^2}{k^2(p^2+m^2)((p-k)^2+m^2)} \mathrm{e}^{-2tp^2} \ J{15} &= \int{p{,k}} \frac{1}{k^2(p^2+m^2)} \mathrm{e}^{-2tp^2} \ J{16} &= \int{p{,k}} \frac{1}{k^2(p^2+m^2)} \mathrm{e}^{-2tp^2} \ J{17} &= \int{p{,k}} \frac{1}{k^2((p-k)^2+m^2)} \mathrm{e}^{-2tp^2} \ J{18} &= \int{p{,k}} \frac{m^2}{k^2((p-k)^2+m^2)} \mathrm{e}^{-2tp^2} \ J{19} &= \int{p{,k}} \frac{m^2}{k^2((p-k)^2+m^2)} \mathrm{e}^{-2tp^2} \ J{19} &= \int{p{,k}} \frac{m^2}{k^2((p-k)^2+m^2)} \mathrm{e}^{-2tp^2} \ J{19} &= \int{p{,k}} \frac{m^2}{k^2((p-k)^2+m^2)} \mathrm{e}^{-2tp^2} \ J{19} &= \int{p_{,k}} \frac{m^2}{k^2((p-k)^2+m^2)^2((p-k)^2+m^2)} \mathrm{e}^{-2tp^2}. \end{split}$$

In the above expressions, m denotes the bare mass and k denotes the gluon loop momentum.

B Choosing v⁰ for the integral J¹⁷

In the Laplace transform of J17, we find singularities at v = −ϵ and v = 0. Both of them give O(1/t² ) terms upon the ϵ-expansion, because a singularity at v corresponds to a term of the form m2^v t −2+v+2ϵ . We need to clarify whether these singularities contribute to either the small- or large-m² t expansion. (This issue is related to which value v⁰ should be set to. In principle, one can determine this by searching for v = v⁰ where every step of the calculation is valid, but this is not always practical.)

First, it turns out that the v = 0 singularity contributes to the small-m² t expansion. This is because the residue at v = 0 is identical to the original loop integral with the mass in the propagator sent to zero. Secondly, as for v = −ϵ, it is useful to consider a bound on J¹⁷ in d dimensions,

$$|J{17}| < \int{p,k} \frac{1}{k^2 (p-k)^2} e^{-2tp^2}$$

$$= C \int_p \frac{1}{(p^2)^{(4-d)/2}} e^{-2tp^2}$$

$$= C't^{2-d} = C't^{-2+2\epsilon}.$$ (B.1)

In the first inequality, we used 1/[(p − k) ² + m² ] < 1/(p − k) 2 . The bound (B.1) is valid as long as 2 < d < 4; the k-integral is UV divergent for d ≥ 4 and infrared divergent for d ≤ 2. Therefore, we assume this range of d and thus ϵ > 0. C and C ′ are positive constants. From this bound, we can see that the v = −ϵ singularity, which gives O(m−2^ϵ t −2+ϵ ), contributes to the large-m² t expansion; if it contributed to the small-m² t expansion, it would exceed the bound for small t. To summarize, we should choose −ϵ < v⁰ < 0 in this calculation.

In the calculation of J17, using eq. (C.4) after the change of momenta p → p + k and then k → −k, leads to a singular expression. One can circumvent this by introducing "the s-parameter in eq. (C.4)" as a regulator and then sending it to zero.

C Formulae

We collect formulae concerning the Laplace transform and loop integrals.

C.1 Laplace transform

$$\int_0^\infty dz \, z^{-v-1/2} \frac{1}{P^2 + z} = \frac{\pi}{\cos(\pi v)} \frac{1}{(P^2)^{v+1/2}},$$

$$\int_0^\infty dz \, z^{-v-1/2} \frac{1}{(P^2 + z)^2} = \frac{1}{2} (1 + 2v) \frac{\pi}{\cos(\pi v)} \frac{1}{(P^2)^{v+3/2}}.$$ (C.1)

P denotes a linear combination of (dimensionless) loop momenta.

C.2 Loop integrals

$$\begin{split} &\int{p} \frac{1}{(p^{2})^{a}} \mathrm{e}^{-sp^{2}} = \frac{1}{(4\pi)^{d/2}} \frac{\Gamma(d/2 - a)}{\Gamma(d/2)} s^{a - d/2} \quad [2], \ &\int{p,k} \frac{\mathrm{e}^{-sp^{2} - uk^{2} - v(k - p)^{2}}}{(k^{2})^{a}} = \frac{1}{(4\pi)^{d}} \frac{\Gamma(d/2 - a)}{\Gamma(d/2)} \frac{1}{(s + v)^{a}} \frac{1}{(uv + su + sv)^{d/2 - a}}, \quad (C.3) \ &\int{p,k} \frac{1}{(k^{2})^{a}(p^{2})^{b}} \mathrm{e}^{-sp^{2} - uk^{2} - v(k - p)^{2}} \ &= \frac{1}{(4\pi)^{d}} \frac{\Gamma(d/2 - b)}{\Gamma(d/2)\Gamma(a)} \frac{1}{(u + v)^{b}} \frac{1}{(uv + su + sv)^{d/2 - b}} \ &\times \left[ \left( \frac{uv + su + sv}{s + v} \right)^{a} \frac{\Gamma(a)\Gamma(d/2 - a - b)}{\Gamma(d/2 - b)} {}{2}F{1}(a; b; 1 + a + b - d/2; \frac{uv + su + sv}{(s + v)(u + v)}) \right. \ &+ \frac{1}{\Gamma(b)} (u + v)^{a} \left( \frac{(s + v)(u + v)}{uv + su + sv} \right)^{b - d/2} \Gamma(a + b - d/2)\Gamma(d/2 - a) \ &\times {}{2}F{1}(\frac{d}{2} - a, \frac{d}{2} - b; \frac{d}{2} + 1 - a - b; \frac{uv + su + sv}{(s + v)(u + v)}) \right], \quad (C.4) \ &\int{k,p} \frac{1}{(k^{2})^{a}[xp^{2} + (1 - x)(k - p)^{2}]^{b}} \mathrm{e}^{-sp^{2} - uk^{2} - v(k - p)^{2}} \ &= \frac{1}{(4\pi)^{d}} \frac{\Gamma(d/2 - a)}{\Gamma(b)\Gamma(d/2)} \ &\times \int_{0}^{\infty} \mathrm{d}\beta \beta^{b - 1} \frac{1}{(s + v + \beta)^{a}} \frac{1}{[u(v + \beta(1 - x)) + (s + \beta x)u + (s + \beta x)(v + \beta(1 - x))]^{d/2 - a}}. \quad (C.5) \end{split}$$

In the last equation, x corresponds to a Fyenmal integral parameter (0 ≤ x ≤ 1). ²F1(a, b; c; z) denotes the hypergeometric function, which is defined by

$${2}F{1}(a,b;c;z) \equiv \sum{n=0}^{\infty} \frac{(a){n}(b){n}}{(c){n}} \frac{z^{n}}{n!}$$ (C.6)

for |z| < 1 and a non-positive integer c. (...)ⁿ is the Pochhammer symbol (a)ⁿ = Γ(a + n)/Γ(a). In our calculation, we often use a formula to change the z variable as z → 1 − z,

$$\begin{split} {2}F{1}(a,b;c;z) &= \frac{\Gamma(c)\Gamma(c-a-b)}{\Gamma(c-a)\Gamma(c-b)} {}{2}F{1}(a,b;a+b+1-c;1-z) \ &+ \frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)} (1-z)^{c-a-b} {}{2}F{1}(c-a,c-b;1+c-a-b;1-z) \end{split} \tag{C.7}$$

when the alternative expression is found more convenient to understand the structure of singularities.

D Perturbative series

We present the two-loop expressions in the small- and large-m² t limits. S for small-m² t is given by

$$\overline{S} = -\frac{N{\rm c}}{8\pi^2} \frac{\overline{m}}{t} \sum{n,k=0}^{\infty} C{n,k}^{S,\ll 1} (\overline{m}^2, t, \mu^2) (\overline{m}^2 t)^k \left( \frac{\alpha{\rm s}(\mu^2)}{\pi} \right)^n \tag{D.1}$$

where

$$C_{0,0}^{S,\ll 1} = 1,$$

$$C_{0,1}^{S,\ll 1} = 2(L_t - L_m),$$

$$C_{0,2}^{S,\ll 1} = 4(L_t - L_m - 1),$$

$$C_{0,3}^{S,\ll 1} = 4L_t - 4L_m - 6,$$

$$C_{0,4}^{S,\ll 1} = \frac{4}{9} (6L_t - 6L_m - 11),$$

$$C_{0,5}^{S,\ll 1} = \frac{1}{9} (12L_t - 12L_m - 25),$$ (D.2)

and

$$\begin{split} &C{1,0}^{S,\ll 1} = C{\mathrm{F}}(1+2\log 2), \ &C{1,1}^{S,\ll 1} = C{\mathrm{F}}\left(3+10\log 2-9\log 3-6\operatorname{Li}{2}\left(\frac{1}{4}\right)-4L{m}-3L{m}^{2}+7L{t}+3L{m}L{t}\right), \ &C{1,2}^{S,\ll 1} = C{\mathrm{F}}\left(-11.5053-15.25L{m}-9L{m}^{2}+9.25L{t}+6L{m}L{t}+3L{t}^{2}\right), \ &C{1,3}^{S,\ll 1} = C{\mathrm{F}}\left(-12.6989-20.4375L{m}-11L{m}^{2}-0.5625L{t}+4L{m}L{t}+7L{t}^{2}\right), \ &C{1,4}^{S,\ll 1} = C{\mathrm{F}}\left(-4.2438-13.3656L{m}-7.5L{m}^{2}-11.9677L{t}-L{m}L{t}+8.5L{t}^{2}\right), \ &C{1,5}^{S,\ll 1} = C{\mathrm{F}}\left(4.2571-2.62156L{m}-2.75L{m}^{2}-16.2118L{t}-4.5L{m}L{t}+7.25L{t}^{2}\right). \end{split}$$ (D.3)

The log independent coefficients are valid to the numbers of significant digits shown. On the other hand, the coefficients of logarithmic terms are calculated to higher precision and numerical values such as 4.000 · · · are replaced by the integer numbers, and also, for instance, 7.5000 · · · of the L 2 ^m coefficient in C S,≪1 1,4 is shown as 7.5.

For large-m² t, the expanion reads

$$\overline{S} = -\frac{N{\rm c}}{16\pi^2} \frac{1}{\overline{m}t^2} \sum{n,k=0}^{\infty} C{n,k}^{S,\gg 1} (\overline{m}^2, t, \mu^2) (\overline{m}^2 t)^{-k} \left(\frac{\alpha{\rm s}(\mu^2)}{\pi}\right)^n \tag{D.4}$$

where

$$C_{0,0}^{S,\gg 1} = 1,$$

$$C_{0,1}^{S,\gg 1} = -1,$$

$$C_{0,2}^{S,\gg 1} = \frac{3}{2},$$

$$C_{0,3}^{S,\gg 1} = -3,$$

$$C_{0,4}^{S,\gg 1} = \frac{15}{2},$$

$$C_{0,5}^{S,\gg 1} = -\frac{45}{2},$$ (D.5)

and

$$\begin{split} C{1,0}^{S,\gg 1} &= -\frac{C{\mathrm{F}}}{4} \bigg( 2 + 3L_m + 3Lt + 6\log 2 - 9\log 3 \bigg), \ C{1,1}^{S,\gg 1} &= \frac{C_{\mathrm{F}}}{16} \bigg( 22 + 27L_m + 21Lt + 42\log 2 - 51\log 3 \bigg), \ C{1,2}^{S,\gg 1} &= -\frac{3C_{\mathrm{F}}}{64} \bigg( 68 + 87L_m + 57Lt + 114\log 2 - 127\log 3 \bigg), \ C{1,3}^{S,\gg 1} &= \frac{C_{\mathrm{F}}}{768} \bigg( 6634 + 8775L_m + 5049Lt + 10098\log 2 - 10719\log 3 \bigg), \ C{1,4}^{S,\gg 1} &= -\frac{C_{\mathrm{F}}}{9216} \bigg( 253556 + 341415L_m + 176985Lt + 353970\log 2 - 364095\log 3 \bigg), \ C{1,5}^{S,\gg 1} &= \frac{5C_{\mathrm{F}}}{36864} \bigg( 746314 + 1012095L_m + 480897L_t + 961794\log 2 - 967383\log 3 \bigg). \end{split}$$

R for small-m² t is given by

$$\overline{R} = -\frac{N{\rm c}}{8\pi^2} \frac{1}{t^2} \sum{n,k=0}^{\infty} C{n,k}^{R,\ll 1} (\overline{m}^2, t, \mu^2) (\overline{m}^2 t)^k \left(\frac{\alpha{\rm s}(\mu^2)}{\pi}\right)^n \tag{D.6}$$

where

$$C_{0,0}^{R,\ll 1} = 1,$$

$$C_{0,1}^{R,\ll 1} = -2,$$

$$C_{0,2}^{R,\ll 1} = 4(L_m - L_t),$$

$$C_{0,3}^{R,\ll 1} = 8(1 + L_m - L_t),$$

$$C_{0,4}^{R,\ll 1} = 12 + 8L_m - 8L_t,$$

$$C_{0,5}^{R,\ll 1} = \frac{8}{9} \left(11 + 6L_m - 6L_t\right).$$ (D.7)

and

$$\begin{split} C{1,0}^{R,\ll 1} &= \frac{C{\mathrm{F}}}{4} \left( 3\log 3 + 4\log 2 - 3L{t} \right), \ C{1,1}^{R,\ll 1} &= -\frac{C{\mathrm{F}}}{2} \left( 7 + 4\log 2 + 3L{t} \right), \ C{1,2}^{R,\ll 1} &= -\frac{C{\mathrm{F}}}{4} \left( 25 + 120\log 2 - 108\log 3 - 36\operatorname{Li}{2} \left( \frac{1}{4} \right) - 42L{m} - 24L{m}^{2} + 66L{t} + 12L{t}L{m} + 12L{t}^{2} \right), \ C{1,3}^{R,\ll 1} &= C{\mathrm{F}} (26.9641 + 29.8333L{m} + 14L{m}^{2} - 11.8333L{t} + 2L{m}L{t} - 16L{t}^{2}), \ C{1,4}^{R,\ll 1} &= C{\mathrm{F}} (22.0263 + 22.1667L{m} + 9.5L{m}^{2} + 28.8333L{t} + 23L{m}L{t} - 32.5L{t}^{2}), \ C{1,5}^{R,\ll 1} &= C{\mathrm{F}} (-11.3406 - 17.6604L{m} - 6L{m}^{2} + 75.6604L{t} + 48L{m}L{t} - 42L_{t}^{2}). \end{split} \tag{D.8}$$

For large-m² t, it is given by

$$\overline{R} = -\frac{N{\rm c}}{8\pi^2} \frac{1}{\overline{m}^2 t^3} \sum{n,k=0}^{\infty} C{n,k}^{R,\gg 1} (\overline{m}^2, t, \mu^2) (\overline{m}^2 t)^{-k} \left(\frac{\alpha{\rm s}(\mu^2)}{\pi}\right)^n \tag{D.9}$$

where

$$\begin{split} C{0,0}^{R,\gg 1} &= 1, \ C{0,1}^{R,\gg 1} &= -\frac{3}{2}, \ C{0,2}^{R,\gg 1} &= 3, \ C{0,3}^{R,\gg 1} &= -\frac{15}{2}, \ C_{0,4}^{R,\gg 1} &= \frac{45}{2}, \end{split} \tag{D.10}$$

and

$$\begin{split} C{1,0}^{R,\gg 1} &= -\frac{C{\mathrm{F}}}{16} \bigg( 4 + 18L_m + 18Lt + 36\log 2 - 45\log 3 \bigg), \ C{1,1}^{R,\gg 1} &= \frac{C_{\mathrm{F}}}{64} \bigg( 94 + 204L_m + 156Lt + 312\log 2 - 351\log 3 \bigg), \ C{1,2}^{R,\gg 1} &= -\frac{C_{\mathrm{F}}}{2304} \bigg( 11686 + 22086L_m + 14202Lt + 28404\log 2 - 30267\log 3 \bigg), \ C{1,3}^{R,\gg 1} &= \frac{C_{\mathrm{F}}}{9216} \bigg( 168494 + 298080L_m + 168480Lt + 336960\log 2 - 347085\log 3 \bigg), \ C{1,4}^{R,\gg 1} &= -\frac{C_{\mathrm{F}}}{12288} \bigg( 891478 + 1510650L_m + 770310L_t + 1540620\log 2 - 1549935\log 3 \bigg). \end{split}$$

As mentioned, the analytic expressions of C S,≪1 1,0 , C S,≪1 1,1 , C R,≪1 1,0 , C R,≪1 1,1 , and C R,≪1 1,2 are obtained by the method explained in sec. 2.6.

We provide these expansions as ancillary files that can be loaded in Mathematica.

Based on the numerical computation explained in sec. 2.5 we provide a grid for dimensionless S¯ and R¯. The format is {m² (µ)t, ^t m(µ) S¯} or {m² (µ)t, t2R¯}, where the perturbative series contain αs(µ 2 ), log(tµ² ), and the color factors as parameters. Also the uncertainties of the perturbative coefficients are labeled by δs[...] or δr[...]. Setting them to 0 (±1) corresponds to the central values (estimated upper/lower bounds) of the perturbative coefficients.

References

• [1] R. Narayanan and H. Neuberger, Infinite N phase transitions in continuum Wilson loop operators, JHEP 03 (2006) 064 [hep-th/0601210].
• [2] M. Lüscher, Properties and uses of the Wilson flow in lattice QCD, JHEP 08 (2010) 071 [1006.4518].
• [3] M. Lüscher, Chiral symmetry and the Yang–Mills gradient flow, JHEP 04 (2013) 123 [1302.5246].
• [4] H. Makino and H. Suzuki, Lattice energy–momentum tensor from the Yang–Mills gradient flow—inclusion of fermion fields, PTEP 2014 (2014) 063B02 [1403.4772].
• [5] Flavour Lattice Averaging Group (FLAG) collaboration, FLAG Review 2024, 2411.04268.
• [6] G. Martinelli, C. Pittori, C.T. Sachrajda, M. Testa and A. Vladikas, A general method for non-perturbative renormalization of lattice operators, Nucl. Phys. B 445 (1995) 81 [hep-lat/9411010].
• [7] Y. Aoki et al., Nonperturbative renormalization of quark bilinear operators and B^K using domain wall fermions, Phys. Rev. D 78 (2008) 054510 [0712.1061].
• [8] C. Sturm, Y. Aoki, N.H. Christ, T. Izubuchi, C.T.C. Sachrajda and A. Soni, Renormalization of quark bilinear operators in a momentum-subtraction scheme with a nonexceptional subtraction point, Phys. Rev. D 80 (2009) 014501 [0901.2599].
• [9] HPQCD collaboration, High-precision charm-quark mass and QCD coupling from current-current correlators in lattice and continuum QCD, Phys. Rev. D 78 (2008) 054513 [0805.2999].
• [10] TUMQCD collaboration, Relations between heavy-light meson and quark masses, Phys. Rev. D 97 (2018) 034503 [1712.04983].
• [11] J. Artz, R.V. Harlander, F. Lange, T. Neumann and M. Prausa, Results and techniques for higher order calculations within the gradient-flow formalism, JHEP 06 (2019) 121 [1905.00882].
• [12] F. Lange, Applications of the perturbative gradient flow at higher orders in quantum chromodynamics, Ph.D. thesis, RWTH Aachen University, 2021. 10.18154/RWTH-2021-07652.
• [13] H. Takaura, R.V. Harlander and F. Lange, Determining the Quark Mass with the Gradient Flow, PoS LATTICE2024 (2025) 288 [2411.13782].
• [14] R.V. Harlander, T. Nellopoulos, A. Olsson and M. Wesle, ftint: Calculating gradient-flow integrals with pySecDec, Comput. Phys. Commun. 306 (2025) 109384 [2407.16529].

• [15] M. Beneke and V.A. Smirnov, Asymptotic expansion of Feynman integrals near threshold, Nucl. Phys. B 522 (1998) 321 [hep-ph/9711391].

• [16] R. Harlander, The gradient flow at higher orders in perturbation theory, PoS LATTICE2021 (2022) 489 [2111.14376].

• [17] M. Beneke, P. Hager and A.F. Sanfilippo, Cosmological correlators in massless ϕ 4 -theory and the method of regions, JHEP 04 (2024) 006 [2312.06766].

• [18] M. Lüscher and P. Weisz, Perturbative analysis of the gradient flow in non-abelian gauge theories, JHEP 02 (2011) 051 [1101.0963].

• [19] F.V. Tkachov, A theorem on analytical calculability of 4-loop renormalization group functions, Phys. Lett. B 100 (1981) 65.

• [20] K.G. Chetyrkin and F.V. Tkachov, Integration by parts: The algorithm to calculate β-functions in 4 loops, Nucl. Phys. B 192 (1981) 159.

• [21] M. Neubert, Scale setting in QCD and the momentum flow in Feynman diagrams, Phys. Rev. D 51 (1995) 5924 [hep-ph/9412265].

• [22] R. Kitano and H. Takaura, Quantum electrodynamics on the lattice and numerical perturbative computation of g−2, PTEP 2023 (2023) 103B02 [2210.05569].

• [23] V.A. Smirnov, Analytical result for dimensionally regularized massless on-shell double box, Phys. Lett. B 460 (1999) 397 [hep-ph/9905323].

• [24] J.B. Tausk, Non-planar massless two-loop Feynman diagrams with four on-shell legs, Phys. Lett. B 469 (1999) 225 [hep-ph/9909506].

• [25] G. Mishima, High-energy expansion of two-loop massive four-point diagrams, JHEP 02 (2019) 080 [1812.04373].

• [26] H. Zhang, Massive two-loop four-point Feynman integrals at high energies with AsyInt, JHEP 09 (2024) 069 [2407.12107].

• [27] T. Huber and D. Maître, HypExp, a Mathematica package for expanding hypergeometric functions around integer-valued parameters, Comput. Phys. Commun. 175 (2006) 122 [hep-ph/0507094].

• [28] T. Huber and D. Maître, HypExp 2, Expanding hypergeometric functions about half-integer parameters, Comput. Phys. Commun. 178 (2008) 755 [0708.2443].

• [29] Wolfram Research, Inc., "Mathematica, Version 13.2."

• [30] T. Binoth and G. Heinrich, An automatized algorithm to compute infrared divergent multi-loop integrals, Nucl. Phys. B 585 (2000) 741 [hep-ph/0004013].

• [31] T. Binoth and G. Heinrich, Numerical evaluation of multi-loop integrals by sector decomposition, Nucl. Phys. B 680 (2004) 375 [hep-ph/0305234].

• [32] S. Borowka, G. Heinrich, S.P. Jones, M. Kerner, J. Schlenk and T. Zirke, SecDec-3.0: Numerical evaluation of multi-scale integrals beyond one loop, Comput. Phys. Commun. 196 (2015) 470 [1502.06595].

• [33] S. Borowka, G. Heinrich, S. Jahn, S.P. Jones, M. Kerner, J. Schlenk et al., pySecDec: A toolbox for the numerical evaluation of multi-scale integrals, Comput. Phys. Commun. 222 (2018) 313 [1703.09692].

• [34] S. Borowka, G. Heinrich, S. Jahn, S.P. Jones, M. Kerner and J. Schlenk, A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec, Comput. Phys. Commun. 240 (2019) 120 [1811.11720].

• [35] G. Heinrich, S.P. Jones, M. Kerner, V. Magerya, A. Olsson and J. Schlenk, Numerical scattering amplitudes with pySecDec, Comput. Phys. Commun. 295 (2024) 108956 [2305.19768].

• [36] M. Lüscher, Future applications of the Yang-Mills gradient flow in lattice QCD, PoS LATTICE2013 (2014) 016 [1308.5598].

• [37] R.V. Harlander, F. Lange and T. Neumann, Hadronic vacuum polarization using gradient flow, JHEP 08 (2020) 109 [2007.01057].

• [38] J. Borgulat, R.V. Harlander, J.T. Kohnen and F. Lange, Short-flow-time expansion of quark bilinears through next-to-next-to-leading order QCD, JHEP 05 (2024) 179 [2311.16799].

• [39] K. Hieda and H. Suzuki, Small flow-time representation of fermion bilinear operators, Mod. Phys. Lett. A 31 (2016) 1650214 [1606.04193].

• [40] E. Mereghetti, C.J. Monahan, M.D. Rizik, A. Shindler and P. Stoffer, One-loop matching for quark dipole operators in a gradient-flow scheme, JHEP 04 (2022) 050 [2111.11449].

• [41] J. Borgulat, R. Harlander, M.D. Rizik and A. Shindler, Two-loop matching of the chromo-magnetic dipole operator with the gradient flow, PoS LATTICE2022 (2023) 313 [2212.09824].

• [42] D.J. Broadhurst, Chiral symmetry breaking and perturbative QCD, Phys. Lett. B 101 (1981) 423.

• [43] V.P. Spiridonov and K.G. Chetyrkin, Nonleading mass corrections and renormalization of the operators mψψ¯ and G² µν, Sov. J. Nucl. Phys. 47 (1988) 522.

• [44] R. Harlander, Quarkmasseneffekte in der Quantenchromodynamik und asymptotische Entwicklung von Feynman-Integralen, Ph.D. thesis, Universität (TH) Karlsruhe, 1998.

• [45] R.V. Harlander, Y. Kluth and F. Lange, The two-loop energy–momentum tensor within the gradient-flow formalism, Eur. Phys. J. C78 (2018) 944 [1808.09837].

• [46] E. Braaten, S. Narison and A. Pich, QCD analysis of the tau hadronic width, Nucl. Phys. B 373 (1992) 581.

• [47] K.G. Chetyrkin and J.H. Kühn, Quartic mass corrections to Rhad, Nucl. Phys. B 432 (1994) 337 [hep-ph/9406299].

• [48] Particle Data Group collaboration, Review of Particle Physics, Phys. Rev. D 110 (2024) 030001.

• [49] P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, Five-Loop Running of the QCD Coupling Constant, Phys. Rev. Lett. 118 (2017) 082002 [1606.08659].

• [50] F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, The five-loop beta function of Yang-Mills theory with fermions, JHEP 02 (2017) 090 [1701.01404].

• [51] T. Luthe, A. Maier, P. Marquard and Y. Schröder, The five-loop beta function for a general gauge group and anomalous dimensions beyond Feynman gauge, JHEP 10 (2017) 166 [1709.07718].

• [52] K.G. Chetyrkin, B.A. Kniehl and M. Steinhauser, Strong Coupling Constant with Flavor Thresholds at Four Loops in the Modified Minimal-Subtraction Scheme, Phys. Rev. Lett. 79 (1997) 2184 [hep-ph/9706430].

• [53] J.A.M. Vermaseren, S.A. Larin and T. van Ritbergen, The 4-loop quark mass anomalous dimension and the invariant quark mass, Phys. Lett. B 405 (1997) 327 [hep-ph/9703284].

• [54] WHOT-QCD collaboration, Exploring N^f = 2+1 QCD thermodynamics from the gradient flow, Phys. Rev. D 96 (2017) 014509 [1609.01417].

• [55] WHOT-QCD collaboration, N^f = 2+1 QCD thermodynamics with gradient flow using two-loop matching coefficients, Phys. Rev. D 102 (2020) 014510 [2005.00251].