The density of states method in Yang-Mills theories and first order phase transitions

Extensions of the standard model that lead to first-order phase transitions in the early universe can produce a stochastic background of gravitational waves, which may be accessible to future detectors. Thermodynamic observables at the transition, such as the latent heat, can be determined by lattice simulations, and then used to predict the expected signatures in a given theory. In lattice calculations, the emergence of metastabilities in proximity of the phase transition may make the precise determination of these observables quite challenging, and may lead to large uncontrolled numerical errors. In this contribution, we discuss as a prototype lattice calculation the first order deconfinement transition that arises in the strong SU(3) Yang-Mills sector. We adopt the novel logarithmic linear relaxation method, which can provide a determination of the density of states of the system with exponential error suppression. Thermodynamic observables can be reconstructed with a controlled error, providing a promising direction for accurate model predictions in the future.


Introduction
Proposals for physics beyond the standard model (BSM) based on strongly interacting non-Abelian gauge sectors can provide an origin for a variety of phenomenologically interesting scenarios, such as dark matter and composite Higgs models (see, e.g., Refs.[1][2][3][4][5] for examples based on the Sp(4) gauge group).Our current understanding of these proposals is not yet detailed and precise, due to their strong-coupling nature.A signature with great potential discovery (or exclusion) reach is the generation of stochastic gravitational wave backgrounds from phase transitions in the early universe.In fact, many strongly interacting models are expected to undergo a change in behaviour between an early universe deconfined phase of quark gluon plasma and the late universe confined phase.Models characterised by a firstorder phase transition lead to bubble nucleation which drives the generation of gravitational waves (see, e.g., Refs.[6][7][8], and [9] in these proceedings).The signature of the gravitational waves produced through bubble nucleation can in principle be predicted looking at thermodynamic observables measured at the phase transition.Of particular importance are the latent heat, which drives the expansion, and the surface tension, acting as a frictional term.
In principle, Lattice Field Theory is a prime tool to compute non perturbatively observables related to the production of relic stochastic gravitational waves in such strongly-coupled models.However, the physical metastabilities that arise near criticality limit the power and reach of commonly used lattice algorithms, and call for the development of new ones.In this contribution, we discuss a novel application of the linear logarithmic relaxation (LLR) algorithm [10], by analysing the thermodynamics of the phase transition in lattice simulations.This algorithm avoids the problems that affect standard Monte Carlo methods in a regions of metastability and hence enable us to calculate thermodynamic observables with robust error estimates.In the next section we discuss the general thermodynamics properties of non-Abelian lattice gauge theories at finite temperature, focusing in particular on the aforementioned algorithmic difficulties associated with metastability.In Sect. 3 the LLR method is reviewed.Then, in Sect.4, some preliminary results are discussed for the SU(3) pure gauge theory (see also Ref. [11] for a review and a recent high-precision calculation).Finally, Sect. 5 contains our summary.We remark that a similar approach is being undertaken for SU(4) (see Ref. [12] and the talk by F. Springer at this conference).

First order phase transitions and metastable dynamics
The Euclidean space-time is discretised onto an isotropic lattice with Ṽ/a 4 = N t × N 3 s sites, where a is the lattice spacing.Since in this work we interpret the lattice field theory as a statistical mechanics system in four dimensions, for simplicity we set a = 1.Periodic boundary conditions are imposed in all directions.The spatial lattice size, N s , is much larger than the temporal one, N t .For fixed N t , the temperature is set by changing the lattice spacing a non-perturbatively, via the lattice coupling β.
The Wilson action for SU(N c ) pure gauge theory is where U µν ( j) denotes a plaquette, which is the parallel transport of the lattice links U µ (i) ∈ SU(N c ) around the elementary square of the lattice in the µν plane, starting at lattice site j.
In standard importance sampling approaches, configurations are generated via Monte Carlo Markov chain methods, with a combination of heat bath and over-relaxation steps, with a probability distribution P β (E) for a value of the action S = E given by where the integral is taken over all SU(N c )-valued link variables U µ (i).Vacuum expectation values (VEV) of observables, ⟨O⟩, are estimated as (ensemble) averages over the sampled configurations.
When focusing on the thermodynamics at the phase transition, the observables of interest are the average plaquette, which we denote as u p , the specific heat (or plaquette susceptibility), the Polyakov loop, and the Polyakov loop susceptibility.The Polyakov loop is defined in terms of parallel transports of the link variables on paths of fixed spatial coordinates that wrap around the temporal direction.The Polyakov loop VEV [13] is the order parameter of the spontaneous breaking of the global centre symmetry associated with the deconfinement phase transition.It is given by where we have introduced the notation i = (n t , ⃗ n s ), with n t being the temporal and ⃗ n s the spatial coordinates of the lattice site i.
In SU(N c ) Yang-Mills gauge theory with N c > 2, the deconfinement phase transition is of first order and therefore exhibits metastable dynamics around the critical point.To accurately measure the thermodynamic quantities of interest for the system with a lattice calculation based on importance sampling, the configuration space must be sampled by an ergodic Markov chain.To ensure ergodicity, in the metastable region the Markov chain must allowed to tunnel between the two phases several times.But standard local Monte Carlo methods only consider small changes in each link, suppressing drastic changes.Consequently, these algorithms struggle to overcome the potential barrier separating the two phases.While for small lattices, and small number of colours N c (between three and six), this is not a problem, in the limit of large volumes (or large N c ) the potential barrier grows.More configurations are therefore required to ensure the system can explore all the configuration space, with the associated computational time growing exponentially with the volume, to compensate for the tunneling probability being exponentially suppressed with the volume.

Linear logarithmic relaxation method
The linear logarithmic relaxation (LLR) method allows, in principle, to avoid the problems in simulating the equilibrium distribution of the metastable dynamics near the critical region of first-order phase transitions.It uses Monte Carlo algorithms with the purpose of recovering microcanonical information, which is then used to reconstruct physical observables.
In the LLR approach to a general statistical system the energy range of interest is broken down into intervals, with energy The interval size δ E should be chosen small enough for a Taylor expansion of physical quantities in δ E to be a good approximation in each interval.The total action of Yang-Mills theories is used as an analogue for the energy1 , E = 6 Ṽ(1 − u p ).For each interval the density of states, is approximated numerically.The partition function of the statistical system with Hamiltonian S (or, equivalently, path integral of the field theory with action S ) is rewritten in terms of ρ(E), and becomes a simple one-dimensional integral over the allowed energy interval: Similarly, the VEV computed at coupling β of an observable O(E) that strictly depends on the energy becomes We approximate ρ(E) in a small interval of amplitude δ E around E n , by Taylor expanding its logarithm, ln(ρ(E)) ≈ a n (E − E n ) + c n , where c n is a constant.By rearranging this equation and using textbook statistical mechanics considerations, the Taylor coefficients a n are identified with the inverse microcanonical temperature at energy E n , as a n = (d ln ρ/dE)| E=E n = (ds/dE)| E=E n ≡ 1/t n , where s is the entropy.
The individual Taylor expansions, repeated over a sequence of intervals E 0 < ... < E n < E n+1 < ... < E N , for E n+1 = E n + δ E , can be combined by imposing continuity, to obtain, for the density of states in an interval This expression can be evaluated for any E in the range The overall factor ρ 0 can be fixed arbitrarily, as it drops in computations of ensemble averages.

Results
We can think of the coefficients a n in terms of the temperature that gives rise to an energy distribution that is symmetric around the centre of the interval E n .This implies ⟨⟨∆E⟩⟩ a n ≡ ⟨⟨E − E n ⟩⟩ a n = 0, where the double brackets denote the sampled mean over configurations restricted to the interval.To solve this equation, we follow Ref. [14] and use the Robbins-Monro (RM) approach [15].We choose an initial guess for a (0)  n ≈ a n .Then, using the equation a (m+1) n = a (m)  n − 12⟨⟨∆E⟩⟩ a (m) n /δ 2 E (m + 1), we iteratively improve our guess of the true value of a n .For a suitable initial value, after a transient, a (m)  n will oscillate around the root a n with an increasingly suppressed amplitude, and lim m→∞ a (m)  n = a n .An example of the calculation of a n for the SU(3) system considered here is shown in Fig. 1.We assess the size of the truncation error of the iteration by repeating the stochastic calculation 20 times, beginning with the same initial guess, and averaging the results.As m increases, all the repeats appear to converge towards the same final value.The values of a n have been calculated for intervals E n that cover the energy range relevant for the physical problem we are interested in.To avoid ergodicity problems associated with the restricted sampling, umbrella sampling was used when applying RM iterations, as described in Ref. [16].
The VEV of observables that strictly depend on the energy at a coupling β are computed by plugging the numerically determined values of a n into Eq.( 6), as follows An analogous expression for Z(β) is left implicit.This equation displays a sum over contributions from all intervals, but notice the suppression of contributions for which the associated inverse micro-canonical temperature is far from the coupling β we are calculating at.
Using Eq. ( 8) we have calculated the VEV of the average plaquette, ⟨u p ⟩ β , and the specific heat, C V ≡ ⟨u 2 p ⟩ β − ⟨u p ⟩ 2 β .These are shown in Fig. 2 plotted against the coupling value they are calculated at.To demonstrate the validity of this method, the plots include the values measured using standard importance sampling methods.
An important observable for gravitational wave physics is the latent heat, L h , of the transition.As shown, e.g., in Ref. [17], the key thermodynamic quantity that enters the calculation of L h is the average plaquette jump at criticality, ∆⟨u p ⟩ β c .For the definition of the critical point, we take the value of β at which both phases are equally probable.The corresponding energy distribution can be approximated with a double Gaussian in which the two peaks corresponding to the two different phases are of equal height [18].The critical point is found by taking the value of β at which the specific heat has a maximum, as an initial guess, and then using a bisection algorithm to refine the estimate until a predefined tolerated different height between the peaks is met.∆⟨u p ⟩ β c is then the energy difference between the position of the confined peak and the position of the deconfined peak at the determined value of β c divided by 6 Ṽ, as shown in Fig. 3.
The presence of metastability can be inferred directly from the non invertibility of a n (E n ), as shown in Fig. 3.The points at which a n = β c correspond to the locations of the extrema of the distribution.The thermodynamics of the system can be further analysed using the micro-canonical ensemble.In analogy with a statistical mechanics system, we can define the temperature t, entropy s, and the free energy of the micro-states as follows: Similarly to the case of a Van der Waals gas below criticality, we expect the free energy to display a swallow tail structure (see, e.g., Ref. [19]), that indicates a first order transition.
A subtlety emerges because the definition of the entropy depends on the constant ρ 0 in the density of states.ρ 0 contributes a term linear in the temperature to the free energy.As we are only interested in seeing the qualitative behavior of the free energy, we subtract a linear term Σt, where Σ is computed as the temperature average of the entropy within the energy range considered.The result is shown in Fig. 4.This clearly shows the telltale signs of a first order transition.Away from the critical region the free energy is single valued.In the critical region the free energy has three values at any given t.The highest such value is associated with the unstable vacuum, the middle one with the false vacuum and the lowest with the true vacuum.The critical point corresponds to the location at which the lines from the true and false vacua meet, resulting in a non-analyticity of the minimum of the free energy as a function of t.Yang-Mills gauge theory on a 4 × 20 3 lattice, computed in the LLR method with energy intervals of size δ E /6 Ṽ = 0.000374281, in the energy range from E 0 /6 Ṽ = 0.439487341 to E 54 /6 Ṽ = 0.459698522.The free-energy is given by F = E − ts, where E = 6 Ṽ(1 − u p ) is the internal energy, s = ln ρ is the entropy and t = 1/a n is the temperature.Σ is a constant, computed as the temperature average of s − ln ρ 0 where ρ 0 originates in the density of states, Eq. ( 7).The inset shows a n and E n , for the corresponding points.

Conclusions
First-order phase transitions in the early universe can give rise to a stochastic background of gravitational waves that can be used to constrain new physics.For strongly coupled new physics models, the relevant observables can in principle be measured using lattice field theory.With the aim to obtain robust results, we have explored the possibility of applying the LLR method, whereby the density of states of the system is determined numerically in the microcanonical ensemble and then used to reconstruct observables performing one-dimensional numerical integrals.We have tested the method for a SU(3) gauge theory close to criticality, on a single lattice volume and at a relatively coarse lattice spacing.The method has been validated using conventional lattice calculations, which are still viable for the chosen lattice parameters.Our approach has also shown clear advantages in terms of the quantities it can give us access to, which include the free energy, a precisely determined probability distribution of the states at criticality and the dependency of the microcanonical temperature on the energy.Those quantities, which cannot be robustly determined in importance sampling approaches, show clear singularities that carry direct information about the phase transition.The material presented here is the first step of a programme that aims to use the LLR method to compute more accurately and efficiently thermodynamic quantities in non-Abelian gauge theories.
puting Wales project) and AccelerateAI A100 GPU system, and on the DiRAC Data Intensive service at Leicester.The Swansea SUNBIRD system and AccelerateAI are part funded by the European Regional Development Fund (ERDF) via Welsh Government.The DiRAC Data Intensive service at Leicester is operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk).The DiRAC Data Intensive service equipment at Leicester was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1.DiRAC is part of the National e-Infrastructure.Open Access Statement -For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

Figure 1 .
Figure 1.An sample of 20 RM trajectories for SU(3) pure gauge theory on a 4 × 20 3 lattice with periodic boundary conditions in all directions and an energy interval of size δ E /6 Ṽ = 0.000374281 centred at E n /6 Ṽ = 0.459324241.The inset shows the obtained a (m) n distribution at the truncation value m = 500.

,Figure 2 :
Figure 2: Thermodynamic observables measured with the LLR method (blue circles), compared to results from standard importance sampling (black triangles), for SU(3) Yang-Mills gauge theory on a 4 × 20 3 lattice.The LLR calculation uses energy intervals of size δ E /6 Ṽ = 0.000374281 in the energy range from E 0 /6 Ṽ = 0.439487341 to E 54 /6 Ṽ = 0.459698522.The blue curves are reconstructed observables from LLR method with a finer resolution in β, restricted to the region around the phase transition.Left panel: average plaquette ⟨u p ⟩ against the coupling β.Right panel: specific heat C V ≡ ⟨u 2 p ⟩ − ⟨u p ⟩ 2 against the coupling β.

,Figure 3 .
Figure 3. Results of the LLR analysis of SU(3) Yang-Mills gauge theory on a 4 × 20 3 lattice with energy intervals of size δ E /6 Ṽ = 0.000374281 in the energy range from E 0 /6 Ṽ = 0.439487341 to E 54 /6 Ṽ = 0.459698522.The top panel shows the values of the micro-canonical inverse temperatures a n against the centre of the energy interval E n , with a linear interpolation between the points.The bottom panel shows the reconstructed probability distribution P βc (E) of the energy E at the critical coupling β c .The horizontal dashed line shows the location of the critical coupling, and the vertical lines are the average plaquette values at which a n = β c , which correspond to the locations of the extrema of the probability distribution.

,Figure 4 .
Figure 4.The free-energy for SU(3) Yang-Mills gauge theory on a 4 × 20 3 lattice, computed in the LLR method with energy intervals of size δ E /6 Ṽ = 0.000374281, in the energy range from E 0 /6 Ṽ = 0.439487341 to E 54 /6 Ṽ = 0.459698522.The free-energy is given by F = E − ts, where E = 6 Ṽ(1 − u p ) is the internal energy, s = ln ρ is the entropy and t = 1/a n is the temperature.Σ is a constant, computed as the temperature average of s − ln ρ 0 where ρ 0 originates in the density of states, Eq. (7).The inset shows a n and E n , for the corresponding points.