$\Lambda$ parameter of the SU(3) Yang-Mills theory from the continuous $\beta$ function

TL;DR

This paper nonperturbatively determines the SU(3) Yang-Mills theory's $eta$ function using the gradient flow scheme, enabling a precise calculation of the $ mar{MS}$ Lambda parameter across various coupling regimes.

Contribution

It provides the first comprehensive nonperturbative calculation of the $eta$ function for SU(3) Yang-Mills theory in the gradient flow scheme, covering a wide coupling range.

Findings

Calculated $eta$ function in the gradient flow scheme for SU(3) Yang-Mills.

Determined $ mar{MS}$ Lambda parameter as 0.622(10).

Results agree with recent direct determinations.

Abstract

Nonperturbative determinations of the renormalization group $\beta$ function are essential to connect lattice results to perturbative predictions of strongly coupled gauge theories and to determine the $\Lambda$ parameter or the strong coupling constant. The continuous $\beta$ function is very well suited for this task because it is applicable both in the weakly coupled deconfined regime as well as the strongly coupled confined regime. Here we report on our results for the $\beta$ function of the pure gauge SU(3) Yang-Mills theory in the gradient flow scheme. Our calculations cover the renormalized coupling range $g^2_{\textrm{GF}} \sim 1.2 - 27$, allowing for a direct determination of $\sqrt{8t_0} \Lambda_{\overline{\textrm{MS}}}$ in this system. Our prediction, $\sqrt{8 t_0} \Lambda_{\overline{\textrm{MS}}}=0.622(10)$, is in good agreement with recent direct determinations of this quantity.

Tables

Table 1. Table 1: The number of gauge field configurations analyzed and HMC acceptance rates for each bare coupling β b subscript 𝛽 𝑏 \beta_{b} and volume ( L / a ) 4 superscript 𝐿 𝑎 4 (L/a)^{4} . Configurations are separated by 20 MDTUs. All ensembles with β b > 5.5 subscript 𝛽 𝑏 5.5 \beta_{b}>5.5 are in the deconfined phase, while most with β b ≤ 5.50 subscript 𝛽 𝑏 5.50 \beta_{b}\leq 5.50 either transition between confined and deconfined or are confined. The values for β b = 5.30 subscript 𝛽 𝑏 5.30 \beta_{b}=5.30 are set in italics because these ensembles are not part of the main analysis and are only used to estimate systematic effects.

	$L/a$
	20		24		28		32		48
${\beta }_b$	Acceptance	No.	Acceptance	No.	Acceptance	No.	Acceptance	No.	Acceptance	No.
4.30	87.9%	451	86.6%	467	85.3%	297	80.7%	165	$\cdots$	$\cdots$
4.35	86.5%	451	84.8%	458	82.0%	277	78.8%	171	$\cdots$	$\cdots$
4.40	86.6%	451	80.8%	460	83.7%	272	82.3%	167	$\cdots$	$\cdots$
4.50	85.1%	451	84.2%	501	86.2%	1391	81.0%	250	$\cdots$	$\cdots$
4.60	86.1%	451	85.2%	490	83.3%	1040	84.9%	202	$\cdots$	$\cdots$
4.70	84.2%	451	84.1%	490	80.5%	681	82.1%	201	$\cdots$	$\cdots$
4.80	86.5%	451	88.0%	469	80.5%	681	78.9%	140	$\cdots$	$\cdots$
4.90	85.0%	451	85.3%	491	82.7%	701	83.4%	163	$\cdots$	$\cdots$
5.00	82.6%	451	85.5%	456	81.0%	772	77.3%	211	80.8%	124
5.30	84.4%	451	88.3%	534	82.9%	911	78.4%	656	81.8%	139
5.50	83.6%	451	87.6%	456	81.8%	701	77.8%	608	78.2%	149
6.00	84.4%	451	84.6%	476	84.6%	661	79.2%	472	76.8%	227
6.50	81.1%	451	80.7%	486	82.8%	661	85.0%	563	77.4%	233
7.00	81.3%	451	79.2%	461	81.7%	701	84.6%	527	74.7%	241
7.50	82.6%	451	81.3%	466	80.5%	661	83.7%	489	73.6%	224
8.00	81.3%	451	78.3%	456	76.1%	701	85.0%	487	73.3%	211
8.50	78.8%	451	77.4%	461	79.5%	661	81.6%	462	74.6%	211
9.00	78.2%	451	76.8%	581	78.0%	524	81.6%	531	71.7%	208
9.50	77.4%	621	77.5%	481	77.7%	547	81.7%	541	69.2%	208

Equations

$$\beta\big{(}g^{2}\big{)}=\mu^{2}\frac{\mathrm{d}g^{2}(\mu)}{\mathrm{d}\mu^{2}},$$

$$\frac{\Lambda^{2}}{\mu^{2}}=\big{(}b_{0}g^{2}(\mu)\big{)}^{-\frac{b_{1}}{b_{0}^{2}}}\exp\bigg{(}-\frac{1}{b_{0}g(\mu)^{2}}\bigg{)}\\ \times\exp\bigg{[}-\int_{0}^{g^{2}(\mu)}\!\!\!\!\mathrm{d}x\;\bigg{(}\frac{1}{\beta(x)}+\frac{1}{b_{0}x^{2}}-\frac{b_{1}}{b_{0}^{2}x}\bigg{)}\bigg{]},$$

$$\mu \propto 1/\sqrt{8t}.$$

$$g_{\text{GF}}^2(t)\equiv \mathcal{N}⟨t^{2}E(t)⟩,$$

$$g^{2}_{\textrm{GF}}\big{(}t;L,g_{0}^{2}\big{)}\equiv\frac{\mathcal{N}}{C(t,L)}\ \left\langle t^{2}E(t)\right\rangle,$$

$$\beta_{\textrm{GF}}\big{(}t;L,g_{0}^{2}\big{)}=-t\frac{\mathrm{d}}{\mathrm{d}t}g^{2}_{\textrm{GF}}\big{(}t;L,g_{0}^{2}\big{)}$$

$$\displaystyle\mathcal{I}_{N}\big{(}g_{\mathrm{GF}}^{2}\big{)}\equiv\frac{-p_{0}g_{\mathrm{GF}}^{4}\Big{(}1+\sum_{i=1}^{N}p_{i}g_{\mathrm{GF}}^{2i}\Big{)}}{1+\sum_{j=1}^{N+1}q_{j}g_{\mathrm{GF}}^{2j}},$$

$$\beta_{4}(g^{2}_{\textrm{GF}})\equiv-g^{4}_{\mathrm{GF}}\Big{(}b_{0}+b_{1}g^{2}_{\mathrm{GF}}+b_{2}g^{4}_{\mathrm{GF}}+b_{p}g^{6}_{\mathrm{GF}}\Big{)},$$

$${\int }_{g_i^2}^{g_f^2}\mathrm{d}x{\beta }^{-1}(x),$$

$$\sqrt{8t_0}{\Lambda }_{\mathrm{GF}}=1.164(19),$$

$$\sqrt{8t_0}{\Lambda }_{\overline{\text{MS}}}=0.622(10).$$

$$\frac{\beta_{\textrm{GF}}\left(\phi(g_{\mathcal{S}}^{2})\right)}{\phi^{\prime}\big{(}g_{\mathcal{S}}^{2}\big{)}}=\beta_{S}\big{(}g_{\mathcal{S}}^{2}\big{)},$$

$$\varphi (x)\approx x+x^2\sum _{n=0}^{N_p-1}{c}_i{x}^i,$$

Full text

††preprint: FERMILAB-PUB-23-096-V   SI-HEP-2023-04

I Introduction

The renormalization group (RG) $\beta$ function is defined as the logarithmic derivative of the renormalized running coupling $g^{2}$ with respect to an energy scale $\mu$

[TABLE]

and describes the scale dependence of the renormalized coupling of a 4-dimensional gauge-fermion system. Precise determination of the $\beta$ function is essential to understand the nonperturbative running of the gauge coupling and to predict, e.g., the $\Lambda$ parameter in quantum chromodynamics (QCD)

[TABLE]

where $g^{2}(\mu)$ is the renormalized coupling at the energy scale $\mu$ and the constants $b_{0}$ and $b_{1}$ denote the universal 1- and 2-loop coefficients of the RG $\beta$ function. Beyond 2-loop the $\beta$ function is scheme dependent, but it is straightforward to connect the $\Lambda$ parameter obtained using a different scheme to the phenomenologically favored $\overline{\mathrm{MS}}$ scheme using the 1-loop perturbative relation between the corresponding running couplings.

In the chiral limit of QCD, the $\Lambda$ parameter sets the physical scale. Its accurate determination allows the prediction of the strong coupling constant $\alpha_{s}(\mu)$ at $\mu\geq M_{Z}$, the mass of the $Z$ boson Aoki et al. (2022); Bruno et al. (2017a). $\Lambda$ parameters of systems with different flavor numbers can be connected by nonperturbative decoupling relations. In particular, it is possible to predict $\alpha_{s}(M_{Z})$ in $N_{f}=3$ QCD from the pure Yang-Mills $\Lambda$ parameter Dalla Brida et al. (2022). A precise calculation of the $\Lambda$ parameter in the pure Yang-Mills system is therefore desirable.

The RG transformation of an asymptotically free system has an ultraviolet (UV) fixed point (FP) on the $g^{2}_{0}=0$ critical surface, where $g_{0}^{2}$ denotes the (marginally) relevant bare coupling. The UVFP and the renormalized trajectory (RT) emerging from the UVFP describe continuum physics. The renormalized gauge coupling in Eq. (1) “measures” the flow along the RT. The gradient flow (GF) transformation Narayanan and Neuberger (2006); Lüscher (2010, 2010) is a particularly promising choice to define a renormalization scheme, determine the running coupling $g^{2}_{\textrm{GF}}(\mu)$, and calculate a nonperturbative $\beta$ function using lattice simulations. While on its own the GF transformation does not describe an RG transformation because it does not reduce the degrees of freedom of the system, it can be modified to define a complete RG scheme if a coarse-graining step is incorporated when defining expectation values Carosso et al. (2018); Carosso (2020). In this setup, the energy scale is set asymptotically by the GF flow time $t$

[TABLE]

For a generic RG flow, there is a one-to-one correspondence between running renormalized coupling and local observables that do not renormalize Makino et al. (2018). In the GF scheme the renormalized coupling can be defined in terms of the continuum flowed energy density $E(t)$ as

[TABLE]

since $E(t)$ has no anomalous dimension Lüscher (2010). The normalization factor $\mathcal{N}={128\pi^{2}}/(3N_{c}^{2}-3)$ in Eq. (4) is chosen such that the renormalized coupling $g^{2}_{\textrm{GF}}$ matches the renormalized coupling of the $\overline{\mathrm{MS}}$ scheme at tree-level.

The GF coupling offers a convenient low-energy hadronic scale $1/\sqrt{8t_{0}}$, with the flow time $t_{0}$ defined where the GF coupling takes a specified value, $g^{2}_{\textrm{GF}}(t_{0})\equiv\mathcal{N}\times 0.3$ Lüscher (2010). The $t_{0}$ scale has been used extensively in simulations to set the lattice scale and its continuum value is known in Yang-Mills systems Lüscher (2010); Dalla Brida and Ramos (2019) and in QCD Aoki et al. (2022) with 2+1 flavors Borsanyi et al. (2012); Blum et al. (2016); Bruno et al. (2017b) and 2+1+1 flavors Dowdall et al. (2013); Bazavov et al. (2016); Miller et al. (2021). Further, we can use $g^{2}_{\textrm{GF}}$ to define the RG $\beta$ function in the GF scheme $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})$. If $\beta_{\textrm{GF}}$ is determined up to $g^{2}_{\textrm{GF}}(t_{0})\approx 15.8$, we can obtain $\sqrt{8t_{0}}\Lambda_{\textrm{GF}}$ by integrating Eq. (2) up to $g^{2}_{\textrm{GF}}(t_{0})$.

The finite volume step-scaling function Lüscher et al. (1994); Fodor et al. (2012); Dalla Brida and Ramos (2019) is a commonly applied method to determine the $\beta$ function. However, it requires that the lattice size is the only dimensionful quantity of the system, preventing its application in the confining, chirally broken regime of QCD. Since the $t_{0}$ scale corresponds to a low energy hadronic scale in the confining regime, the step-scaling method cannot determine the $\beta$ function at strong enough gauge couplings to directly predict the $\Lambda$ parameter at the $t_{0}$ scale. An alternative approach is to determine the infinite volume continuous $\beta$ function as described in Refs. Hasenfratz and Witzel (2020, 2019); Fodor et al. (2018a). While this method requires separate infinite volume and continuum limit extrapolations, it is applicable in the confining regime.

In this paper, we report on our findings on the continuous $\beta$ function up to and even beyond the renormalized gauge coupling $g^{2}_{\textrm{GF}}(t_{0})$ in the pure gauge SU(3) Yang-Mills theory. Using the continuous $\beta$ function method we are able to determine the scale dependence of the running coupling in the confining regime. We do not need any other hadronic observable to determine the $\Lambda$ parameter at the $t_{0}$ scale. We compare our value of the $\Lambda$ parameter to other determinations from Refs. Dalla Brida and Ramos (2019); Wong et al. (2023) and values reported in the FLAG 2021 review Aoki et al. (2022). Preliminary results of our work were reported in Ref. Peterson et al. (2022). A similar calculation of the continuous $\beta$ function of SU(3) Yang-Mills theory reports preliminary results in Ref. Wong et al. (2023).

This paper is organized as follows. In Sec. II we summarize the details of our numerical simulations and explain our determination of the continuous $\beta$ function in Sec. III. In the subsequent section we discuss the determination of the $\Lambda$ parameter and close by discussing our results in Sec. V.

II Numerical Details

Our study is based on simulations performed using the tree-level improved Symanzik (Lüscher-Weisz) gauge action Lüscher and Weisz (1985a, b). We consider nineteen bare gauge couplings and generate configurations with periodic boundary conditions (BC) in all four directions using the hybrid Monte Carlo (HMC) update algorithm Duane et al. (1987) as implemented in GRID Boyle et al. (2015a, b). We set the trajectory length to $\tau=2$ molecular dynamics units (MDTU) and save configurations every 10 trajectories (20 MDTU). Subsequently, we use QLUA Pochinsky (2008); Pochinsky et al. (2008) to perform gradient flow measurements.

Table 1 lists the number of thermalized configurations analyzed for each bare coupling and volume as well as the HMC acceptance rates which all range between 70% and 90%. In the strong coupling regime111To better distinguish between the RG $\beta$ function and the bare gauge coupling, we refer to the latter as $\beta_{b}\equiv 6/g_{0}^{2}$. ($\beta_{b}\leq 4.9$) we have four volumes, ($(L/a)^{4}=20^{4},24^{4},28^{4},32^{4}$), while at weaker couplings ($\beta_{b}\geq 5.0$) we add a fifth, $48^{4}$ volume.222Although $L/a=20$ ensembles were generated for all $\beta_{b}$, our analysis uses $L/a=20$ only for $\beta_{b}\leq 5.00$.

The Polyakov loop expectation value in Fig. 1 shows that most ensembles for $\beta_{b}\geq 5.50$ are in the deconfined phase, while most $\beta_{b}\leq 5.0$ ensembles are confining.333Note that this apparent phase transition is only a finite volume effect. Our simulations are performed at zero temperature using symmetric $(L/a)^{4}$ volumes. The transition of the Polyakov loop indicates when the system size becomes comparable to the deconfinement length scale and transits from the small volume $\varepsilon$-regime to the large volume $p$-regime. In the limit of infinite volume, the pure Yang-Mills system is confining at all values of the bare gauge coupling. Since $\beta_{b}=5.30$ sits in the transition region, we discard it from our main analysis but use it to check for systematic effects later on. We study the autocorrelation of the renormalized coupling $g^{2}_{\textrm{GF}}$ using the $\Gamma$-method Wolff (2004) and typically find integrated autocorrelation times of less than four measurements (80 MDTU). However, near the transition from the deconfined to the confined regime, the integrated autocorrelation times increase up to fifteen measurements (300 MDTU).

While simulations in the deconfined regime have zero topological charge, nonzero topological charges are expected in the confined region. As we decrease $\beta_{b}$ for a fixed volume, we indeed observe that nonzero topological charges arise and that their fluctuations increase. We study the effect of nonzero topological charge on $g^{2}_{\textrm{GF}}$ by “filtering” configurations according to topological sectors and compare the corresponding values of $g^{2}_{\textrm{GF}}$. Our data suggest that for the fast-running pure gauge system, the $\beta$ function is sufficiently large that the impact of nonzero topological charges is not statistically resolved.444The situation is quite different for (near)-conformal simulations with many dynamical flavors, where the $\beta$ function runs slow and any topology can significantly alter the final result Hasenfratz and Witzel (2021).

III GF coupling and the RG $\beta$ function

III.1 Definition of the GF coupling and $\beta$ function

We define the finite volume gradient flow coupling on the lattice Lüscher (2010); Fodor et al. (2012) as

[TABLE]

where $g^{2}_{0}=6/\beta_{b}$ is the bare gauge coupling. The term $C(t,L)=1+\delta(t,L)$ corrects for the zero modes caused by the periodic BC of the gauge fields Fodor et al. (2012). The energy density $\left\langle E(t)\right\rangle$ can be approximated by local gauge observables. In general, $g^{2}_{\textrm{GF}}(t;L,g^{2}_{0})$ also depends on the gauge action, the gradient flow, and the operator chosen to approximate $\langle E(t)\rangle$. In this work we use Symanzik improved gauge action and consider two different gradient flow transformations, Wilson (W) and Zeuthen (Z) flow. We determine the Wilson plaquette (W), the clover (C), and the tree-level improved Symanzik operator (S) to estimate $\langle E(t)\rangle$. It is possible to calculate the tree-level lattice corrections to $g^{2}_{\textrm{GF}}(t;L,g_{0}^{2})$ for a given action-flow-operator combination and include those in $C(t,L)$ Fodor et al. (2014). In our analysis, we consider $C(t,L)$ with and without tree-level corrections, and we refer to the former as “tree-level normalization” (tln). We use a shorthand notation to distinguish the different flow-operator combinations, e.g. ZS refers to Zeuthen flow and Symanzik operator. When using tln in the definition of $g^{2}_{\mathrm{GF}}(t;L,g_{0}^{2})$ we prepend “n”, i.e. nZS refers to tln improved Zeuthen flow and Symanzik operator. In the continuum limit, the RG $\beta$ function is independent of the bare action, GF transformation or operator choice, and the comparison of the different combinations can serve to estimate systematic effects. Continuum extrapolations on the lattice are more stable when lattice artifacts are small, and we chose the nZS combination as our preferred analysis.

From the GF coupling in Eq. (5), we derive the GF $\beta$ function

[TABLE]

in finite volume by discretizing the flow-time derivative $\mathrm{d}/\mathrm{d}t$ with a 5-point stencil. The interval of the time derivative is set by the time step $\epsilon$ used to integrate the gradient flow equation for the gauge field. We choose $\epsilon=0.04$ and explicitly verified that this choice has no impact on our analysis by repeating the GF measurements using $\epsilon=0.01$ for selected ensembles.

Our numerical analysis starts by determining the renormalized couplings $g^{2}_{\textrm{GF}}(t;L,g_{0}^{2})$ at all flow times for our three operators. We determine $g^{2}_{\textrm{GF}}(t;L,g_{0}^{2})$ for both data sets obtained with Zeuthen and Wilson flow, respectively, and for our entire set of gauge field ensembles. Next, we numerically calculate the derivative as specified in Eq. (6). The uncertainties are propagated using standard correlated error propagation techniques implemented in the software packages gvar Lepage (2015) and lsqfit Lepage (2014).

The analysis of the RG $\beta$ function proceeds by first taking the infinite volume limit $a/L\rightarrow 0$ of both $g^{2}_{\textrm{GF}}(t;L,g_{0}^{2})$ and $\beta_{\textrm{GF}}(t;L,g_{0}^{2})$ independently at fixed $\beta_{b}$ and $t/a^{2}$. This is followed by an interpolation of $\beta_{\textrm{GF}}(t;g_{0}^{2})$ in $g^{2}_{\textrm{GF}}(t;g_{0}^{2})$ at fixed $t/a^{2}$. The last step of our analysis is to take the $a^{2}/t\to 0$ continuum limit at fixed $g^{2}_{\textrm{GF}}$. These steps are detailed in the rest of this section.

III.2 Infinite volume extrapolation

In a 4-dimensional gauge-fermion system the volume is a relevant parameter, and the RG equation in finite volume includes a term describing its effect on the running of the renormalized coupling. We prefer to avoid complications arising from this a priori unknown quantity and define the RG $\beta$ function in the infinite volume limit. This requires first extrapolating to infinite volume before considering the continuum limit.

Since the energy density $\left\langle E(t)\right\rangle$ is a dimension-4 operator, the finite volume corrections of $g^{2}_{\textrm{GF}}$ are expected to be $\mathcal{O}((a/L)^{4})$ at leading order. We independently extrapolate both the GF coupling $g^{2}_{\textrm{GF}}(t;L,g_{0}^{2})$ and GF $\beta$ function \beta_{\textrm{GF}}\big{(}t;L,g_{0}^{2}\big{)} linearly in $(a/L)^{4}\rightarrow 0$ for each fixed bare gauge coupling $\beta_{b}$ and lattice flow time $t/a^{2}$. This analysis strategy was first outlined in Ref. Peterson et al. (2022). Alternative methods are discussed, e.g., in Refs. Fodor et al. (2018b); Hasenfratz and Witzel (2019, 2020); Kuti et al. (2022).

In Fig. 2 we show typical infinite volume extrapolations for $g^{2}_{\textrm{GF}}(t;L,g_{0}^{2})$ and $\beta_{\textrm{GF}}(t;L,g_{0}^{2})$ at relatively weak coupling ($\beta_{b}=6.0$), at intermediate couplings ($\beta_{b}=4.9$ and 5.50), and in the strong coupling regime ($\beta_{b}=4.35$). Each panel shows the $(a/L)^{4}$ extrapolation at five different lattice flow time values that cover the range we use in the continuum limit extrapolation (cf. Sec. III.5). In the strong coupling regime, we observe very mild volume dependence. This is consistent with the expectation that below the confinement scale $\mu_{\text{conf}}$, i.e. at renormalized couplings stronger than $g^{2}_{\textrm{GF}}(\mu_{\text{conf}})$, the confinement scale provides an infrared cutoff and the volume dependence is suppressed. Thus we find it sufficient to use volumes with $L/a=\{20,\,24,\,28,\,32\}$ in this regime. At renormalized couplings that correspond to energy scales above the confinement scale, the volume dependence is more significant, as the only infrared cutoff is due to the finite volume. Here we drop the $L/a=20$ ensembles and add $48^{4}$ lattices in the infinite volume extrapolation. At $\beta_{b}=5.5$ and 6.0 volumes $L/a\leq 32$ appear deconfined, while $L/a=48$ is transitioning to deconfined at $\beta_{b}=6.0$ and is confined at $\beta_{b}=5.5$. Similarly the $L/a=32$ at $\beta=5.5$ ensemble is transitioning and exhibits very long autocorrelation times. Likely these long autocorrelations result in underestimated statistical errors causing the $2\sigma$ deviation. At $\beta_{b}=5.0$ the Polyakov loop expectation value shows that volumes $L/a\geq 28$ are confining, while $L/a=20$ and 24 are in the transition region. $L/a=24$ shows especially large autocorrelation times. Nevertheless, at $\beta_{b}=5.0$ we use all five volumes in the infinite volume extrapolation, though it is dominated by the larger volumes. We always perform the infinite volume extrapolation using a linear fit ansatz in $(a/L)^{4}$. For flow times $2.0\lesssim t/a^{2}\lesssim 4.0$ contributing to the continuum extrapolations in Sec. III.5, most fits have good $p$-values ($\gtrsim 10\%$). Notable exceptions are $\beta_{b}=5.0$ and 5.5, where large autocorrelations make it difficult to estimate the errors correctly. When adding $\beta_{b}=5.3$ to our analysis, the infinite volume extrapolations have vanishing $p$-values. The situation did not improve despite extending the Monte Carlo simulations and using larger thermalization cuts. We therefore conclude that $\beta_{b}=5.3$ suffers from sitting in the transition region and discard it from our main analysis.

III.3 Different flows, operators and tree-level improvement

In Figs. 3 and 4 we show the infinite volume extrapolated GF coupling $g^{2}_{\textrm{GF}}(t;g_{0}^{2})$ as a function of the gradient flow time $t/a^{2}$ with and without tln improvement at $\beta_{b}=6.0$ and 4.35, respectively. We consider both Zeuthen and Wilson flow and all three operators. In the continuum limit the different flows and operators must agree, so any difference between them at finite bare coupling points to cutoff effects. Figure 3 showcases the improvement due to tln at weak coupling. As the left panels show, there are significant differences between the three operators even at flow time $t/a^{2}=4.0$. Comparing the top and bottom panels on the left, large differences for $g^{2}_{\textrm{GF}}(t;g_{0}^{2})$ between Zeuthen (top) and Wilson (bottom) flow can be seen. Both the operator and flow dependence are largely removed by tln normalization, which is demonstrated in the panels on the right.

Figure 4 shows the same quantities at a strong bare coupling, $\beta_{b}=4.35$. The gauge coupling runs fast and the same flow time range covers nearly 20 times the range of $g^{2}_{\mathrm{GF}}$ at $\beta_{b}=4.35$ compared to $\beta_{b}=6.0$. While the different operators and flows without tln correction may appear to be closer at $\beta_{b}=4.35$ than at $\beta_{b}=6.0$, the absolute difference is much larger. Tree-level correction works similarly in the strong coupling, removing most flow and operator dependence. The somewhat larger spread after tln suggests that the perturbatively motivated improvement is not as effective in the strong coupling regime as at weak coupling. The effect of tln improvement at strong coupling was recently studied in the context of gradient flow scale setting for the case of 2+1 flavor gauge field configurations generated with Iwasaki gauge action and Shamir domain-wall fermions Schneider et al. (2023). That work also showed that tln removed most flow and operator dependence, but the improved predictions still had significant cutoff effects. However, here we study the pure-gauge system where no effect due to fermion loops can spoil the improvement and, furthermore, we use the perturbatively improved Symanzik gauge action in comparison to Iwasaki gauge action. In Sec. III.5 we will show that in our case tln improvement not only reduces the difference between different flows and operators but also reduces cutoff effects.

Based on the improvement from tln, we designate the tree-level improved Zeuthen flow with Symanzik operator combination (nZS) as our preferred analysis. In the remainder of this paper we therefore concentrate on analyzing this combination. Other tln improved flow-operator combinations are consistent within uncertainties and considered as part of our discussion on systematic uncertainties in III.7. Moreover, we repeat our analysis without tln improvement, which increases discretization errors. Within these more sizeable uncertainties, the unimproved combinations are consistent with our preferred nZS prediction.

III.4 Interpolation of $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})$ vs. $g^{2}_{\textrm{GF}}$

The last task before we can take the $a^{2}/t\to 0$ continuum limit is the identification of \{\beta_{\textrm{GF}}\big{(}t;g^{2}_{0}),\,t/a^{2}\} pairs at fixed g^{2}_{\textrm{GF}}\big{(}t;g^{2}_{0}). We achieve this by interpolating $\beta_{\textrm{GF}}(t;g^{2}_{0})$ in terms of $g^{2}_{\textrm{GF}}(t;g^{2}_{0})$ at fixed values of lattice flow times $t/a^{2}$.

Our interpolating form has to be able to describe the weak coupling regime where we expect perturbative behavior \beta(g_{\mathrm{GF}}^{2})\sim-b_{0}g_{\mathrm{GF}}^{4}+\mathcal{O}\big{(}g_{\mathrm{GF}}^{6}\big{)} as well as the strong coupling regime where it appears that $\beta_{\mathrm{GF}}(g_{\mathrm{GF}}^{2})\approx c_{0}+c_{1}g_{\mathrm{GF}}^{2}$. A polynomial is not appropriate to cover both regimes. Instead, we choose a form that is similar to the ratio of polynomials used in Padé approximations

[TABLE]

such that we enforce the leading order $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})\sim g_{\mathrm{GF}}^{4}$ behavior. Since cutoff effects can change even the asymptotic behavior of the $\beta$ function, this ansatz could constrain the lower limit on the flow time $t/a^{2}$ used in the analysis. However with $N=4$ we find that \mathcal{I}_{N}\big{(}g^{2}_{\textrm{GF}}\big{)} provides a good description with $p$-values between 17% – 32% for flow times $2.0\lesssim t/a^{2}\lesssim 4.0$.

The two panels of Fig. 5 illustrate the interpolation at various flow times in the range $t/a^{2}\in[2.0,4.0]$. The colored data points on both panels correspond to infinite volume limit $(g^{2}_{\textrm{GF}},\,\beta_{\textrm{GF}})$ pairs at each bare coupling and five selected flow time values. For reference we show the universal 1- and 2-loop and the GF 3-loop Harlander and Neumann (2016) perturbative lines in gray using dashed, dotted, and dash-dotted lines, respectively. The colored bands describe the interpolation according to Eq. (7) with $N=4$. The top panel of Fig. 5 shows $\beta_{\textrm{GF}}$ vs. $g^{2}_{\textrm{GF}}$ in the entire range of our data. The data points at different bare couplings and flow times form a smooth curve, indicating that cutoff effects are mild. The corresponding interpolating curves basically overlap, creating the purple-hued band. This plot also illustrates that the $\beta$ function is approximately linear in the strong coupling regime. We will discuss this feature further in Sec. III.8. In the insert of the top panel we zoom in on the weak coupling regime. To enhance the small $g^{2}_{\textrm{GF}}$ behavior, we plot $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})/g_{\mathrm{GF}}^{4}$ vs. $g^{2}_{\textrm{GF}}$ on the bottom panel. At large flow times and small couplings, $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})$ has to be consistent with the perturbative predictions. As the panel shows, both the data and the interpolating forms are close to the perturbative values in our selected flow time range, though we do not constraint the intercept of $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})/g^{4}_{\textrm{GF}}$.

III.5 Continuum limit extrapolation

To obtain the continuous $\beta$ function in the continuum limit, we need to perform one last step, the $a^{2}/t\to 0$ continuum extrapolation. We use a linear extrapolation in $a^{2}/t\to 0$ and the flow time $t/a^{2}$ has to be chosen large enough for the RG flow to reach the renormalized trajectory where a linear extrapolation describes the remnant cutoff effects (within statistical uncertainties). We also have to limit the maximum flow time because the infinite volume extrapolation is reliable only when the finite volume corrections follow the leading order $t^{2}/L^{4}$ scaling form. In practice, we vary $t_{\mathrm{min}}/a^{2}$ and $t_{\mathrm{max}}/a^{2}$ and monitor the quality of the linear extrapolation. In Fig. 6 we show examples of the continuum extrapolation from our weakest available coupling of $g^{2}_{\textrm{GF}}$ up to the $t_{0}$ scale of $g^{2}_{\textrm{GF}}(t_{0})=0.3\,\mathcal{N}\approx 15.8$. In all cases, we use $t/a^{2}\in[2.0,4.0]$ but show additional data points both at smaller and larger flow times to illustrate the linear behavior in $t/a^{2}$.

It is worth pointing out that the flow time is technically a continuous variable and the corresponding $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})$ values are highly correlated. We obtain the final continuum limit prediction by selecting a discrete set of $t/a^{2}$ values. For this set we perform an uncorrelated fit and estimate its statistical uncertainty by repeating this fit using the central values shifted by $\pm 1\sigma$. This avoids the complication of inverting a poorly conditioned correlation matrix and ensures we are not underestimating the statistical uncertainty.

III.6 The nonperturbative $\beta_{\textrm{GF}}$ function

Figure 7 shows our result for the nonperturbative $\beta_{\textrm{GF}}$ function for pure gauge SU(3) Yang-Mills in the gradient flow scheme with statistical errors only. In addition to the continuum limit prediction shown with a salmon-colored band, the plot also shows the infinite volume extrapolated nZS lattice data at flow times $t/a^{2}\in[2.0,4.0]$ where, for better visibility, we “thin” the data and only show every fifth data point, i.e. flow time values are separated by $\Delta t/a^{2}=0.2$. The nZS combination shows very little cutoff dependence and the raw lattice data sits on top of the continuum extrapolated value. Overall, our results span the coupling range from $g^{2}_{\textrm{GF}}\gtrsim 1.2$ up to $g^{2}_{\textrm{GF}}\lesssim 27$, well into the confining regime of the system.

III.7 Systematic uncertainties of the $\beta_{\textrm{GF}}$ function

In addition to the statistical uncertainties shown for the final result of our $\beta_{\textrm{GF}}$ function in Fig. 7, we check for systematic effects by considering variations of our preferred nZS analysis. We discuss the different variations below and show the outcome as relative changes w.r.t. our nZS analysis in Fig. 9.

• •

For the interpolation in $g^{2}_{\textrm{GF}}$ we use a ratio of polynomials similar to Padé approximation. We alter the functional form by changing the order of the polynomials in Eq. (7) from $N=4$ to $N=2$ or 6, but find this is only a subleading effect on our resulting $\beta_{\textrm{GF}}$ function. In Fig. 9 we show the (larger) effect by taking the absolute difference between interpolations based on $N=4$ and $N=2$.

• •

Our simulations at $\beta_{b}=5.30$ sit right at the deconfined/confined transition and we observe poor $p$-values when performing the infinite volume extrapolation at this bare gauge coupling. We therefore discarded $\beta_{b}=5.30$ from our main analysis, but we use it now to check for systematic effects by including it as an additional data point in a sensitive region of the $\beta_{\textrm{GF}}$ function. As expected, the impact of adding $\beta_{b}=5.30$ is largest around $g^{2}_{\textrm{GF}}\sim 5$.

• •

To validate our infinite volume extrapolation, we consider two variations:

We drop the smallest volume and perform a linear fit to our three largest volumes. 2. 2.

We note that in Fig. 2 our largest volume $L/a=32$ ($L/a=48$) for $\beta_{b}<5.00$ ($\beta_{b}\geq 5.00$) is very close to the extrapolated infinite volume value. Therefore we simply repeat our analysis using only the largest volumes.

Both variations result in comparable uncertainties, but just using the largest volumes has a slightly larger effect. Hence we use the latter to estimate finite volume effects. As can be seen in Fig. 9 this could be our dominant systematic uncertainty for strong couplings $g^{2}_{\textrm{GF}}\gtrsim 11$.

• •

We test the continuum limit extrapolation by varying the range of the flow times entering the linear fit in $a^{2}/t$ at fixed values of $g^{2}_{\textrm{GF}}$. Keeping $t_{\text{max}}/a^{2}=4.0$ fixed, we vary $t_{\text{min}}/a^{2}$ from 1.52 to 2.0. Similarly, we vary $t_{\text{max}}/a^{2}$ from 4.0 to 5.0 while $t_{\text{min}}/a^{2}=2.0$ is kept fixed. Comparing these to our preferred analysis we see at most a variation of $\mathcal{O}(0.3\%)$ in the central value of our $\beta_{\textrm{GF}}$ function. We show the maximum of these variations as “fit range” uncertainty in Fig. 9. Further, we repeat the fits using our preferred fit range but reduce the number of data points fitted by increasing the separation in flow time. As expected this variation results only in minuscule changes and can be neglected.

• •

In the continuum limit different flows and operators should predict the same renormalized $\beta_{\textrm{GF}}$ function. Qualitatively this is the case, as Fig. 8 demonstrates. There we show six different tln improved flow-operator combinations to determine $\beta_{\textrm{GF}}$ which all sit on top of each other forming nearly a single band. Looking at the relative changes in Fig. 9 we do, however, see deviations of $\mathcal{O}(0.8\%)$. Since for couplings in the range $7\lesssim g^{2}_{\textrm{GF}}\lesssim 11$ that effect is larger than other systematic effects, we conservatively include these variations when obtaining our systematic uncertainty.

• •

We further check for consistency by analyzing our data without using tln improvement. We find that removing the tln improvement increases the discretization errors noticeably, but within the larger uncertainties, the results are consistent with our preferred analysis.

Using the information compiled in Fig. 9, we obtain the total uncertainty of our $\beta_{\textrm{GF}}$ function by adding our statistical error and the largest systematic effect at each $g^{2}_{\textrm{GF}}$ value in quadrature. Our final result for the $\beta_{\textrm{GF}}$ function for the range of $g^{2}_{\textrm{GF}}$ relevant to determine $\sqrt{8t_{0}}\Lambda_{\textrm{GF}}$ (see below) is shown in Fig. 10.555Our final result for $\beta_{\textrm{GF}}$ is also provided as an ASCII file sup .

III.8 $\beta_{\textrm{GF}}$ in the strong coupling regime

As can be seen in particular in Fig. 8, the $\beta_{\textrm{GF}}$ function in the strong coupling regime ($g^{2}_{\textrm{GF}}\gtrsim 13$) is approximately linear. Studying the derivative $\mathrm{d}\beta_{\textrm{GF}}/\mathrm{d}g^{2}_{\textrm{GF}}$, we observe a plateau within errors in the range $20\lesssim g^{2}_{\textrm{GF}}\lesssim 27$. In that range we can approximate $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})\approx c_{0}+c_{1}g^{2}_{\textrm{GF}}$ with $c_{0}=5.8(3)$, $c_{1}=-1.32(1)$. This is very different from the perturbative prediction that would suggest a polynomial form with terms of $g_{\mathrm{GF}}^{4}$ and higher order. Clearly, nonperturbative effects are at play here. We may use the above observation to predict the flow time dependence of the flowed energy density \langle t^{2}E(t)\rangle\sim a_{0}+a_{1}\big{(}t/\tilde{t})^{-c_{1}}, where $a_{0}=-c_{0}/c_{1}\approx 4.4(2)$ and $\tilde{t}$ is an integration constant. We do not have a physical explanation for this behavior, though it is tempting to think that the form of $\langle E(t)\rangle$ as $t\rightarrow\infty$ is related to instantons, the only objects in the vacuum after UV fluctuations are removed by the flow. We mention that Refs. Nakamura and Schierholz (2023); Schierholz (2022); Ryttov and Sannino (2008); Chaichian and Frasca (2018) hypothesize that in confining systems the RG $\beta$ function is linear in the confining, strongly coupled regime. References Nakamura and Schierholz (2023); Schierholz (2022) attempt to solve the strong CP problem by connecting $\langle E(t)\rangle$ to the instanton density at large flow time. However, this work considers only one bare gauge coupling, does not take the continuum limit and uses much larger flow times than we do on similar volumes. While our results are consistent with a linear $\beta$ function in the strong coupling range, our slope $c_{1}=-1.32(1)$ is significantly smaller than -1 obtained in Refs. Nakamura and Schierholz (2023); Schierholz (2022). Moreover, we observe a nonzero intercept. It would be interesting to understand the nonperturbative origin of the linearity of the $\beta_{\textrm{GF}}$ function in the strong coupling confining region both in pure Yang-Mills systems and, if it persists, with dynamical fermions.

IV The $\Lambda$ parameter

IV.1 Matching the perturbative regime

In order to determine the $\Lambda$ parameter according to Eq. (2), we need to know $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})$ down to $g^{2}_{\textrm{GF}}=0$. In the weak coupling regime we expect to recover the perturbatively predicted $\beta$ function. In Fig. 11 we show $\beta_{\textrm{GF}}/g^{4}_{\textrm{GF}}$ to emphasize the weak coupling behavior of our nonperturbative result. The numerical data, shown by the salmon-colored band, is close to the 3-loop GF perturbative curve at our weakest gauge coupling $g^{2}_{\textrm{GF}}\approx$ 1.2, but it does not yet connect smoothly. To remedy this limitation we extend our numerically determined $\beta_{\textrm{GF}}$ function to the region below our weakest coupling data point $(g^{2}_{\textrm{GF}}\approx$ 1.2) using the parameterization

[TABLE]

where $b_{0}$, $b_{1}$, $b_{2}$ are the 1-, 2- and 3-loop GF perturbative coefficients from Ref. Harlander and Neumann (2016) and $b_{p}$ is free. We determine $b_{p}$ by integrating the inverse $\beta$ function

[TABLE]

using both our numerically determined $\beta_{\textrm{GF}}$ function and also using $\beta_{4}$ parametrized by $b_{p}$ for $\beta(x)$. We determine $b_{p}$ by equating these two values. The integration limits are set to cover the regime where we attempt to match $\beta_{\textrm{GF}}$ and $\beta_{4}$. In the example shown in Fig. 11, we choose to match in the range $g^{2}_{i}=1.4$, $g^{2}_{f}=1.8$. In Fig. 12 we demonstrate that varying $g_{i}^{2}$ and $g_{f}^{2}$ has negligible impact on our value of $\sqrt{8t_{0}}\Lambda_{\overline{\text{MS}}}$. To account for the combined statistical and systematic uncertainty of our nonperturbative $\beta_{\textrm{GF}}$ function, we repeat the matching shifting the central values by $\pm 1\sigma$. That way we obtain the purple bands resulting in the upper and lower end of the blue band connecting our nonperturbative result to $g^{2}_{\textrm{GF}}=0$.

IV.2 Calculating the $\Lambda$ parameter

The final step of this analysis is to calculate $\sqrt{8t_{0}}\Lambda_{\textrm{GF}}$ by integrating Eq. (2) using the combination of our nonperturbative $\beta_{\textrm{GF}}(g^{2}_{\textrm{GF}})$ for $g^{2}_{\textrm{GF}}\geq 1.4$ and the matched $\beta_{4}(g^{2}_{\textrm{GF}})$ function for $g^{2}_{\textrm{GF}}<1.4$. The upper integration limit is set according to $g^{2}_{\textrm{GF}}(t_{0})\approx 15.8$. Our prediction is

[TABLE]

where the error accounts for the uncertainty of our $\beta_{\textrm{GF}}$ function as well as the uncertainty encountered due to the matching procedure. As a last step, we convert our prediction to the $\overline{\text{MS}}$ scheme using the exact 1-loop relation and find

[TABLE]

IV.3 Comparison of $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}$ determinations

The Yang-Mills $\Lambda$ parameter has been studied previously using different approaches. Only the recent gradient flow studies in Refs. Dalla Brida and Ramos (2019); Wong et al. (2023) directly predict the combination $\sqrt{8t_{0}}\Lambda_{\mathrm{\overline{\textrm{MS}}}}$. In addition, we compare our value to determinations of $r_{0}\Lambda_{\overline{\textrm{MS}}}$ listed by the flavor lattice averaging group (FLAG) Aoki et al. (2022) to meet the quality criteria to enter the average. These determinations are obtained using Schrödinger functional step-scaling methods Capitani et al. (1999); Ishikawa et al. (2017), Wilson loops Gockeler et al. (2006); Kitazawa et al. (2016), or the short distance potential Brambilla et al. (2010). We use the values quoted by FLAG 2021 for $r_{0}\Lambda_{\overline{\textrm{MS}}}$ and convert them to $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}$ using $\sqrt{8t_{0}}/r_{0}=0.948(7)$ Lüscher (2010) (open symbols) or $\sqrt{8t_{0}}/r_{0}=0.9414(90)$ Dalla Brida and Ramos (2019) (filled symbols). Following the FLAG convention, we refer to the different results in Fig. 13 using either the name of the first author or, if applicable, the name of the collaboration and the two-digit year.

Given the spread in the values of $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}$, further scrutiny and understanding are needed before obtaining an average. We note, however, that the three most recent predictions are all mutually consistent. The high-precision result of Ref. Dalla Brida and Ramos (2019) was re-affirmed in Ref. Nada and Ramos (2021) using an alternative approach with better control over the continuum extrapolation. The estimate given in Ref. Schierholz (2022) is also consistent with these predictions. A possible source of difference to the older determinations is the conversion of $r_{0}$ to $\sqrt{8t_{0}}$.

IV.4 Nonperturbative matching of different schemes

In our analysis, a considerable systematic uncertainty arises from the weak coupling limit $g^{2}_{\textrm{GF}}\lesssim 1.4$. Our gradient flow setup is not efficient at weak coupling. It would be more economical to use data from existing calculations, e.g. the high precision Schrödinger functional data of Ref. Dalla Brida and Ramos (2019) in the $0<g^{2}_{\textrm{GF}}<1.4$ regime and match it nonperturbatively to our data.

Such a matching requires finding the relation between our $g^{2}_{\textrm{GF}}$ coupling and the coupling $g^{2}_{\mathcal{S}}$ of another scheme $\mathcal{S}$ expressed as g^{2}_{\textrm{GF}}=\phi\big{(}g^{2}_{\mathcal{S}}\big{)}. The relation of the corresponding $\beta$ functions can be obtained using the chain rule applied to the derivative of $g^{2}_{\textrm{GF}}$ with respect to $\mu^{2}$, which leads to the simple relation

[TABLE]

where \phi^{\prime}\big{(}g^{2}_{\mathcal{S}}\big{)}\equiv\mathrm{d}\phi\big{(}g^{2}_{\mathcal{S}}\big{)}/\mathrm{d}g^{2}_{\mathcal{S}}. Parametrizing $\phi$ as a polynomial

[TABLE]

turns Eq. (12) into a straightforward fitting problem with $N_{p}$ undetermined coefficients. The only constraint is to identify and use the renormalized coupling range in the fit where the two schemes overlap. Such a nonperturbative matching and combination of different schemes could lead to a significantly improved prediction. Although we do not explore this method in the present analysis, it is worth considering in the future.

V Discussion

In this paper, we present a nonperturbative determination of the renormalization group $\beta$ function for the pure gauge Yang-Mills action. Using the gradient flow based continuous RG $\beta$ function, we present results for a wide range of values of the renormalized running coupling. Our results span the range of the perturbative weak coupling region $g^{2}_{\textrm{GF}}\approx 1.2$ up to the strongly coupled regime at $g^{2}_{\textrm{GF}}\approx 27$. This showcases the advantage of the continuous RG $\beta$ function because the continuous infinite volume $\beta$ function can be extended without limitation to the confining region. We also demonstrate the effectiveness of tree-level improvement of the gradient flow even in the strong coupling regime.

We investigate various sources of systematical uncertainties. For most of the $g^{2}_{\textrm{GF}}$ range covered, our statistical uncertainties are around 0.6%. In the weak coupling region, statistical and systematic errors are of similar size. At intermediate to strong coupling, we observe an increase in the systematic errors to 1.5% due to enhanced finite volume effects.

While in the weak coupling our results are close to the perturbative values, we observe in the confining regime that the GF $\beta$ function depends approximately linearly on the running coupling, implying a scaling relation of the flowed energy density \langle t^{2}E(t)\rangle\sim a_{0}+a_{1}\big{(}t/\tilde{t}\big{)}^{-c_{1}} with exponent $c_{1}\approx-1.32(1)$. This observation could be related to the topological structure of the vacuum, a possibility that warrants further investigation.

In the weak coupling regime, we are able to match our numerical results to the 3-loop GF $\beta$ function by extending the perturbative expression with a single $g^{10}_{\mathrm{GF}}$ term. This matching allows us to predict the $\Lambda$ parameter in the GF scheme. Using the perturbatively determined relation of the GF coupling $g^{2}_{\textrm{GF}}$ and the $\overline{\textrm{MS}}$ coupling, we obtain $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}=0.622(10)$, where the error combines statistical and systematic uncertainties. This value is in good agreement with recent direct determinations of $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}$ Wong et al. (2023); Dalla Brida and Ramos (2019).

A significant source of the systematic uncertainties in determining $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}$ arises from the weak coupling regime. We outline a nonperturbative matching procedure to combine existing high-precision data in the weak coupling and our nonperturbative $\beta_{\textrm{GF}}$ function that extends into the confining regime even beyond $g^{2}_{\textrm{GF}}(t_{0})$. Such a combined determination could lead to a sizable reduction of the uncertainties on $\sqrt{8t_{0}}\Lambda_{\overline{\textrm{MS}}}$.

Acknowledgements.

We are very grateful to Peter Boyle, Guido Cossu, Anontin Portelli, and Azusa Yamaguchi, who develop the GRID software library providing the basis of this work and who assisted us in installing and running GRID on different architectures and computing centers. We benefited from many comments and discussions during the “Gradient Flow in QCD and Other Strong Coupled Gauge Theories” Workshop at the European Center for Theoretical Studies in Nuclear Physics and Other Related Areas (ECT*), Trento, Italy, March 20-24, 2023. A.H. acknowledges support by DOE Grant No. DE-SC0010005. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 2040434. Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy, the RMACC Summit supercomputer Anderson et al. (2017), which is supported by the National Science Foundation (awards No. ACI-1532235 and No. ACI-1532236), the University of Colorado Boulder, and Colorado State University, and the OMNI cluster of the University of Siegen. We thank BNL, Fermilab, the University of Colorado Boulder, the University of Siegen, the NSF, and the U.S. DOE for providing the facilities essential for the completion of this work.