Strong coupling from non-equilibrium Monte Carlo simulations

We compute the running coupling of non-Abelian gauge theories in the Schr\"odinger-functional scheme, by means of non-equilibrium Monte Carlo simulations on the lattice.


Introduction
During the past few years there has been significant progress towards the understanding of quantum systems out of equilibrium and of the interplay between quantum and thermodynamics effects. Research combining theoretical tools from statistical mechanics, conformal field theory, the theory of integrable systems, and quantum information has led to a deeper comprehension of the connection between entanglement entropy and thermodynamic entropy in stationary states [1], as well as a clarification of the mechanism determining the time evolution of entanglement in many-body quantum systems out of equilibrium [2].
At the same time, powerful fluctuation theorems were discovered and extensively studied in classical statistical mechanics (see refs. [3] for reviews), that encode analytical relations among quantities characterizing systems driven out of thermodynamic equilibrium. These include the transient fluctuation theorem describing the probability of violations of the second law of thermodynamics in non-equilibrium steady states [4] and the Jarzynski identity, relating the free-energy difference between two equilibrium states of a system to the exponential average of the work done on the system to drive it out of equilibrium [5].
In the present work, we show how the latter theorem can be applied to study the renormalized coupling in non-Abelian, non supersymmetric gauge theory. This quantity is of major relevance in elementary particle theory: in particular, the gauge coupling g of quantum chromodynamics (QCD) is one of the fundamental parameters in the Standard Model and plays a central rôle in theoretical predictions relevant for the physics probed in high-energy experiments 1 like those at the CERN LHC [8].
Specifically, we study the scale dependence of the gauge coupling in a non-equilibrium generalization of a Monte Carlo calculation in the lattice regularization [9] by defining the theory in a four-dimensional box of finite linear extent L, with boundary conditions enforcing a nonzero minimal-Euclidean-action configuration, and monitoring the response of the system under a sequence of quantum quenches (in Monte Carlo time) that deform the boundary conditions driving the system out of equilibrium. The formalism rests directly on the definition of the definition of the coupling in the Schrödinger-functional scheme [10,11], whereby the inverse squared physical coupling at distance L is defined to be proportional to the quantum effective action of the system with boundary conditions in Euclidean time. Here, we determine the strength of the coupling from the "stiffness" of the theory with respect to changes in the parameters that define the boundary field values, which is computed by means of Jarzynski's equality. For the sake of simplicity, the calculation is carried out in the pure-glue sector, for SU (2) and SU(3) gauge groups, and we show that the results obtained are fully compatible with previous calculations in the conventional (equilibrium) setting [12,13]. We remark that the generalization to include dynamical matter fields and/or to other non-Abelian gauge groups is straightforward.

Numerical implementation
Jarzynski's theorem [5] states that when a thermodynamic system, initially in thermal equilibrium at temperature T , is driven out of equilibrium by a time-dependent variation protocol for the parameters λ (such as couplings, etc.) of its Hamiltonian during a finite time interval [t in , t fin ], the exponential average of the work W done on the system is equal to the ratio of the partition functions (denoted by Z) for equilibrium states of the system with parameters λ(t fin ) and λ(t in ): The quantity on the left-hand side of eq. (1) is a statistical average over all possible evolutions of the system, when its parameters are modified according to a protocol λ(t), which is fixed and arbitrary.
In our calculations, we regularized the SU(N ) Yang-Mills theory (with N = 2 and 3) on a hypercubic lattice of spacing a and linear extent L in each direction. The degrees of freedom of the theory (matrices U in the defining representation of the gauge group) are associated with the oriented lattice bonds. Periodic boundary conditions are assumed along the three spatial directions, whereas fixed boundary conditions are imposed at the initial and final Euclidean time, where the fields are set to fixed, spatially uniform, Abelian matrices. The dynamics is governed by the action S = −(1/g 2 0 ) p w(p)Re Tr U p [9], where g 0 denotes the bare coupling, U p is the path-ordered product of the matrices on the a × a square ("plaquette") labelled by p. w(p) = 1 in the bulk of the system, while it equals 1/2 for spatial plaquettes on the three-dimensional slices at the initial (x 0 = 0) and final (x 0 = L) Euclidean time, and it equals c t (g 0 ) for plaquettes parallel to the Euclidean-time direction and touching the boundaries. For consistency with previous works we compare our results with, we set the "improvement coefficient" c t (g 0 ) to 1 for N = 2 [12], whereas c t (g 0 ) = 1 − 0.089g 2 0 for N = 3 [13]. For later convenience, we also define β = 2N/g 2 0 . The reformulation of eq. (1) in the Euclidean quantum-field-theory setting relevant for our Monte Carlo simulations is straightforward, with W/T replaced by the total Euclidean-action variation ∆S during each non-equilibrium trajectory of the field configuration [14]. In our calculations, λ is identified with the angle η that defines the field configurations for spatial link matrices at the boundaries, viz U = exp(iaC x 0 ) with for N = 2 and for N = 3 (in the following, we set ν = 0). Classically, this induces a spatially uniform Abelian gauge field configuration with Euclidean action for N = 2, and for N = 3. We define the evolution of λ(t) as a sequence of n qq quantum quenches in Monte Carlo time, in which η is varied from an initial value η(t in ) (equal to π/4 for N = 2, or to 0, for N = 3) to a final value η(t fin ) = η(t in ) + ∆η; for simplicity, the amplitude of these quenches is taken to be constant, ∆η/n qq . After each quench, the field configuration is changed by a Monte Carlo step (which consists of one heat-bath [15] and three to ten over-relaxation updates [16] on SU(2) subgroups [17] for all U matrices): this is done without allowing the field to thermalize, thus driving the configuration progressively out of equilibrium. We verified that a "reverse" implementation of this non-equilibrium evolution, from η(t fin ) to η(t in ), always yield consistent results: in view of the non-symmetric rôles of the initial and final states, this is a non-trivial check of the robustness of our calculation. We compute the Z λ(t fin ) /Z λ(t in ) ratio using eq. (1). Setting Γ = − ln Z, the physical coupling at the length scale L is then defined as for the SU(2) theory and in the SU(3) theory. It is worth remarking that the quality of the numerical estimate of the average on the lefthand side of eq. (1) depends crucially on how far from equilibrium the field configurations are driven during the Monte Carlo trajectories, and on the statistics of trajectories that are sampled. In a nutshell, the exponential average in eq. (1) implies that arbitrarily large deviations from equilibrium would require prohibitively large statistics to probe the tail of the ∆S distribution. The present calculation, however, does not require to probe deep out-of-equilibrium dynamics, as the physical coupling is obtained in the limit of small ∆η (and, consequently, small deviations from equilibrium). The bounds on the number of trajectories required to achieve a given level of precision in experimental or numerical sampling of out-of-equilibrium distributions are mathematically well understood [18] and are always satisfied in our Monte Carlo simulations.
Following the procedure outlined in refs. [12,13], the evolution of the physical coupling as a function of the momentum scale O(1/L) is then defined in an iterative way, in terms of the continuum-extrapolated step-scaling function σ(s, g 2 (L)) = g 2 (sL). Note that σ can be thought of as an integrated version of the β function of the theory, as it describes the evolution of the coupling between the length scales L and sL. We used s = 2 and s = 3/2.

Results and analysis
3.1 Results for the SU(2) theory We first discuss the SU(2) theory. The first step in the analysis of our numerical results consists in studying the distribution of Euclidean-action variations along the non-equilibrium trajectories. As an example, figure shows the results obtained from simulations with N = 2, L = 5a at β = 4/g 2 0 = 2.7124, for different values of ∆η and n qq = 200 quenches. We note that the numerical results can be are approximately modelled by Gaußian distributions centered at −0.156800(39) (for ∆η = 0.015), at −0.104812(26) (for ∆η = 0.01), at −0.052544(13) (for ∆η = 0.005), at −0.021054(5) (for ∆η = 0.002), and at −0.0105310(26) (for ∆η = 0.001). The width of these distributions decreases to zero with ∆η, as expected at fixed n qq . This is simply a consequence of the fact that, for very small values of ∆η/n qq , the field configurations remain close to equilibrium in every trajectory: for ∆η/n qq = 0, the simulation would reduce to a conventional equilibrium Monte Carlo. We also note that the distributions of Euclidean-action variations in reverse trajectories, from η(t in ) = π/4 + ∆η to η(t fin ) = π/4, are approximately symmetric with respect to those observed in direct trajectories.
A summary of a larger sample of our data for the SU(2) gauge theory, for different values of ∆η, and from direct and reverse implementations of our non-equilibrium Monte Carlo simulations, is reported in table 1. Here, the values of the bare coupling and system size are those in the last series of ref. [12, table 2], and the table shows the effective action difference ∆Γ = − ln[Z η(t fin ) /Z η(t in ) ]. The table reveals clearly that direct and reverse simulations yield compatible results, and that ∆Γ scales linearly with ∆η. This makes the evaluation of g 2 (L) from eq. (6) robust and unambiguous. In view of these results, in order to reduce the computational cost of several statistically independent computations at different values of ∆η, we then proceeded to the calculation of g 2 (L) from the results for ∆Γ obtained at ∆η = 0.0001, increasing n qq to 1000.
The results of this computation are reported in tables 2 and 3. The bare couplings and system sizes reproduce those that were analyzed in ref. [12] and the results obtained from our non-equilibrium Monte Carlo calculations are fully compatible with those reported in that work. We observe that results obtained from "direct" and "reverse" implementations of our algorithm are consistent with each other: a non-trivial check that our algorithm is correctly sampling the distribution of Euclidean action differences along the non-equilibrium trajectories. Our final results for the squared coupling at different lattice spacings are then obtained from the average of the two.
In (2). The histograms show the distribution of ∆S induced by a "direct" non-equilibrium transformation in which η is varied from η = π/4 to η = π/4 + ∆η through a sequence of n qq = 200 quenches, for different values of ∆η. The larger inset shows the same distributions using a logarithmic scale for the vertical axis. In the smaller inset, the results from "reverse" transformations, from η = π/4 + ∆η to η = π/4.   (7) 3.7019 (7) 3    Continuum extrapolation of the results in tables 2 and 3 reveals that the value of the squared coupling in the Schrödinger-functional scheme evaluated on the largest lattices is g 2 = 4.85(4), to the continuum limit by a constant-plus-linear-term fit in a/L, and are compatible with those reported in ref. [12]. and that the values of g 2 (2L) obtained in the a → 0 limit from each of the four data sets are close to the value of g 2 (L) in the next set. As a consequence, the ratio of the corresponding length scales is close to unity, and can be reliably estimated using the perturbative β function truncated at two loops. This procedure allows one to determine the evolution of the coupling from the "hadronic" scale down to the microscopic scale, where perturbation theory becomes reliable, using the step-scaling function defined as σ s, g 2 (L) = g 2 (sL) (8) and constructing the β function of the theory in a recursive way. This procedure allows one to follow the evolution of the coupling in the Schrödinger-functional scheme non-perturbatively over a wide range of scales, without having to perform simulations on lattices with a prohibitively large number of sites. The last step of the analysis consists, then, in the explicit construction of the β function of the theory. To this purpose, we focus on the low-energy regime. We run an additional set of simulations for (β, L/a) combinations yielding a value of g 2 (L) sufficiently close to the one extrapolated from the largest lattices listed in table 3 i.e. g 2 = 4.85 (4) We can then make contact with a physical low-energy scale of the theory, such as the string tension σ 0 , i.e. the asymptotic force between fundamental probe charges at large distance, 3 using the data reported in refs. [12,[20][21][22] in the range β ∈ [2.2, 2.85]. In particular, using the value σ 0 a 2 = 0.00830(6) at β = 2.74 from ref. [21], one obtains that the length at which g 2 equals 4.85 is L = 0.843(3)/ √ σ 0 , or 0.3781 (14) fm. Thus, taking the momentum scale to be defined as µ = 1/L and using the continuum extrapolations of the results listed in tables 2 and 3, one obtains the results for α s = g 2 /(4π) plotted in figure 4. Also shown are the analytical predictions from perturbation theory at one and two loops [6,23]: which yields ln Our numerical results are in very good agreement with those reported in ref. [12], and confirm the accuracy of the two-loop perturbative β function down to µ ∼ 1 GeV. At the three lowest energy scales the two-loop perturbative prediction systematically underestimates the non-perturbative Monte Carlo results. As a curiosity, extrapolating our results to high energies using the two-loop perturbative β function, at the pole mass of the physical Z 0 boson of the Standard Model we obtain α s (m Z 0 ) = 0.1081(6) in the Schrödinger functional scheme for the SU(2) purely gluonic Yang-Mills theory.

Results for the SU(3) theory
Our study of the SU(3) theory follows closely the one presented in subsection 3.1 for the N = 2 case. In this case, we compare our results with those reported in ref. [13]. The main difference with respect to the SU(2) case is that the rank of the algebra of group generators is 2 (instead of 1) and the fundamental domain is specified by the two parameters η and ν (instead of just η). We run our non-equilibrium simulations around η = 0 (instead of η = π/4), at fixed ν = 0. In addition, as mentioned earlier, for the SU(3) theory we set c t , the improvement coefficient for plaquettes parallel to the Euclidean-time direction touching the boundaries, to 1 − 0.089g 2 0 (instead of 1).
The results reported in tables 6, 7, and 8 are plotted against sL/a in figure 5, which also shows their extrapolation to the continuum limit, and the comparison with the two-loop perturbative predictions.
As for the analysis of the N = 2 theory, we note that the value of the squared coupling on the largest system that we simulated is g 2 = 3.467 (15). We therefore run an additional set of simulations that approximately correspond to this value of the squared coupling, for the (β, L/a) combinations listed in table 9. Then, we estimate the actual values of β that would yield g 2 = 3.467 using the two-loop approximation in eq. (10) and the parametrization of the lattice spacing as a function of β reported in ref. [24, eq. (2.6)], which holds in the whole range of β values reported in table 9, and take the difference with respect to the simulated values of β as an estimate of the uncertainty on the results. This leads to the values listed in table 10, which can be fitted to the functional form β = 4.797(6) + 0.798(4) · ln(L/a); both the data points and the fitted curve are plotted in figure 6.
Finally, we can then match our low-energy results with a physical scale of the theory. In this case, we choose Sommer's parameter r 0 , defined as the distance r at which the force F between fundamental probe charges satisfies r 2 F (r) = 1.65 [25]. Using the high-precision parametrization of the lattice spacing in units of r 0 that was reported in ref. [24]: ln a r 0 = −1.6804 − 1.7331 · (β − 6) + 0.7849 · (β − 6) 2 − 0.4428 · (β − 6) 3 , for 5.7 ≤ β ≤ 6.92, (11) we deduce that the value of the lattice spacing at β = 6.08169 equals a = 0.1625157 · r 0 , and that, as a consequence, L = 5a = 0.8125786 · r 0 . We can convert this into physical units by taking the estimate for the physical value of r 0 in QCD r 0 = 0.468(4) fm from ref. [26] 4 , which leads to L = 0.3803(33) fm and q = 1/L = 0.519(4) GeV. Proceeding as for the SU(2) theory, we finally obtain the results plotted in figure 7. As in the N = 2 case, one can see that our numerical results reproduce those obtained from "conventional" (equilibrium) Monte Carlo calculations [13]. Our results are in very good quantitative agreement with the analytical prediction from the two-loop perturbative β function, eq. (10): this holds for all values of µ, down to µ ∼ 1 GeV. In this theory, an extrapolation results for mesons [19], is approximately (440 MeV) 2 . β type L/a n traj (L) g 2 (L) sL/a n traj (sL) g 2 (sL)  Table 6: Results for g 2 (L) and g 2 (sL) from direct and reverse transformations with ∆η = 0.0001 and n qq = 1000 in SU(3) Yang-Mills theory, and their average. β type L/a n traj (L) g 2 (L) sL/a n traj (sL) g 2 (sL)     L/a n traj (L) g 2 (L) sL/a n traj (sL) g 2 (sL)  of α s to the pole mass of the Z 0 boson of the Standard Model yields α s (m Z 0 ) = 0.07297 (19). When converted to the modified minimal-subtraction (MS) scheme via the one-loop relation α MS s = α s + 1.25563(4)α 2 s [13], this corresponds to approximately 0.07966 (22). For comparison, note that at this scale the value of α s in the modified minimal-subtraction scheme recently reported in ref. [29] for QCD with five quark flavors is α MS s = 0.11852(84) (for further technical details, see also refs. [30] and the references therein).
Finally, a comparison of the results for α s (m Z 0 ) that we obtained in the theories with N = 2 and N = 3 shows that N α s (m Z 0 ) is independent of N , to a sub-percent level of precision. This relation has a natural interpretation in terms of the 't Hooft coupling in the large-N limit of QCD [31] (see also refs. [32]), and it is unsurprising that it holds even for values of N as small as N = 2 [33].

Conclusions
In this work, we presented the results of a non-perturbative study of the running coupling of non-Abelian gauge theories, by means of a non-equilibrium Monte Carlo algorithm that implements a numerical realization of Jarzynski's theorem. Specifically, we evaluated the response in effective action induced by a deformation of the boundary conditions at the initial and final Euclidean time, and extracted the running coupling in the Schrödinger-functional scheme [11].
The latter scheme provides a well-defined formulation of the theory in a finite system, with Dirichlet boundary conditions along Euclidean time and periodic (or periodic up to a constant phase, for fermionic fields) boundary conditions along the three spatial directions, and allows one β type L/a n traj (L) g 2 (L)   to define the renormalized coupling at a momentum scale defined as the inverse of the linear extent of the system in each direction, µ = 1/L. This formulation is amenable to lattice regularization and, through the iterative procedure that we discussed in section 3, it allows one to study the evolution of a renormalized quantity when the momentum scale varies by orders of magnitude, by recursively matching the (continuum-extrapolated) value of the coupling obtained on lattices of the same physical extent, but different lattice spacing. The formalism has been successfully applied to study the renormalized gauge coupling and a number of other physical quantities in QCD [34] and, more recently, has also been used to investigate the dynamics of other strongly coupled non-supersymmetric gauge theories coupled to fermionic fields [35]. The goal of this work consisted in showing that the Schrödinger-functional formalism can be directly implemented in Monte Carlo calculations out of equilibrium, using the powerful fluctuation theorems that have been recently developed in statistical mechanics. As discussed in section 2, the central idea is to drive the system out of equilibrium through a sequence of "quantum quenches in Monte Carlo time": Jarzynski's theorem [5] implies that the exponential average of the Euclidean-action variation induced in this process is directly related to the exponential of the difference in effective action between the initial and the final states of the system.
We emphasize that, while in our computation we evaluate a quantity (the discretized derivative of the effective action with respect to η, the parameter that specifies the Dirichlet boundary conditions of the system at the initial and final Euclidean times) which is directly related, and can be made arbitrarily close, to the one that is evaluated in conventional simulations of the Schrödinger functional, the approach is intrinsically radically different, as our calculation does not rest on the standard formalism of equilibrium equilibrium-Monte Carlo. In particular, our approach closely "mimics" the non-trivial dynamics that is induced in physical statistical systems that are experimentally driven out of equilibrium and the corresponding measurements that can be performed on them. Experimental applications of this type are diverse, and include, for example, irreversible mechanical stretching of ribonucleic acid molecules [36].
We focused on SU(2) and SU(3) Yang-Mills theories, and showed that the results obtained in our non-equilibrium Monte Carlo simulations are fully compatible with those from standard (equilibrium) lattice simulations [12,13]. While we presented results for purely bosonic theories, the generalization of this calculation to include dynamical fermions poses no additional computational challenge, and could be easily carried out with the same techniques common to lattice QCD, e.g. through (a non-equilibrium version of) the hybrid Monte Carlo algorithm [37].
Finally, we mention some other recent, and very interesting, articles that present applications of non-equilibrium statistical-mechanics theorems in a context relevant for quantum field theory [38][39][40][41]; in particular, refs. [38,41] focus on the calculation of the entanglement entropy from Jarzynski's equality, whereas refs. [40,41] discuss the implications of non-equilibrium theorems for quantum field theories. We expect that many more such studies, at the interface between modern statistical mechanics and quantum field theory, will appear in the near future, and that they may lead to new insights into open problems both in condensed matter theory and in elementary particle theory.