The equation of state with non-equilibrium methods

Jarzynski's equality provides an elegant and powerful tool to directly compute differences in free energy in Monte Carlo simulations and it can be readily extended to lattice gauge theories to compute a large set of physically interesting observables. In this talk we present a novel technique to determine the thermodynamics of strongly-interacting matter based on this relation, which allows for a direct and efficient determination of the pressure using out-of-equilibrium Monte Carlo simulations on the lattice. We present results for the equation of state of the $\mathrm{SU}(3)$ Yang-Mills theory in the confined and deconfined phases. Finally, we briefly discuss the generalization of this method for theories with fermions, with particular focus on the equation of state of QCD.


Introduction
The determination of the equation of state of strongly-interacting matter represents a crucial endeavour in theoretical physics, with many applications in different fields such as nuclear physics and cosmology. It serves as an input for the analysis of thermal systems such as those created in heavy-ion collision experiments or for the study of the early phases of the Universe itself. From the theoretical side, the most significant and reliable contribution to this effort comes undoubtedly from the lattice formulation of QCD, which provides a tool for first-principles numerical predictions with ever-increasing precision and accuracy. In the last few years there have been major advancements in the computation of equilibrium thermodynamics for full QCD with 2 + 1 (or more) dynamical quark flavors; still, such calculations require an impressive numerical effort and many systematic effects have to be taken into account. Thus, recently there has been renovated interest in the study of new ways of computing the equation of state in addition to standard techniques such as the integral method [1]: among the latest advancements, we mention studies in a moving reference frame [2] and with the gradient flow [3].
The purpose of this paper is to present a novel method for the calculation of the pressure in lattice gauge theories exploiting a well-known result by C. Jarzynski in out-of-equilibrium statistical mechanics. The so-called Jarzynski's equality [4,5] relates the difference in free energy ∆F between two equilibrium states with the average of the exponential of the work done on the system of interest during a transformation between the two states: crucially, such transformations in general will not be performed with the system in thermodynamic equilibrium. Even more importantly, the average of the exponential of the work must be taken on an ensemble of realizations of this transformation.
Speaker, e-mail: anada@to.infn.it In Ref. [6] Jarzynski's equality was succesfully tested in the context of lattice gauge theories in two different benchmark studies: the first concerning the free energy associated to an interface in the Z 2 gauge model, and the second focusing on the pressure in the SU(2) gauge theory on a small range of temperatures.
In this work, which is a natural prosecution of the work done in Ref. [6], we present a preliminary study of the equation of state of the SU(3) Yang-Mills theory obtained with non-equilibrium methods based on Jarzynski's equality: the theory without quarks represents a perfect testing ground for new techniques, since one can avoid the complications related to dynamical fermionic fields. At the same time, the seminal work on SU(3) thermodynamics of Ref. [7] has been improved in recent years by high-precision determinations obtained with different methods [2,8], that give us the possibility to test the reliability of our technique.

Jarzynski's equality
In this section we will state Jarzynski's equality precisely, and analyse in detail its practical implementation in Monte Carlo simulations. In order to discuss the non-equilibrium work relation, we start from the second law of thermodynamics in the form of the Clausius inequality where S is the entropy, T is the temperature, and Q is the heat exchanged with the environment during a transformation between macrostates A and B; for an isothermal transformation it can be rewritten as using the first law (∆E = Q + W) and the definition of free energy F = E − S T , E being the internal energy. Let us now consider a system whose dynamics is described by a Hamiltonian H λ that depends explicitly on a certain parameter λ, e.g. the coupling of the model; the corresponding partition function Z will be Z(T, λ) = dΓe −H(Γ,λ)/T where Γ indicates a microstate of the system and k B = 1. In this framework, we can think the transformation A → B to be driven by a change in this parameters, i.e. λ A → λ B . Moreover, we know that for a microscopic system Eq. (2) is valid only statistically, i.e.
where the ... from now on will denote an average over all possible realizations of λ A → λ B transformation, during which the total work W spent to perform the switch in λ is measured. We can now state the non-equilibrium work relation by C. Jarzynski [4,5] exp that puts in relation the average of the exponential of the work performed in an isothermal transformation λ A → λ B with the difference in free energy between the initial and final states or, equivalently, the corresponding ratio of partition functions Z. Using Jensen's inequality, i.e. e x ≥ e x , valid for a real variable x, it is easy to show that Eq. (4) is a generalization for microscopic systems of the second law of thermodynamics.

The non-equilibrium work relation for Monte Carlo simulations
Jarzynski's equality has been derived also in the context of stochastic processes (see for example [5,9]) and in particular for Markov chains: thus, the implementation for Monte Carlo simulations is rather straightforward. However, before using Eq. (4) it is crucial to understand what precisely W is and how in practice the non-equilibrium transformation is performed. Firstly, the transformation has to be discretized into N intervals, so that at each intermediate step the λ parameter changes: where λ 0 corresponds to the initial macrostate previously denoted as A and λ N to B; the nonequilibrium relation does not depend on the specific protocol used to switch λ. We will also have the corresponding set of intermediate configurations [φ n ] of the system where, crucially, [φ 0 ] must be a thermalized configuration. The work W is defined quite naturally as the sum of the difference in the Hamiltonian at each step of the Markov process: Let us analyze how the entire non-equilibrium transformation is implemented in practice during a Monte Carlo simulation. These are the steps that must be followed: 1. the non-equilibrium work relation requires the system to be at equilibrium at the beginning of each trajectory; 2. we switch the parameters from λ 0 to λ 1 , following the chosen protocol; 3. we compute the work done on the system to perform this first change of λ, simply taking the difference of the Hamiltonian note that the Hamiltonians are evaluated using the same configuration but different values of λ; 4. we update the system with the algorithm of choice to the new configuration [φ 1 ] keeping λ fixed to λ 1 5. we repeat steps 2, 3 and 4 until the transformation is completed. At each step n, the parameters are changed following the given protocol λ n → λ n+1 , the work performed on the system is computed as and the system is then updated using the new parameter 6. at the end of each trajectory, the total work defined in Eq. (5) is computed; 7. a new equilibrium configuration [φ 0 ] is generated by thermalizing the system again with λ 0 , and a new trajectory can begin.
It is extremely important to stress that one has to perform several independent realizations of the transformation so that the exponential average of Eq. (4) provides reliable results. The interplay between the number of such realizations, denoted as n R , and the number N of intervals in λ is crucial to improve the efficiency of this technique. We conclude this section by noting that the relation can be extended to non-isothermal transformations (see, for example, Ref. [10]).

The equation of state with Jarzynski's equality
Following the work of Ref. [6], in this section we will review how to compute the pressure using non-equilibrium transformations in finite-temperature lattice simulations. We start by considering a model with a given partition function Z(T ) and a free energy density f = −T (ln Z)/V, defined on an hypercubic lattice Λ of sizes aN t × (aN s ) 3 , with N t and N s representing the number of lattice sites in the temporal and spatial directions. The spatial volume corresponds to V = (aN s ) 3 , while the temperature is defined as usual as T = 1/(aN t ). Our primary observable is the pressure p, that in the thermodynamic limit can be written as Other relevant thermodynamical quantities are the energy density and the trace of the energy-momentum tensor ∆ = − 3p, which can be conveniently expressed as The dimensionless ratio p(T )/T 4 can be written as and if we compute differences in p/T 4 between two temperatures T and T 0 , then we have and it is the Z(T )/Z(T 0 ) ratio that can be computed directly with Jarzynski's equality using nonequilibrium transformations in which the temperature T is varied. In practice, the temperature is changed throughout each trajectory by changing the lattice spacing a, i.e. by tuning the inverse coupling β to the desired values. Indeed, Eq. (4) becomes where ∆S is the total change in the Euclidean action and ... indicates the average on a set of n R nonequilibrium trajectories performed as described in section 2.1. The protocol to change the parameter (in this case the inverse coupling β) in the transformation is chosen to be linear in the index n of the intermediate steps, so that of course, in this notation T 0 ≡ 1/(a(β 0 )N t ) and T ≡ 1/(a(β N )N t ). ∆S is the quantity computed for each trajectory and it is the equivalent of the "work" defined in Eq. (5), in which the action S has taken the place of H/T . Before being able to compute the physical value of the pressure, we need to take care of the quartic divergence in a: a simple way to do so is to compute the same quantity of Eq. (11) but at T = 0, i.e. on a symmetric lattice Λ of hypervolume (a N) 4 . Then we can subtract it from the result of the finite-temperature simulation so that where the exponent γ = N 3 s × N t / N 4 is necessary since the lattice hypervolumes at T = 0 and T 0 are in general different.

Numerical results for the SU(3) pure-glue theory
For the Euclidean Yang-Mills action of the lattice theory we chose the standard Wilson action [11], Re Tr U µν (x) where g denotes the (bare) lattice coupling and is the plaquette. The partition function of our theory is where dU µ (x) denotes the Haar measure of the SU(N) variable at the site x and direction µ.
The simulations were performed on lattices with N t = 6, 7, 8, 10, while keeping the aspect ratio N s /N t > 12 and in most cases around 16, in order to avoid finite-size effects. The scale is controlled using the values of the Sommer scale r 0 /a determined in [12]; moreover we used for the critical temperature the value T c r 0 = 0.7457(45) computed in [13]. The finite lattice spacing results were first interpolated with spline functions and then the continuum extrapolation was performed with a linear fit in (1/N t ) 2 . Preliminary results for the pressure are illustrated in Fig. 1.
As it can be seen, the data for the pressure obtained using the technique based on Jarzynski's equality are in good agreement with previous high-precision determinations by the Wuppertal-Budapest collaboration [8] and more recently by L. Giusti and M. Pepe [2]. We remark however that these two last computations showed a small but clearly visible discrepancy in the deconfining phase, in particular in the region between T c and 1.5T c ; this disagreement becomes gradually smaller as T increases and disappears for T > 3T c . We will not attempt a discussion of the possible origin of this problem here: we limit ourselves to note that our preliminary data for p/T 4 seem to agree well with those of  Ref. [2] in the aforementioned region. In order to investigate this issue more in detail, we also present preliminary data for the trace anomaly ∆ in Fig. 2 and the energy density in Fig. 3.
First, we note that in order to obtain ∆ it is necessary to numerically derive the pressure p with respect to T : this is to be contrasted with the integral method [1] and computations in a moving frame, in which a secondary observable such as the pressure requires numerical integration from the quantity that is directly measured on the lattice, i.e. the trace anomaly ∆ or the entropy density s = ( + p)/T . In order to do so we performed first a Padé interpolation and then the derivation of the resulting fit: preliminary results in Fig. 2 and Fig. 3 confirm that our method seems to favor the results of Ref. [2]. We report that the most delicate temperature region to probe is the one immediately above T c , where a strong dependence on the lattice spacing was observed, leading to larger uncertaintes in the a → 0 extrapolation.
As previously discussed in Ref. [6], the exponential average of Eqs. (4) and (13) requires a sample of n R trajectories that is large enough, in order to converge to the correct result. The correct way to ensure this is to repeat the transformation in the reverse direction and to check if the results agree with the "direct" one; this has been performed in this work, and whenever the agreement was not satisfactory, the transformation was repeated with a new combination of N intermediate steps and n R realizations.
A possible issue related to this novel technique concerns the use of non-equilibrium transformations that cross the deconfinement phase transition, as in such cases Jarzynski's equality cannot be used. This fact is confirmed by the strong disagreement between results of trajectories going from the confined to the deconfined phase, and trajectories which perform the same transformation in the reverse direction. In order to avoid this issue, our non-equilibrium transformations never crossed T c .

Conclusions
In this work a new determination of the SU(3) equation of state has been performed, using a technique based on Jarzynski's equality: preliminary results obtained with out-of-equilibrium Monte Carlo simulations show good agreement with past determinations. The new method proved to be very efficient, since only the starting configuration has to be at equilibrium, but also highly reliable, as each transformation has to be in agreement with the same transformation performed in the opposite direction. Jarzynski's equality has a wide range of possible applications in lattice gauge theory, as the problem of computing differences in free energy is a very general and common one; however the most natural prosecution of this work is the application of this technique to the equation of state in full QCD. A particularly intriguing idea is to set up non-equilibrium transformations in which both the temperature T (via the inverse coupling β) and the bare masses of the fermions are changed simultaneously.