The equation of state with nonequilibrium methods
Abstract
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 stronglyinteracting matter based on this relation, which allows for a direct and efficient determination of the pressure using outofequilibrium Monte Carlo simulations on the lattice. We present results for the equation of state of the YangMills 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.
The equation of state with nonequilibrium methods
Alessandro Nada^{}^{}Speaker, email: anada@to.infn.it, Michele Caselle, and Marco Panero
1Department of Physics, University of Turin & INFN, Turin
Via Pietro Giuria 1, I10125 Turin, Italy
\@textsuperscript2ArnoldRegge Center, University of Turin
Via Pietro Giuria 1, I10125 Turin, Italy

Abstract.
1 Introduction
The determination of the equation of state of stronglyinteracting 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 heavyion 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 firstprinciples numerical predictions with everincreasing precision and accuracy. In the last few years there have been major advancements in the computation of equilibrium thermodynamics for full QCD with (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 wellknown result by C. Jarzynski in outofequilibrium statistical mechanics. The socalled Jarzynski’s equality [4, 5] relates the difference in free energy 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.
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 gauge model, and the second focusing on the pressure in the 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 YangMills theory obtained with nonequilibrium 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 thermodynamics of Ref. [7] has been improved in recent years by highprecision determinations obtained with different methods [8, 2], that give us the possibility to test the reliability of our technique.
2 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 nonequilibrium work relation, we start from the second law of thermodynamics in the form of the Clausius inequality
(1) 
where is the entropy, is the temperature, and is the heat exchanged with the environment during a transformation between macrostates and ; for an isothermal transformation it can be rewritten as
(2) 
using the first law () and the definition of free energy , being the internal energy. Let us now consider a system whose dynamics is described by a Hamiltonian that depends explicitly on a certain parameter , e.g. the coupling of the model; the corresponding partition function will be
where indicates a microstate of the system and . In this framework, we can think the transformation to be driven by a change in this parameters, i.e. . Moreover, we know that for a microscopic system Eq. (2) is valid only statistically, i.e.
(3) 
where the from now on will denote an average over all possible realizations of transformation, during which the total work spent to perform the switch in is measured. We can now state the nonequilibrium work relation by C. Jarzynski [4, 5]
(4) 
that puts in relation the average of the exponential of the work performed in an isothermal transformation with the difference in free energy between the initial and final states or, equivalently, the corresponding ratio of partition functions . Using Jensen’s inequality, i.e. , valid for a real variable , it is easy to show that Eq. (4) is a generalization for microscopic systems of the second law of thermodynamics.
2.1 The nonequilibrium 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 is and how in practice the nonequilibrium transformation is performed. Firstly, the transformation has to be discretized into intervals, so that at each intermediate step the parameter changes:
where corresponds to the initial macrostate previously denoted as and to ; the nonequilibrium relation does not depend on the specific protocol used to switch . We will also have the corresponding set of intermediate configurations of the system
where, crucially, must be a thermalized configuration. The work is defined quite naturally as the sum of the difference in the Hamiltonian at each step of the Markov process:
(5) 
Let us analyze how the entire nonequilibrium transformation is implemented in practice during a Monte Carlo simulation. These are the steps that must be followed:

the nonequilibrium work relation requires the system to be at equilibrium at the beginning of each trajectory;

we switch the parameters from to , following the chosen protocol;

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 ;

we update the system with the algorithm of choice to the new configuration keeping fixed to

at the end of each trajectory, the total work defined in Eq. (5) is computed;

a new equilibrium configuration is generated by thermalizing the system again with , 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 , and the number 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 nonisothermal transformations (see, for example, Ref. [10]).
3 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 nonequilibrium transformations in finitetemperature lattice simulations. We start by considering a model with a given partition function and a free energy density , defined on an hypercubic lattice of sizes , with and representing the number of lattice sites in the temporal and spatial directions. The spatial volume corresponds to , while the temperature is defined as usual as . Our primary observable is the pressure , that in the thermodynamic limit can be written as
(6) 
Other relevant thermodynamical quantities are the energy density
(7) 
and the trace of the energymomentum tensor , which can be conveniently expressed as
(8) 
The dimensionless ratio can be written as
(9) 
and if we compute differences in between two temperatures and , then we have
(10) 
and it is the ratio that can be computed directly with Jarzynski’s equality using nonequilibrium transformations in which the temperature is varied. In practice, the temperature is changed throughout each trajectory by changing the lattice spacing , i.e. by tuning the inverse coupling to the desired values. Indeed, Eq. (4) becomes
(11) 
where is the total change in the Euclidean action and indicates the average on a set of 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 of the intermediate steps, so that
(12) 
of course, in this notation and . is the quantity computed for each trajectory and it is the equivalent of the “work” defined in Eq. (5), in which the action has taken the place of . Before being able to compute the physical value of the pressure, we need to take care of the quartic divergence in : a simple way to do so is to compute the same quantity of Eq. (11) but at , i.e. on a symmetric lattice of hypervolume . Then we can subtract it from the result of the finitetemperature simulation so that
(13) 
where the exponent is necessary since the lattice hypervolumes at and are in general different.
4 Numerical results for the pureglue theory
For the Euclidean YangMills action of the lattice theory we chose the standard Wilson action [11],
(14) 
where denotes the (bare) lattice coupling and
(15) 
is the plaquette. The partition function of our theory is
(16) 
where denotes the Haar measure of the variable at the site and direction .
The simulations were performed on lattices with , while keeping the aspect ratio and in most cases around , in order to avoid finitesize effects. The scale is controlled using the values of the Sommer scale determined in [12]; moreover we used for the critical temperature the value 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 . 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 highprecision determinations by the WuppertalBudapest 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 and ; this disagreement becomes gradually smaller as increases and disappears for . We will not attempt a discussion of the possible origin of this problem here: we limit ourselves to note that our preliminary data for 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 with respect to : 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 . 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 , where a strong dependence on the lattice spacing was observed, leading to larger uncertaintes in the extrapolation.
As previously discussed in Ref. [6], the exponential average of Eqs. (4) and (13) requires a sample of 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 intermediate steps and realizations.
A possible issue related to this novel technique concerns the use of nonequilibrium 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 nonequilibrium transformations never crossed .
5 Conclusions
In this work a new determination of the equation of state has been performed, using a technique based on Jarzynski’s equality: preliminary results obtained with outofequilibrium 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 nonequilibrium transformations in which both the temperature (via the inverse coupling ) and the bare masses of the fermions are changed simultaneously.
Acknowledgements
The simulations were run on GALILEO and MARCONI supercomputers at CINECA. We thank M. Hasenbusch and R. Sommer for helpful comments and insightful discussions.
References
 [1] J. Engels, J. Fingberg, F. Karsch, D. Miller, M. Weber, Phys. Lett. B252, 625 (1990)
 [2] L. Giusti, M. Pepe, Phys. Lett. B769, 385 (2017), 1612.00265
 [3] M. Kitazawa, T. Iritani, M. Asakawa, T. Hatsuda, H. Suzuki, Phys. Rev. D94, 114512 (2016), 1610.07810
 [4] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)
 [5] C. Jarzynski, Phys. Rev. E56, 5018 (1997)
 [6] M. Caselle, G. Costagliola, A. Nada, M. Panero, A. Toniato, Phys. Rev. D94, 034503 (2016), 1604.05544
 [7] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lütgemeier, B. Petersson, Nucl. Phys. B469, 419 (1996), heplat/9602007
 [8] S. Borsányi, G. Endrődi, Z. Fodor, S. Katz, K.K. Szabó, JHEP 07, 056 (2012), 1204.6184
 [9] G.E. Crooks, Journal of Statistical Physics 90, 1481 (1998)
 [10] C. Chatelain, J. Stat. Mech. 0704, P04011 (2007), condmat/0702044
 [11] K.G. Wilson, Phys. Rev. D10, 2445 (1974)
 [12] S. Necco, R. Sommer, Nucl. Phys. B622, 328 (2002), heplat/0108008
 [13] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus, H. Ohno, Phys. Rev. D91, 096002 (2015), 1503.05652