The Wilson flow and the finite temperature phase transition

We consider the determination of the finite temperature phase transition in the Yang-Mills SU(3) gauge theory. We compute the difference of the spatial and temporal energy density at a physical Wilson flow time. This difference is zero in the confined phase and becomes non zero in the deconfined phase. We locate the phase transition by using a new technique based on an exponential smoothing spline. This method is an alternative to the determination of the phase transition based on the Polyakov loop susceptibility and can also be used with dynamical fermions.


Introduction
The pure SU(3) gauge theory at finite temperature exhibits a first order phase transition separating the confined phase at low temperature from the deconfined phase at high temperature. The standard way to locate the critical temperature by Monte Carlo simulations is based on the computation of the susceptibility of the Polyakov loop. The position of the peak of the susceptibility determines the temperature and the scaling of the peak value with the spatial volume identifies the order of the transition.
In this work we study an alternative method to locate the phase transition based on the energy density computed using the Wilson flow. We separate the spatial and temporal components of the energy density. Their difference is zero in the confined phase and becomes non zero in the deconfined phase. This behaviour allows to locate the phase transition. Technically we propose to use an exponential smoothing spline to fit the data of the energy difference. We give explicit formulae for the construction of such a spline. We define the critical temperature to correspond to the maximum slope of the spline. We verify by Monte Carlo simulations that the new method agrees with the results of the standard method based on the Polyakov loop susceptibility.
The motivation to look at an alternative method to locate the phase transition is that in QCD with dynamical fermions one usually uses the peak of the chiral susceptibility to locate the phase transition or cross over (depending on the quark mass), which is a computationally expensive quantity. The method based on the energy difference which we describe in this paper can be applied also in presence of dynamical fermions.

JHEP10(2016)061
2 Finite temperature phase transition The finite temperature can be investigated on lattices where the time extension N t is much smaller than the spatial one N s . This means, it should at least hold N s ≥ 4·N t as described in [1]. The temperature T is related to the time extension by the equation with lattice spacing a(β). This implies for a given N t , there is a critical coupling β c such that the critical temperature can be computed via T c (β c ). Finally, the critical temperature can be computed in physical units, for example, as with c = 197.3 MeV fm and r 0 = 0.49 fm [2]. The formula for the determination of the scale r 0 /a(β) is for 5.7 ≤ β ≤ 6.92 and explained in [3].

Polyakov loop susceptibility
Usually, the critical coupling β c is determined using the Polyakov loop susceptibility. The lattice average P of the Polyakov loop is computed as which is the mean of the product of links U ( x, n t ) in time direction n t for each 3d-spatial vector x. As the lattice average P is a complex number, the Polyakov loop susceptibility χ p is the variance of its absolute value |P |, For a given lattice of size (N t , N 3 s ), the coupling constant β has to be chosen and through Monte Carlo simulation an ensemble of gauge configurations is generated. The expectation value in the ensemble of the Polyakov loop susceptibility χ p (β) can be computed. This has to be done for several values of β around the critical coupling β c such that the curve (β, χ p ) k for k = 1, . . . , n β contains a peak around β c . The desired value β c is found by a fit of the data, for example, with a smoothing spline s(β). The critical coupling β c is computed as the value of β at which the fit s(β) reaches its maximum value. Fixing N t and varying N s leads to several values β c (χ p , N t , N s ) of the critical coupling. Figure 1 shows the results of the simulations which we will describe in section 4. At the end, a linear extrapolation in 1/N 3 s leads to the desired critical coupling β c (χ p , N t , ∞) for infinite volume.

The Wilson flow
The Wilson flow [4][5][6] V (t) is a flow of lattice gauge fields and it is defined as the solution of the ordinary Lie group differential equatioṅ with link variables V x,µ (t) being elements of the special unitary Lie group SU(3) and a function Z x,µ (V (t)) which takes values in the Lie algebra su(3). Particularly, the function Z x,µ is the Lie derivative Z x,µ V (t) = −∂/∂V x,µ (t)S W of the Wilson action such that it depends not only on the link itself but also on its staples. S W is the sum over all oriented plaquettes p and V p (t) is the product of link variables around one of these plaquettes p, see [7]. The Wilson flow (2.3) can be integrated numerically with initial values V (t 0 ) taken from an ensemble of configurations generated for a particular value of the coupling constant β. The numerical integration is performed up to a certain flow time and during this computation gauge invariant observables of interest are measured. Thereby, we focus on the energy density with lattice spacing a and symmetric field strength tensor G µν as described in [7].

Difference in the Wilson energy
We investigate the mean energy density E described in equation (2.4) on a four-dimensional lattice with fixed temporal lattice size N t = 8. Particularly, we focus on the mean spatial and temporal part of the energy E ss and E st and its difference ∆E := E ss − E st . The splitting in the temporal and spatial part is done as follows numbering the dimensions 0, 1, 2, 3 with 0 being the time dimension: the spatial planes are (1, 2), (1, 3), and (2, 3) and the temporal ones (0, k) with k = 1, 2, 3. Here, the exact formulae for the spatial part of the energy E ss and the temporal one E st are The critical coupling β c divides the confined phase and the deconfined phase. In the confined phase, the values of β are smaller than β c , in the deconfined one larger. It is known from [1] that the critical coupling for lattices with N t = 8 is at β c =6.0625. For a first test, we simulated lattices with β = 6.03 < β c and β = 6.07 > β c and computed E ss , E st and its difference ∆E as a function of the flow time. It can be seen in figure 2 that in the confined phase (blue) with values of β smaller than β c there is no difference in the spatial and the temporal part of the energy. On the other hand, in the deconfined phase (red) with values of β larger than β c , there is a difference in both parts of the energy. This implies, the spatial and temporal mean energy densities E ss and E st coincide in the confined phase and differ in the deconfined one. In the following paragraphs, we fix a certain flow time √ t/r 0 in units of the scale r 0 such that ∆E means

JHEP10(2016)061
For a fixed flow time √ t/r 0 , it is shown in figure 2 (in this case √ t/r 0 = 0.4) that the energy difference ∆E is approximately zero in the confined phase and grows suddenly in the deconfined one. This advises that the critical coupling can be found using the energy difference at a certain flow time.

The energy difference method
We developed a new method 1 for the detection of the critical temperature and called it the energy difference method. Therefore, we consider the Wilson flow and focus on the difference in the spatial and temporal energy density at a certain flow time √ t/r 0 . The critical coupling β c can be computed by a fit through an exponential smoothing spline developed in [10]. The idea is to combine smoothing splines, which approximate the data in a spline setting, with an exponential spline, such that undesired oscillations not given in the data are avoided (see section 3 for more details). For the detection of the critical coupling, we proceed as follows similarly to the standard Polyakov loop susceptibility approach: first, we fix a lattice size (N t , N s ) and the values of β k for k = 1, . . . , n β . Then, we need the results of a HMC or heat bath simulation for these values as input for the Wilson flow. Moreover, the Wilson flow has to be computed up to a certain simulation time which has to be determined such that the statistical errors are minimized. We fixed the time √ t/r 0 = 0.15 as described in paragraph 4.2. Additionally, the values for the spatial and temporal energy density have to be measured for the flow time √ t/r 0 = 0.15. After the simulation, an exponential smoothing spline s(β) is determined to fit the data (β, ∆E ). The critical coupling β c (∆E, N t , N s ) for the specific lattice size (N t , N 3 s ) is determined as value of β where the steepest gradient of the spline s(β) occurs. This has to be repeated for a few spatial lattice sizes N s such that the values β c (∆E, N t , N s ) can be extrapolated towards an infinite space dimension N s = ∞ in order to compute β c (∆E, N t , ∞) for the finite temperature phase transition in infinite volume. Figure 3 shows the results of the simulations which we will describe in section 4.
The energy difference ∆E is equivalent (up to discretization effects) to the difference of the temporal and spatial plaquettes, cf. [7], and it corresponds to the sum of the pressure and the energy density, see e.g. [9]. It is expected to develop a discontinuity in the thermodynamic limit 2 since the energy density is discontinuos in the Yang-Mills SU(3) theory, see e.g. [1]. Presently our data shown in figure 3 do not allow to verify this expectation.

Exponential smoothing splines
The class of exponential smoothing splines approximates data with uncertainties avoiding undesired oscillations. It couples the ideas of smoothing splines [11] with exponential splines [12,13]. Here, we explain the idea of exponential smoothing splines and give the necessary formulae for its implementation. The mathematical essentials can be found in [10] based on [14] .

Idea
We start with data (x i , y i ) with errors w i , i = 1, . . . , n. The data should be approximated by a spline s(x) within the region of their errors and, moreover, no artificial oscillations should be added. The spline can be found by minimizing the energy function Here, Λ(x) avoid undesired oscillations, the weights w i are correlated to the errors of y i and the smoothing parameter S defines how much the approximated values f (x i ) may differ from y i . There are three limiting cases which coincide with already known kinds of splines: • λ = 0, S = 0 describes the well-known cubic spline, • λ = 0, S = 0 leads to the smoothing spline which approximates the data; it may lead to oscillations not given in the data, • and λ = 0, S = 0 results in the exponential spline (also known as spline under tension) interpolating the data exactly avoiding oscillations not given in the data.

The shape of the exponential smoothing spline
Finally, the exponential smoothing spline is parameterized as ) and coefficients s i , s i+1 , d i , d i+1 , µ i and λ i for all intervals i = 1, . . . , n − 1.
For convenience, the variables h i , t and µ i are used. They are specified as and q i := 1 h i for i = 1, . . . , n − 1, j = 1, . . . , n − 2 and k = 2, . . . , n − 2 such that the matrices read Thus, U and D are of size n × n, T of size (n − 2) × (n − 2) and and Q of size n × (n − 2).

JHEP10(2016)061
These formulae still depend on the unknowns λ i and µ i = λ i · h i . According to [13], the tension parameters λ i can be chosen as uniformly distributed random values in the interval [4h i , 15h i ].
Finally, the Lagrange parameter p is still unknown and can be computed by This can be done, for example, by a Newton iteration or interval nesting. Thereby, it must be taken into account that the starting value p (0) = 0 should be avoided since it would lead to vectors s = 0, d = 0 and therefore to a spline s(x) = 0, see [10].

Monte Carlo results
We simulate the Wilson plaquette gauge action [15] for gauge group SU(3) using the Hybrid Monte Carlo algorithm [16]. We let HMC simulations run for varying β (6.03 ≤ β ≤ 6.10) and the lattice sizes N t × N 3 s with N t = 8 and N s = {32, 40, 48}. Taking the gauge configurations of the HMC simulations as initial values for the Wilson flow, we computed the Wilson flow up to √ t/r 0 = 0.15 and measured its spatial and temporal energy density. Statistical errors and autocorrelation times are computed with the method of [17]. Besides, we computed the Polyakov loop susceptibility along the Wilson flow for comparison reasons.

Determination of the critical temperature
For the determination of the critical temperature, we computed the mean energy difference ∆E for all β at the flow time √ t/r 0 = 0.15 and for each lattice separately. Then, the data (β, ∆E ) are fitted by exponential smoothing splines s(β), see figure 3. Here, the smoothing parameter S is set to S = 2 and the weights w i used in the matrix D are set to the statistical errors of ∆E i . Note that a small variation of S does not change the result; if S is chosen too large, the exponential smoothing spline does not fit each value of ∆E in the region of its error δ ∆E . The maximum of the slope s (β) of this spline leads to the critical coupling β c (∆E, N t , N s ) which is the value of β at which s (β) = 0 holds. The errors δβ c (∆E, N t , N s ) are computed by a small variation of the single values of ∆E using the Gaussian error propagation law: The partial derivatives ∂β c /∂ ∆E k are computed using the symmetric difference quotient with perturbed value ∆E k ; we use 10 percent of the error of ∆E k . At the end, a linear extrapolation in 1/N 3 s leads to the final value β c (∆E, N t = 8, N s = ∞) = 6.0601 (11)  value 6.0624 (12) computed in [18] and with the newest result 6.06239(38) of ref. [19]. Using equation (2.2) this result translates to T c r 0 = 0.7426 (14)(37), where the first error is the uncertainty in β c and the second one is a 0.5% error for r 0 /a computed from equation (2.2). For comparison reasons, we determined the critical coupling from our data with the standard Polyakov loop susceptibility approach described in paragraph 2.1 and shown in figure 1. Also this method leads to a critical coupling β c (χ p , N t = 8, ∞) = 6.0602 (7) close to the value given in the literature [18]. Both methods lead to very consistent values of the critical coupling and also the size of the errors are of the same magnitude.

Statistical errors at different flow times
We investigated the statistical errors of both methods at different flow times. For the Polyakov loop susceptibility, the relative statistical errors do not change during the flow. The Wilson flow has just an effect on its absolute value, see figure 5. Therefore, it would be no advantage to combine the Wilson flow and the Polyakov loop susceptibility. On the other hand, the statistical errors for the energy difference method have a minimum at √ t/r 0 = 0.15. We checked that this flow time is large enough such that cut-off effects (measured by comparing different definitions of the energy density) are small and therefore the best choice for the computation of the critical temperature.
The statistical errors already include autocorrelation effects shown in figure 6 which are quite large close to the critical coupling. For this reason, long simulations have to be run for these values of the coupling constant β such that the number of independent  configurations is large enough. Surprisingly, the number of independent configurations for the energy difference at flow time √ t/r 0 = 0.15 is for almost all values of β (except at β c ) larger than the one of the Polyakov loop susceptibility at √ t/r 0 = 0, mentioned in table 1. This means, the number of configurations needed for the energy difference method is smaller than for standard approach.

Conclusion and outlook
For pure gauge theory, it is possible to detect the critical coupling β c via the energy difference of the Wilson flow. Our results agree with the standard method as well as with the reference value given in [1,18] and [19]. For the energy difference method, the Wilson flow has to be computed in addition to the HMC as it leads to a reduction of the statistical error. This is a disadvantage compared to the standard approach but just a small one, since the energy difference can be evaluated at small flow times. For most values of β the statistical errors of the energy difference method are smaller than for the standard approach such that the simulation has to be run for less configurations to reach the same amount of independent configurations. Moreover, our results are promising for simulations with fermions. Our method can be used as alternative to the expensive computation of the chiral susceptibility, which is usually taken in this case. Our technique based on the exponential smoothing spline is also suitable for other applications, as it allows to fit a smooth function through a data set and compute its derivatives. Additionally, the integration of the Wilson flow can be improved by using better integrators. For example, adaptive step size methods developed in [20] and [21] reduce at the same time the computational cost and are able to control the statistical errors.