Strong coupling from non-equilibrium Monte Carlo simulations

We compute the running coupling of non-Abelian gauge theories in the Schrödinger-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,2], as well as a clarification of the mechanism determining the time evolution of entanglement in many-body quantum systems out of equilibrium [3].
At the same time, powerful fluctuation theorems were discovered and extensively studied in classical statistical mechanics (see refs. [4][5][6] 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 [7][8][9][10] and Jarzynski's 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 [11,12].
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 [16][17][18][19][20]. 1 The most striking feature of the physical QCD coupling is its dependence on the momentum scale µ: the dimensionless parameter αs = g 2 /(4π) is a decreasing function of µ [13,14], so that QCD becomes a free theory at asymptotically high energies, while its behavior at low energies is non-perturbative. Note that the logarithmic running of the strong coupling is such that QCD remains a self-consistent theory for arbitrarily high energies [15]: a behavior remarkably different from other theories, like quantum electrodynamics, which break down at some high, but finite, energy scale.

JHEP07(2020)233
Specifically, we study the scale dependence of the gauge coupling in a non-equilibrium generalization of a Monte Carlo calculation in the lattice regularization [21] by defining the theory in a four-dimensional box of finite linear extent L, with boundary conditions enforcing a non-trivial 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 terminology that we are using here is inspired by condensed-matter theory, where quantum quenches are a convenient tool to study systems driven out of equilibrium. A quantum quench is defined as a sudden change of the Hamiltonian of a quantum manybody system [22,23]; before the quench, the system is in the ground state of the initial Hamiltonian, and the dynamical evolution after the quench is unitary. The system response to a quantum quench allows one to study many interesting aspects of its dynamics, including those related to localization [24], thermalization [25], the interplay between integrable and non-integrable dynamics [26], and entanglement entropy [27]. What makes quantum quenches particularly interesting is that, beside their theoretical interest, they can also be realized experimentally in certain condensed-matter systems [28].
Here, we apply an analogous idea in the context of numerical simulations of a non-Abelian gauge theory regularized on the lattice; instead of changing the (bulk) Hamiltonian, we change the Dirichlet boundary conditions that are imposed on the fields along one of the four Euclidean directions of the system, and, instead of studying the real-time evolution induced by this change, we study its evolution in Monte Carlo time. For technical reasons (related to the algorithm efficiency, to be discussed below), we apply not only one, but a sequence of such quenches. We then extract the physical gauge coupling at the length scale defined by the system size from the response of the system (as encoded in its quantum effective action) under the sequence of quenches that drives it out of equilibrium. The formalism rests directly on the definition of the coupling in the Schrödinger-functional scheme [29,30], whereby the inverse squared physical coupling at distance L is given (up to a normalization factor) by the derivative of the quantum effective action with respect to the parameter, to be denoted as η, specifying the Dirichlet boundary conditions. In other words, the Schrödinger-functional coupling is defined as a coefficient that encodes the response of the theory to a variation of the field enforced by the boundary conditions. It is important to note that, while the quantum effective action of the theory cannot be directly accessed in Monte Carlo calculations (neither by conventional algorithms, nor by the one we discuss in this work), its derivative with respect to η is a "measurable" quantity in a simulation, and, as we will discuss in detail below, our algorithm estimates numerically precisely this quantity, 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 [31,32]. For the SU(2) theory, another relevant work was reported in ref. [33], which studied the approach to the continuum with and without boundary improvement terms. We remark that the generalization to include dynamical matter fields and/or to other non-Abelian gauge groups is straightforward.

JHEP07(2020)233 2 Numerical implementation
Jarzynski's theorem [11,12] 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 H during a finite time interval [t in , t fin ], the exponential average of the work W done on the system in units of the temperature 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. (2.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.
To clarify the meaning of the average appearing on the left-hand side of eq. (2.1), and to describe how our algorithm works, it is convenient to review the proof of eq. (2.1) in a setup that is relevant for our calculation. Note that the time t that appears in Jarzynski's theorem can be either real time or Monte Carlo time; in the following, we identify t with Monte Carlo time. We remark that the identity encoded in eq. (2.1) is valid under general conditions; in particular, it holds when the starting configurations in the evolution of the system are at equilibrium and the dynamics of the system preserves the equilibrium distribution if the parameters λ are fixed. To prove eq. (2.1) for a statistical system undergoing Monte Carlo evolution which satisfies the stronger detailed-balance condition (like the lattice version of the gauge theories considered in this study), let φ denote the degrees of freedom of the system, let be the normalized equilibrium distribution corresponding to a given, fixed value of λ (and Z λ the corresponding partition function), and let P λ (φ → φ ) be the normalized transition probability density from a given configuration φ to another configuration φ at fixed λ. The detailed-balance condition reads In general, the work done on the system when the parameters are varied according to the λ(t) protocol (and, as a consequence, the fields undergo a non-equilibrium evolution φ(t)) is whereλ is the derivative of λ with respect to t. In the rest of this article, following the terminology of refs. [11,12], we call each non-equilibrium evolution of the fields in Monte Carlo time φ(t) a non-equilibrium Monte Carlo trajectory; as stated above, during their evolution the fields do not remain in equilibrium, because the parameters λ are varied JHEP07(2020)233 as a function of Monte Carlo time. Note that, in this context, the concept of trajectory is distinct from the one of a trajectory of the hybrid Monte Carlo algorithm [34][35][36][37][38], which is widely used in the lattice QCD literature (even though, as will be briefly mentioned in section 4, the non-equilibrium algorithm discussed in this work could also be implemented in a non-equilibrium version of the hybrid Monte Carlo algorithm). Let us assume that each Monte Carlo trajectory is made of n steps, and that at each step along a trajectory, λ is first updated, then the field configuration is let evolve according to the dynamics defined by the new value of λ. Let i label Monte Carlo time (with i = 0 corresponding to the initial time t in and i = n corresponding to the final time t fin ), so that φ i denotes the field configuration at Monte Carlo time labeled by i. The work done on the system when λ is varied, say, from λ i to λ i+1 (and before the Monte Carlo update of the field), is the difference Then, using eq. (2.2), the exponential of minus the work divided by T , which appears in the average on the left-hand side of eq. (2.1), can be rewritten as On the left-hand side of eq. (2.1), this quantity is averaged over all trajectories that φ can span in its Monte Carlo evolution: this corresponds to integrating over all possible configurations at each step in the Monte Carlo trajectory. The initial configurations are distributed according to the equilibrium distribution π λ 0 , while those at later Monte Carlo times are further weighted by products of the transition probabilities P λ i+1 (φ i → φ i+1 ). Thus, the left-hand side of eq. (2.1) can be expressed as exp (−W/T ) = dφ 0 dφ 1 dφ 2 · · · dφ n π λ 0 (φ 0 ) (2.6) All partition functions appearing in the product cancel against each other, except for the one appearing in the denominator of the fraction in the first term, and the one in the numerator in the last factor, leaving Z λn /Z λ 0 . Using eq. (2.3), the previous equation can be rewritten as (2.7) so that also the π distributions cancel against each other, except for π λn (φ n )/π λ 0 (φ 0 ). Multiplying this quantity by the π λ 0 (φ 0 ) term in front of the product, one is left with In the latter expression, φ 0 appears only in the P λ 1 (φ 1 → φ 0 ) factor, and the normalization of the transition probabilities implies dφ 0 P λ 1 (φ 1 → φ 0 ) = 1. The same argument can then be repeated for the integrals over φ 1 , φ 2 , . . . , φ n−1 . Finally, the last integral is JHEP07(2020)233 dφ n π λn (φ n ), which also equals one, because the equilibrium distributions are normalized, and one arrives at eq. (2.1).
Let us now discuss how we implemented eq. (2.1) in our algorithm for the numerical evaluation of the running coupling in the Schrödinger-functional scheme. Following refs. [30][31][32], 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 links. Periodic boundary conditions are assumed along the three spatial directions, whereas fixed boundary conditions are imposed at the initial (x 0 = 0) and final (x 0 = L) Euclidean time, where the spatial links are set to fixed, spatially uniform, Abelian matrices defined below, while no boundary conditions are imposed on the temporal links between sites on the boundaries and sites in the bulk of the lattice (and there are no positively oriented temporal links from the sites on the boundary at Euclidean time x 0 = L, nor negatively oriented temporal links from the sites on the boundary at Euclidean time x 0 = 0). The dynamics is governed by the action S = −(1/g 2 0 ) p w(p) Re Tr U p [21], where g 0 denotes the bare coupling, and U p is the path-ordered product of the matrices on the a × a square ("plaquette") labeled 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 Euclidean times x 0 = 0 and x 0 = L, and it equals c t (g 0 ) for plaquettes parallel to the Euclidean-time direction and touching the boundaries. For consistency with the previous works we compare our results with, we set the "improvement coefficient" c t (g 0 ) to 1 for N = 2 [31], whereas c t (g 0 ) = 1 − 0.089g 2 0 for N = 3 [32]. For later convenience, we also define β = 2N/g 2 0 . The reformulation of eq. (2.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 [39,40]. Note that, since each non-equilibrium trajectory is decomposed into n steps, so is the Euclidean-action variation ∆S: more precisely, (2.9) 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 S cl = 24L 4 g 2 0 a 4 sin 2 a 2 2L 2 (π − 2η) (2.12)

JHEP07(2020)233
for N = 2, and We define the evolution of λ(t) as a sequence of n qq = n 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 [41,42] and three to ten over-relaxation updates [43,44] on SU(2) subgroups [45] 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. (2.1). Setting Γ = − ln Z, the physical coupling at the length scale L is then defined as the ratio between the derivative of g 2 0 S cl with respect to η and the derivative of Γ with respect to η. In turn, the latter is given by the limit of the difference quotient ∆Γ/∆η for ∆η → 0, so that one obtains for the SU(2) theory (having set η = π/4) and It is worth remarking that the quality of the numerical estimate of the average on the left-hand side of eq. (2.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. (2.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 [46][47][48][49][50] and are always satisfied in our Monte Carlo simulations.
Following the procedure outlined in refs. [31,32], 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) that was introduced in ref. [51]. 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 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 1 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 approximately modeled 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 . As will be discussed in detail in subsection 3.3, 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 more detailed analysis of the results displayed in figure 1, providing information about the efficiency of our numerical algorithm, is presented in subsection 3.3, together with a comparison with the computational costs of lattice calculations of the running coupling in the Schrödinger-functional scheme by means of conventional, equilibrium Monte Carlo algorithm.
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. [31, 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 non-equilibrium transformations yield compatible results, and that ∆Γ scales linearly with ∆η. Figure 2 shows the results for ∆Γ/∆η against ∆η, as obtained from direct transformations at the different β values reported in table 1, corresponding to different values of the lattice spacing, for approximately constant physical linear size of the lattice: the plot reveals that the difference quotient remains essentially constant for all values of ∆η. The figure shows the consistency of the results at different β values, and a mild trend towards values of ∆Γ/∆η that are slightly more negative (i.e. larger in modulus) when ∆η is reduced towards zero, albeit only by an amount (with respect to the smallest ∆η considered here, i.e. ∆η = 0.001) that is comparable with our statistical uncertainties. This makes the evaluation of g 2 (L) from eq. (2.14) 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 (a value ten times smaller than the smallest one used to produce the data sets reported in table 1 and shown in figure 2), increasing n qq to 1000.
The results of this computation are reported in tables 2 and 3. The bare couplings and system sizes are the same that were analyzed in ref. [31] and the results obtained from JHEP07(2020)233  Figure 1. Distribution of the Euclidean action difference ∆S in non-equilibrium simulations of the SU(2) gauge theory in a hypercubic box of linear size L = 5a at β = 2.7124, with the boundary fields specified in eq. (2.10). The histograms show the distribution of ∆S induced by a "direct" nonequilibrium 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, are displayed. our non-equilibrium Monte Carlo calculations are fully compatible with those reported in that work, which, for the reader's convenience, we reproduce in table 4. We also 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 figure 3 we show our results for g 2 (2L) at different values of the lattice spacing. The data, displayed by red squares, fall on four nearly horizontal bands, corresponding to four values 2 of g 2 (L) = 2.059(11), 2.353(4), 2.871 (14), and 3.546(16), i.e. to four values of L, and are plotted as a function of the lattice spacing divided by 2L. For comparison, the figure also shows the results reported in ref. [31] as black circles. Since leading discretization 2 These values are obtained from the averages reported in the fifth column of tables 2 and 3 after extrapolation to the continuum limit by a constant-plus-linear-term fit in a/L, and are compatible with those reported in ref. [ effects in the lattice formulation of the Schrödinger functional are expected to be of order a, we fit each of the four data sets to the sum of a constant plus a linear function of the lattice spacing, obtaining the continuum-extrapolated results displayed by the red squares on the vertical axis of the plot: following the analysis carried out in ref. [31], we find that all of them are very close to the two-loop perturbative predictions (horizontal blue segments on the vertical axis). One could also compare these results with threeloop perturbative predictions, which have since become available [52,53] and which are discussed below, but at this stage we limit ourselves to note that, as was observed in JHEP07(2020)233 β type L/a n traj (L)  JHEP07(2020)233  (7) 3.7019 (7) 3   ref. [31], the two-loop perturbative predictions already provide a good approximation for the continuum-extrapolated non-perturbative results shown in figure 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), 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 very close to unity, and, following ref. [31], can be reliably estimated using the perturbative β function truncated at two loops. As mentioned above, in principle, the estimate of this length-scale ratio could now be refined using the three-loop perturbative β function (or, alternatively, in a fully non-perturbative way, by additional sets of simulations), but, due to the smallness of the differences between the length scales, this would not yield significantly different results. Hence, at this step we applied exactly the same procedure that was carried out in ref. [31] (which is justified, given the exploratory nature of the present study, that does not attempt to produce results of direct phenomenological relevance), postponing the detailed discussion of the three-loop perturbative prediction and its comparison with our lattice results to the final part of this subsection.
The procedure outlined above, i.e. following the variation of the coupling through a sequence of lattices whose linear sizes are in a ratio s, allows one to determine the evolution JHEP07(2020)233  of g 2 from the "hadronic" scale down to the microscopic scale, where perturbation theory becomes reliable, using the step-scaling function first introduced in ref. [51] and defined as σ s, g 2 (L) = g 2 (sL) (3.1) and to extract the β function of the theory from it. It is then possible to obtain 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 first focus on the low-energy regime and 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): the results are shown in table 5. Note that, while the number of lattice points in each direction in these lattices varies from 8 to 12, the corresponding lattice spacings decrease, and the Schrödinger-functional coupling remains nearly constant, at a value that is the largest among those that we considered for this theory. As a consequence, these lattices have  nearly the same physical size, which is the largest among those that we studied for the SU(2) theory in this work. More precisely, from the values in table 5, we obtain the (β, L/a) combinations corresponding to g 2 = 4.85(4) that are listed in table 6. They can be fitted to the functional form β = 1.866(22)+0.3928(9)·ln(L/a), with reduced χ 2 0.15, as shown in figure 4.
For comparison, in table 7 we also reproduce the analogous values obtained in ref. [31] for the lattice of largest physical size studied in that work, corresponding to g 2 = 4.765 (a value close to the squared coupling obtained from continuum extrapolation of the stepscaling function evaluated on the largest set of lattices considered therein, which reads g 2 = 4.76 (12) and is fully compatible with our result).
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. [31,[55][56][57] in the range β ∈ [2.2, 2.85]. In particular, using the value σ 0 a 2 = 0.00830(6) at β = 2.74 from ref. [56], one obtains that JHEP07(2020)233 β type L/a n traj (L) g 2 (L)      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 5. Also shown are the analytical predictions from perturbation theory at one, two, and three loops [13,14,52,53,58,59]. In particular, the two-loop perturbative prediction is obtained from dα s d(ln µ) which yields Note that the comparison with perturbation theory in the Schrödinger-functional scheme can now be pushed to three loops (which was not yet possible at the time of publication of ref. [31]): this can be done by combining the two-loop relation between the Schrödinger-functional coupling and the bare lattice coupling that was worked out in JHEP07(2020)233 ref. [52] with the one relating the bare lattice coupling to the coupling in the MS-scheme [53], from which one can derive that, for the SU(2) Yang-Mills theory, where α s and α MS can be defined at two different momentum scales, respectively µ 1 and µ 2 , and c 1 and c 2 are functions of their ratio r = µ 2 /µ 1 : Note that eq. (3.5) can be inverted as In turn, eqs. (3.5) and (3.7) can be combined with the three-loop perturbative expression for the β function in the MS scheme [60] (see also refs. [61,62]), which, for a purely gluonic SU(N ) gauge theory, reads to obtain the three-loop perturbative β function for the SU(2) Yang-Mills theory in the Schrödinger-functional scheme: Integrating eq. (3.9) (assuming r = 1), one obtains again an expression like the one given in eq. (3.3), but with eq. (3.4) replaced by (3.11) Our numerical results are in very good agreement with those reported in ref. [31], and confirm the accuracy of the three-loop (and two-loop) perturbative β function down to µ ∼ 1 GeV. For a comparison, in table 8 we reproduce the values of the running coupling in the Schrödinger-functional scheme, as a function of L/L 8 (with L 8 denoting the length scale at which g 2 (L) reaches the largest value studied in that work), that were obtained in ref. [31]. Note that the largest value of g 2 that was considered in ref. [31] (namely 4.765) is close but not exactly equal to ours (which is 4.85(4)). As remarked above, this mismatch arises at an intermediate step in the calculation, in particular at the level of the continuum extrapolation of the results for g 2 (2L), which (as shown for instance in figure 3) are affected by different statistical uncertainties in the two studies. Note, however, that the g 2 value JHEP07(2020)233 that was obtained from the continuum extrapolation of the step-scaling function on the set of largest lattices in ref. [31] is 4.76 (12), in perfect agreement with our value within one standard deviation. The overall agreement between the two works remains very good.
At the three lowest energy scales that we studied, the two-loop perturbative prediction systematically underestimates the non-perturbative Monte Carlo results, while the threeloop perturbative prediction is in agreement with them. As a curiosity, extrapolating our results to high energies using the two-loop (or the three-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 purely gluonic SU(2) 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. [32]. 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). As before, we study the running coupling in the Schrödinger-functional scheme through simulations on lattices of physical linear sizes L and sL, where s = 2 throughout, except for the last set of four simulations, corresponding to g 2 (L) = 2.77, for which s = 3/2. We emphasize that this choice is simply due to our decision of reproducing the exact choice of parameters as in ref. [32], for an easier comparison with that reference work, and is not due to any inherent limitation of our algorithm.
Our results, reported in tables 9, 10, and 11, are plotted against a/(sL) in figure 6, which also shows their extrapolation to the continuum limit, and the comparison with the two-loop perturbative predictions.
For comparison, we also reproduce the results obtained with a conventional Monte Carlo algorithm in ref. [32]   As for the analysis of the N = 2 theory, we focus on the value of the squared coupling on the largest system that we simulated, which in this case 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 13. Then, we estimate the actual values of β that would yield g 2 = 3.467 using the two-loop approximation in eq. For comparison, in table 15 we also reproduce the analogous results from ref. [32, table 3], which correspond to a similar value of the squared coupling: g 2 (L) 3.48.
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 [64]. Using the high-precision JHEP07(2020)233 β type L/a n traj (L) g 2 (L) sL/a n traj (sL) g 2 (sL)  Table 9. 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. parametrization of the lattice spacing in units of r 0 that was reported in ref. [63]: 12) 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. [65], 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 can then compare our lattice results with perturbative predictions.

JHEP07(2020)233
β type L/a n traj (L) g 2 (L) sL/a n traj (sL) g 2 (sL)     For the theory with N = 3 color charges, the perturbative three-loop β function in the Schrödinger-functional scheme was worked out in ref. [68] and reads Integrating eq. (3.13), we finally obtain the curves plotted in figure 8 alongside our simulation results. Also in this case, we can compare our results with those from conventional Monte Carlo simulations in ref. [32]: table 16 reproduces the values for g 2 (L) that were obtained in that work at different length scales L, in units of L max , defined as the box size at which g 2 3.48.
While the values of the squared Schrödinger-functional coupling corresponding to the lattice with the largest physical size obtained in our work and in ref. [32] are slightly different (a mismatch due to the different level of precision of the extrapolations involved),

JHEP07(2020)233
β type L/a n traj (L) g 2 (L)     we see that the behavior of the running coupling obtained in the two works with two radically different numerical algorithms is quantitatively consistent. We conclude that, as in the N = 2 case, our numerical results reproduce those obtained from "conventional" (equilibrium) Monte Carlo calculations [32]. Our results are in very good quantitative agreement with the analytical predictions from the three-loop perturbative β function, eq. (3.13) (and also with its truncation at two loops): this holds for all values of µ, down to µ ∼ 1 GeV. In this theory, an extrapolation 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 + 1.25563(4)α 2 s [32], 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. [69] for QCD with five quark flavors, combining the three-flavor running with a perturbative matching across the thresholds corresponding to the charm-and bottom-quark masses, where perturbation theory works well, is α MS = 0.11852(84) (for further technical details, see also refs. [70,71] and the references therein). For a review of lattice results on α s , see ref. [72].

JHEP07(2020)233
Finally, a comparison of the results for α s (m Z 0 ) that we obtained in the theories with N = 2 and N = 3, when taken at face value, 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 [73] (see also refs. [74,75]), and it is unsurprising that it holds even for values of N as small as N = 2 [76]. However, it is worth remarking that, even if one only compares purely gluonic theories, this statement may still depend on the way the physical scale is set.

Computational efficiency analysis
In this subsection, we assess the computational efficiency of our algorithm, and compare it with conventional (equilibrium) Monte Carlo algorithms for the calculation of the running coupling in the Schrödinger-functional scheme. We first discuss general features and expectations, then we present quantitative details obtained from the analysis of our results.
The first important observation is that, at least superficially, our numerical determination of the running coupling follows quite closely the computation that is carried out in conventional simulations, in the sense that it is based on the same quantity, i.e. the derivative of the effective action Γ with respect to η. In conventional algorithms, this derivative (which is nothing but a sum of chromoelectric field components evaluated on the Euclidean-time-slices x 0 = 0 and x 0 = L: see, e.g., ref. [31, eq. (3.2)]) is evaluated directly, on an ensemble of equilibrium configurations. The calculation boils down to computing a sum of local objects (traces of plaquette operators with the insertion of a generator of the gauge-group algebra, as explicitly detailed in ref. [31, eq. (3.3)]) on the two boundaries of the system at the initial and final Euclidean times. Correspondingly, in our calculation, we estimated the effective-action derivative from the ∆η → 0 limit of the ∆Γ/∆η quotient ratio, extracting ∆Γ from eq. (2.1), which requires the numerical evaluation of the variation in Euclidean action ∆S in each non-equilibrium trajectory. In turn, the variation in Euclidean action induced by a sequence of changes in η is simply expressed in terms of the difference in the plaquette values at the boundaries.
Note, however, that, as we remarked in section 2, the total Euclidean-action variation during each trajectory is decomposed into the sum of n qq terms (each of them corresponding to the Euclidean-action variation induced by a single "quantum quench"), hence one might be tempted to conclude that our non-equilibrium algorithm becomes exceedingly more timeconsuming than a conventional Monte Carlo simulation for large values of n qq (as discussed above, we used n qq = 1000 for most of our production runs). This argument, however, is misleading, because it does not take the distribution of values of the observable into account.
To clarify this point, it is enlightening to discuss our non-equilibrium algorithm in two different, and particularly interesting, limits.
Firstly, we note that, for n qq = 1, η would be immediately switched from its initial to its final value: then, according to eq. (2.9), our algorithm would simply use a set of equilibrium configurations generated at η = η(t in ) to compute the exponential average of the difference in Euclidean action obtained by changing the boundary fields from η = η(t in ) to η = η(t fin ). Then, eq. (2.1) would reduce to 14)

JHEP07(2020)233
with the average on the left-hand side evaluated on the equilibrium configurations generated at η = η(t in ). This would then correspond to considering the exponential of the Euclidean-action difference between the two ensembles as an observable to be evaluated on the ensemble with η = η(t in ), i.e. to an implementation of reweighting [77]. Reweighting is a computational technique that is widely used in Monte Carlo simulations of statisticalmechanics systems, whose derivation does not involve any non-equilibrium assumptions, and whose computational efficiency may be limited by the existence of a (possibly severe) overlap problem: when the typical configurations contributing to Z η(t in )+∆η and Z η(t in ) are very different, the evaluation of the left-hand side of eq. (3.14) requires exceedingly large statistics. Note that one of the works we compare our results with, namely ref. [32], pointed out the existence of a long tail in the distribution of the derivative of Γ with respect to η at the largest coupling considered there, leading to long auto-correlation times for this observable, and presented a detailed discussion of how this numerical problem was tackled using the reweighting technique. Secondly, we note that in the opposite limit n qq → ∞, the η parameter would be varied infinitely slowly. As a consequence, the system would remain in equilibrium throughout the Monte Carlo evolution: from a statistical-mechanics viewpoint, this would then imply that the work W done on the system during each trajectory would be exactly equal to the difference between the free energies of the final and of the initial ensemble. In particular, this would also imply that the width of the distribution of W/T values contributing to the average on the left-hand side of eq. (2.1) would tend to zero (so that, in principle, one could determine the Z λ(t fin ) /Z λ(t in ) ratio from just one trajectory).
We now proceed to a more quantitative study of these aspects, presenting a detailed analysis of our lattice results. We start from the simulation ensembles whose results are plotted in figure 1 and summarized in table 1: as discussed in subsection 3.1, they are obtained from non-equilibrium simulations of SU(2) Yang-Mills theory on a lattice with L/a = 5 at fixed β = 2.7124 and for fixed n qq = 200, from η(t in ) = π/4 to η(t fin ) = η(t in )+∆η, for the runs labeled as "direct" transformations (or vice versa, for those denotes as "reverse" transformations), for different values of ∆η. As mentioned in subsection 3.1, the distributions of Euclidean-action differences ∆S observed along the non-equilibrium trajectories spanned in those runs are approximated well by Gaußian distributions. This is made manifest by the analysis of the first few moments and cumulants of these distributions (up to the fourth order): the mean, the variance, the standardized skewness and the excess kurtosis. For a real-valued variable X with normalized probability distribution P(X), they are respectively defined as   We estimated these measures for the Euclidean-action differences ∆S observed in our nonequilibrium simulations, obtaining the results reported in table 17, where the values that are statistically compatible with zero are denoted by a star.
The table shows very clearly that, as expected, the width of the distribution of ∆S values obtained along non-equilibrium trajectories shrinks to zero when ∆η/n qq tends to zero at fixed n qq : as anticipated, this is the limit in which the field configurations spanned during the trajectories do not depart strongly from equilibrium, hence all values of the Euclidean-action difference measured numerically by our algorithm are very close to each other, and, according to eq. (2.1), to the effective-action difference that is induced by that ∆η. The fact that these trajectories are indeed quite close to equilibrium is also confirmed by another observation, namely that the Euclidean-action-difference distributions observed in "reverse" transformations (with statistical measures on the same lines of table 17) are exactly symmetric, within the uncertainties on the parameters. Note that, in general, this is not always the case: for simulations in which the system is driven to deviate strongly JHEP07(2020)233 from equilibrium, the distribution of ∆S obtained in "direct" transformations is not equal to the distribution of −∆S measured in "reverse" transformations -on the contrary, the two distributions are expected to cross each other at one point ∆S = ∆Γ: in the setup that we are considering, this is nothing but the statement of Crooks' theorem [78], from which eq. (2.1) can be immediately derived.
Coming back to the discussion of the results listed in table 17, we also note that all of the distributions analyzed there are characterized by values of the skewness and of the excess kurtosis compatible with zero, within their uncertainties. For example, for the data set reported in the first cell of the table, corresponding to "direct" transformations with ∆η = 0.015, we foundμ 3 0.02, with an uncertainty (estimated by jackknife binning) of a few units, and κ − 3 0.03, with an even larger uncertainty, and similar results persist for the other data sets. On the one hand, the fact that both the skewness and the excess kurtosis vanish, confirms that in the (near-equilibrium) regime probed by those simulations, the Euclidean-action distributions measured along the trajectories are approximated very well by Gaußians (as suggested by figure 1). Note, also, that the large uncertainties oñ µ 3 and on κ − 3 are not due to limited statistics (the number of independent trajectories analyzed for each data set reported in table 17 is of the order of 4 · 10 4 ), but to the very limited amount of data that deviate significantly from normal distributions; it is also worth remembering that bothμ 3 and κ − 3 are suitably normalized parameters, whose definitions encode fine cancellations among terms related to different moments of the distribution they characterize.
We conclude that the data in table 17 support the expectation that the distributions of values of ∆S can be modeled by normal distributions, whose mean values scale linearly with ∆η (see also figure 2). The width of these distributions tends to zero when the simulations approach the equilibrium limit, which happens when n qq is large for a given ∆η: in that limit, the distributions of ∆S from direct transformations become symmetric with respect to those from reverse transformations, and tend to sharp distributions centered exactly at ∆Γ.
For completeness, in table 18 we report the results of the cumulant analysis for a set of simulations of the SU(3) gauge theory. In this case, the data were obtained from nonequilibrium simulations with ∆η = 0.0001 and n qq = 1000 fixed. As before, the results show clearly that the distributions of ∆S are very sharply peaked, that the results obtained from direct and reverse transformations are fully consistent with each other, and the parameters describing deviations from a Gaußian shape are always consistent with zero.
Finally, it is also interesting to note that, for a fixed physical size of the system, the shape of the normalized distribution of ∆S values depends only very mildly on the linear size of the system in units of the lattice spacing: this is shown in figure 9, where the main plot displays the results that we obtained for the SU(3) theory from the set of simulations summarized in the first block of table 9. The plot refers to the Euclidean-action difference obtained in direct non-equilibrium transformations, on lattices corresponding to g 2 (L) 1.247, i.e. for a fixed physical size L, with ∆η = 0.0001, n qq = 1000, and for L/a ranging from 5 to 8. The four distributions show a remarkable collapse to a common curve (up to very small deviations in regions far from the peak, which however are not JHEP07(2020)233  significant within the statistical uncertainties, and a corresponding slight increase of s 2 with L/a), despite a nearly sevenfold increase in the number of degrees of freedom (N 2 − 1 for each of the link matrices not fixed by the Dirichlet boundary conditions discussed in section 2). Similar conclusions are obtained for the SU(2) theory: in this case, the inset in figure 9 shows the distributions obtained from direct non-equilibrium transformations, again with ∆η = 0.0001 and n qq = 1000, on lattices of approximately fixed physical size L at which g 2 (L) 2.04, for L/a values from 5 to 10 (corresponding to an increase in the number of degrees of freedom by a factor larger than 17), that are listed in the first block of table 2. Also in this case, the distributions obtained numerically, from our non-equilibrium simulations, exhibit a remarkable collapse to the same curve (up to slight deviations in regions far from the maximum, where the observed distributions are, however, smaller by more than three orders of magnitude with respect to the region close to the peak, and not significant within the uncertainties). In addition, we note that the differences among the various distributions in the maximum region are perhaps slightly more visible than in the SU(3) case. This, however, is likely to be at least partially due to the slightly larger relative difference in the values of g 2 (L) for this set of SU(2) simulations: comparing the results shown in the first block of table 2 and those in the first block of table 9, we observe that the relative variations in the g 2 (L) values are always well below 1%, but are slightly larger for the SU(2) theory than for the SU(3) theory.  These features can be compared and contrasted with those characterizing the distribution of values that are obtained by directly computing ∆Γ/∆η in ordinary, equilibrium Monte Carlo simulations. In that case, there is no departure from equilibrium at all, and hence no parameter analogous to our n qq ; accordingly, the distribution of values of ∆Γ/∆η has its own, fixed, width and shape, which has to be efficiently sampled by the simulation algorithm. In many cases, this task can be carried out without any challenges by ordinary Monte Carlo algorithms. For the SU(2) theory, in 1992 the authors of ref. [31] were able to achieve the significant level of precision of the results reported in that work using an amount of CPU time on CRAY YMP processors ranging from approximately 35 hours for the lattices of linear size in units of the lattice spacing L/a ≤ 14, to about 140 hours for the lattices with L/a = 20. For comparison, for the same theory the amount of CPU time that in this work we used to generate 3347 trajectories on a lattice with L/a = 20 and β = 3.7425 on a single core of CINECA machines equipped with Intel Xeon Phi 7250 CPU (Knights Landing) processors at 1.40 GHz is approximately 100 days. While this comparison may appear very humbling for our non-equilibrium algorithm, it is important to remark that these trajectories were produced with a large n qq value (n qq = 1000), which, JHEP07(2020)233 as can be appreciated from our final results, leads to enhanced precision for the value of the coupling that we could extract from them.
The efficiency of the non-equilibrium algorithm over standard equilibrium Monte Carlo, however, becomes particularly manifest in cases when the latter has to sample long tails in the distribution of ∆Γ/∆η, which lead to large autocorrelation times (as discussed, for example, in ref. [32] for simulations of the SU(3) theory at g 2 3.48). In a conventional equilibrium Monte Carlo algorithm, this problem should necessarily be addressed by applying reweighting methods, or some other similar technique. It is under these circumstances that our non-equilibrium algorithm proves particularly competitive in terms of CPU costs, since it allows one to bypass the computational overhead that is typically associated with reweighting and other analogous techniques. More precisely, the fact that our algorithm drives the field configurations to evolve out of equilibrium allows them to efficiently probe regions of the phase space of the system that would be exponentially difficult to access by standard reweighting methods.
In other words, the non-equilibrium algorithm discussed in the present work can be considered as one that simultaneously generalizes both standard Monte Carlo simulations and reweighting algorithms (recovering them in two particular limits, as discussed above). For quantities that can be efficiently estimated by conventional simulation algorithms, our code fares no better than them: although its design relies on a different approach, the physical observable it computes is, in fact, very similar to the one that is directly accessed by standard algorithms used for the determination of the Schrödinger-functional coupling, and the drawback of discretizing the trajectories into a sufficiently large number of steps is offset by the smaller fluctuations affecting the Euclidean-action differences measured along them. Under conditions where the dynamics of the theory is such, that conventional Monte Carlo algorithms are hampered by long autocorrelation times, difficulties in sampling the configuration space, or ensemble-overlap problems, however, non-equilibrium simulations can provide a very efficient alternative.
Finally, it is worth remarking that the conclusions we discussed above should hold regardless of the precise functional form of the time-dependent variation protocol for the parameters of the system that are modified during the non-equilibrium transformations (η in this case). While we restricted our analysis only to one class of such protocols, i.e. to sequences of "quantum quenches" of the same amplitude in Monte Carlo time, it is easy to guess that alternative choices for η(t) can have a significant impact on the computational efficiency of this type of non-equilibrium simulations. Interesting examples include a protocol η 1 (t) that remains constant to η(t in ) until the last step, when η is suddenly driven to η(t fin ), and another, η 2 (t), in which, conversely, η is immediately switched to η(t fin ) at the beginning of the trajectory, but then remains constant for the rest of the trajectory. The former would not drive the field configurations out of equilibrium at all during the whole trajectory except at the last step, and it would then determine ∆Γ simply by reweighting the last "measured" configuration to η = η(t fin ). By contrast, a protocol like η 2 (t) would closely mimic one sudden, large, initial "quench", which drives the field configurations violently out of equilibrium at the beginning of each trajectory, and then allows them to (tend to) relax towards the equilibrium state corresponding to η = η(t fin ). More interesting ex-JHEP07(2020)233 amples could include, for instance, protocols η(t) in which η is neither constant nor linear: while their investigation is beyond the scope of the present work, they may be optimized to improve the computational efficiency of the non-equilibrium algorithm, particularly under dynamical conditions in which it is competitive with respect to conventional equilibrium Monte Carlo simulations.

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ödingerfunctional scheme [30].
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 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 (continuumextrapolated) value of the coupling obtained on lattices of the same physical extent, but different lattice spacing. The application of this technique to study the non-perturbative renormalization in lattice QCD was pioneered in ref. [79], where the renormalization of the axial current, the running coupling in the chiral limit, and the momentum-scale evolution of the renormalized axial density were discussed, as well as the non-perturbative determination of the coefficients of improvement terms to reduce lattice artifacts. Since then, the formalism has been successfully applied to study the renormalized gauge coupling and a number of other physical quantities in QCD [80][81][82][83][84][85][86][87][88] and, more recently, has also been used to investigate the dynamics of other strongly coupled non-supersymmetric gauge theories with dynamical fermionic fields [89][90][91][92][93][94][95][96][97].
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 [11,12] 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 JHEP07(2020)233 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 Monte Carlo calculations. In particular, our approach closely "mimics" the non-trivial dynamics 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 [98].
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 [31,32]. While we presented results for purely bosonic theories, the generalization of this calculation to include dynamical fermions poses no additional conceptual 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 [34][35][36][37][38].
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 [99][100][101][102]; in particular, refs. [99,102] focus on the calculation of the entanglement entropy from Jarzynski's equality, whereas refs. [101,102] 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.