Exceptional thermodynamics: The equation of state of G(2) gauge theory

We present a lattice study of the equation of state in Yang-Mills theory based on the exceptional G(2) gauge group. As is well-known, at zero temperature this theory shares many qualitative features with real-world QCD, including the absence of colored states in the spectrum and dynamical string breaking at large distances. In agreement with previous works, we show that at finite temperature this theory features a first-order deconfining phase transition, whose nature can be studied by a semi-classical computation. We also show that the equilibrium thermodynamic observables in the deconfined phase bear striking quantitative similarities with those found in SU(N) gauge theories: in particular, these quantities exhibit nearly perfect proportionality to the number of gluon degrees of freedom, and the trace anomaly reveals a characteristic quadratic dependence on the temperature, also observed in SU(N) Yang-Mills theories (both in four and in three spacetime dimensions). We compare our lattice data with analytical predictions from effective models, and discuss their implications for the deconfinement mechanism and high-temperature properties of strongly interacting, non-supersymmetric gauge theories. Our results give strong evidence for the conjecture that the thermal deconfining transition is governed by a universal mechanism, common to all simple gauge groups.


Introduction
Due to its highly non-linear, strongly coupled dynamics, analytical understanding of the strong nuclear interaction remains incomplete [1]. Essentially, the fact that the spectrum of physical states is determined by non-perturbative phenomena (confinement and chiral symmetry breaking) restricts the theoretical toolbox for first-principle investigation of QCD at low energies to numerical simulations on the lattice-while the applicability of weakcoupling expansions is limited to high-energy processes.
At present, one of the major research directions in the study of QCD (both theoretically and experimentally) concerns the behavior of the strong interaction under conditions of finite temperature and/or density. Asymptotic freedom of non-Abelian gauge theories suggests that, at sufficiently high temperatures, ordinary hadrons should turn into a qualitatively different state of matter, characterized by restoration of chiral symmetry and liberation of colored degrees of freedom, which interact with each other through a screened long-range force [2]: the quark-gluon plasma (QGP). After nearly twenty years of dedicated experimental searches through relativistic heavy-nuclei collisions, at the turn of the millennium the QGP was eventually discovered at the SPS [3] and RHIC [4][5][6][7] facilities.
These results indicate that the theoretical investigation of the QGP requires nonperturbative tools, such as computations based on the gauge/string correspondence (whose applications in QCD-like theories at finite temperature are reviewed in refs. [46,47]) or lattice simulations [48]. The lattice determination of the deconfinement crossover and chiral transition temperature, as well as of the QGP bulk thermodynamic properties (at vanishing quark chemical potential µ) is settled [49][50][51] and accurate results are being obtained also for various parameters describing fluctuations, for the QGP response to strong magnetic fields, et c. [52]. Due to the Euclidean nature of the lattice formulation, the investigation of phenomena involving Minkowski-time dynamics in the QGP is more challenging, but the past few years have nevertheless witnessed a lot of conceptual and algorithmic advances, both for transport properties [53] and for phenomena like the momentum broadening experienced by hard partons in the QGP [54][55][56][57][58][59][60][61][62][63]: it is not unrealistic to think that in the near future the results of these non-perturbative calculations could be fully integrated in model computations that provide a phenomenological description of experimentally observed quantities (for a very recent, state-of-the-art example, see ref. [64]).
Notwithstanding this significant progress towards more and more accurate numerical predictions, a full theoretical understanding of QCD dynamics at finite temperature is still missing. From a purely conceptual point of view, the problem of strong interactions in a thermal environment can be somewhat simplified, by looking at pure-glue non-Abelian gauge theories. This allows one to disentangle the dynamics related to chiral symmetry breaking from the problem of confinement and dynamical generation of a mass gap, retaining-at least at a qualitative or semi-quantitative level-most of the interesting features relevant for real-world QCD. As the system is heated up, these theories will interpolate between two distinct limits: one that can be modeled as a gas of massive, non-interacting hadrons (glueballs) at low temperature, and one that is described by a gas of free massless gluons at (infinitely) high temperature. These limits are separated by a finite-temperature region, in which deconfinement takes place.
In SU(N ) gauge theories, the phenomenon of deconfinement at finite temperature can be interpreted in terms of spontaneous breaking of the well-defined global Z N center symmetry [65] (see also ref. [66] for a very recent work on the subject), and is an actual phase transition: a second-order one for N = 2 colors, and a discontinuous one for all N ≥ 3 (see also refs. [67][68][69]). While this is consistent with the interpretation of confinement in non-supersymmetric gauge theories as a phenomenon due to condensation of center vortices [70,71] (see also ref. [72] for a discussion), it begs the question, what happens in a theory based on a non-Abelian gauge group with trivial center? In this respect, it is particularly interesting to consider the G 2 gauge theory: since this exceptional group is the smallest simply connected group with a trivial center, it is an ideal toy model to be studied on the lattice.
Note that, even though smaller continuous non-Abelian groups with a trivial center do exist, strictly speaking what actually counts is the fundamental group of the compact adjoint Lie group associated with the Lie algebra of the gauge group. For example, the SO(3) group has a trivial center, but (contrary to some inaccurate, if widespread, claims) this property, by itself, does not make the SO(3) lattice gauge theory a suitable model for studying confinement without a center, 1 nor one to be contrasted with the SU(2) gauge group (which has the same Lie algebra). Indeed, the first homotopy group of SO(3), which is the compact adjoint Lie group associated with the B 1 Lie algebra, is Z 2 , i.e. the same as the first homotopy group of the projective special unitary group of degree 2 (whose associated Lie algebra is A 1 = B 1 ). This leaves only G 2 , F 4 and E 8 as compact simply-connected Lie groups with a trivial center; of these, G 2 , with rank two and dimension 14, is the smallest and hence the most suitable for a lattice Monte Carlo study. In fact, numerical simulations of this Yang-Mills theory have already been going on for some years [74][75][76][77][78][79][80][81][82][83][84]. Besides numerical studies, these peculiar features of G 2 Yang-Mills theory have also triggered analytical interest [85][86][87][88][89][90][91][92][93][94][95].
Note that the question, whether G 2 Yang-Mills theory is a "confining" theory or not, depends on the definition of confinement. If one defines confinement as the absence of non-color-singlet states in the physical spectrum, then G 2 Yang-Mills theory is, indeed, a confining theory. On the other hand, if one defines confinement as the existence of an asymptotically linear potential between static color sources, then the infrared dynamics of G 2 Yang-Mills theory could rather be described as "screening". Indeed, previous lattice studies indicate that, at zero and low temperatures, the G 2 Yang-Mills theory has a confining phase, in which static color sources in the smallest fundamental irreducible representation 7 are confined by string-like objects, up to intermediate distances. At very large distances, however, the potential associated with a pair of fundamental sources gets screened. This is a straightforward consequence of representation theory (and, ultimately, of the lack of a non-trivial N -ality for this group): as eq. (A. 15) in the appendix A shows, the representation 7 appears in the decomposition of the product of three adjoint representations 14, thus a fundamental G 2 quark can be screened by three gluons.
One further reason of interest for a QCD-like lattice theory based on the G 2 group is that it is free from the so-called sign problem [96]: with dynamical fermion fields in the fundamental representation of the gauge group, the introduction of a finite quark chemical potential does not make the determinant of the Dirac matrix complex, thus the theory can be simulated at finite densities [97,98]. As compared to another well-known QCD-like theory which shares this property, namely two-color QCD [99][100][101][102][103], one advantage is that "baryons" in G 2 QCD are still fermionic states, like in the real world.
Finally, the possibility that gauge theories based on exceptional gauge groups may be relevant in walking technicolor scenarios for spontaneous electro-weak symmetry breaking was studied (via a perturbative analysis) in ref. [104].
In this work, we extend previous lattice studies of G 2 Yang-Mills theory at finite temperature [82][83][84] by computing the equation of state in the temperature range T 3T c , where T c denotes the critical deconfinement temperature. 2 After introducing some basic definitions and the setup of our simulations in section 2, we present our numerical results in section 3; then in section 4 we compare them with the predictions of some analytical calculations, pointing out qualitative and quantitative analogies with SU(N ) gauge theories. Finally, in section 5 we summarize our findings and list possible extensions of the present work. Some general properties of the G 2 group and of its algebra are reported in the appendix A.

Setup
Our non-perturbative computation of the equation of state in G 2 Yang-Mills theory is based on the standard Wilson regularization of the theory on a four-dimensional, Euclidean, hypercubic lattice Λ of spacing a [106]. Throughout this article, we denote the Euclidean time direction by the index 0 (or by a subscript t ), and the spatial directions by 1, 2 and 3 (or by a subscript s ). Periodic boundary conditions are imposed along the four directions. Using natural units c = = k B = 1, the physical temperature is given by the inverse of the length of the shortest side of the system (which we take to be in the direction 0), T = 1/(aN t ), while the other three sides of the hypertorus have equal lengths, denoted by L s = aN s . In order to avoid systematic uncertainties caused by finite-volume effects, we always take L s T 1; in practice, previous studies of SU(N ) Yang-Mills thermodynamics have shown that, at the temperatures of interest for this work, finite-volume effects are negligible for L s T 4 [107][108][109].
To extract vacuum expectation values at low temperature (in the confining phase), we carry out simulations on lattices of sizes L 4 s : for the parameters of our simulations, this choice corresponds to temperatures which are sufficiently "deep" in the confining phasemeaning temperatures, at which the values of the bulk thermodynamic quantities, that we are interested in, are well below the statistical precision of our data.
The partition function of the lattice system is defined by the multiple group integral where dU α (x) denotes the Haar measure for the generic U α (x) matrix, which represents the parallel transporter on the oriented bond from x to x + aα. The U α (x) matrices take values in the representation of the G 2 group in terms of real 7 × 7 matrices, 3 while is the gauge-invariant Wilson lattice action [106]. Here, g is the bare lattice coupling and denotes the plaquette stemming from site x and lying in the oriented (µ, ν) plane. In the following, we also introduce the Wilson action parameter β, which for this theory can be defined as β = 7/g 2 .
As usual, expectation values of gauge-invariant quantities O are then defined as and can be computed numerically, via Monte Carlo integration. To this purpose, we generated ensembles of matrix configurations using an algorithm that performs first a heathbath update, followed by five to ten overrelaxation steps, on an SU(3) subgroup of G 2 (in turn, both the heathbath and overrelaxation steps are based on three updates of SU(2) subgroups [111]). Finally, a G 2 transformation is applied, in order to ensure ergodicity. The parameterization of G 2 that we used is described in refs. [112,113]. After a thermalization transient, we generate the ensemble to be analyzed by discarding a certain number (which depends on the physical parameters of the simulation-in particular, on the proximity to the deconfinement temperature, which affects the autocorrelation time of the system) of intermediate configurations between those to be used for our analysis. Typically, the number of configurations for each combination of parameters (β, N s and N t ) is O(10 4 ). This leads to an ensemble of (approximately) statistically independent configurations, allowing us to bypass the problem of coping with difficult-to-quantify systematic uncertainties due to autocorrelations. Throughout this work, all statistical errorbars are computed using the gamma method; a comparison on a data subset shows that the jackknife procedure gives roughly equivalent results.
We computed expectation values of hypervolume-averaged, traced Wilson loops at zero temperature as well as of volume-averaged, traced Polyakov loops at finite temperature (where V t=0 denotes the spatial time-slice of Λ at t = 0) and of hypervolume-averaged plaquettes (both at zero and at finite temperature). From the expectation values of Wilson loops, which we computed with the multilevel algorithm [114], the heavy-quark potential V (r) can be extracted via In practice, since we are bound to use loops of finite sizes, in order to avoid possible contamination from an excited state, we perform both a single-(k 1 = 0) and a two-state (k 1 free) fit W (r, L) = e −LV (r) + k 1 e −LV 1 (r) , (2.9) extracting our results for V from fits in the L range where the results including or neglecting the second addend on the right-hand side of eq. (2.9) are consistent, within their uncertainties.
We set the scale by carrying out a two-parameter fit of V (r) to the Cornell form which is appropriate at the distances probed in this work (gluon screening sets in at much longer distances), and enables one to extract the values of the string tension in lattice units σa 2 , at each value of β. Note that the coefficient of the 1/r term in eq. (2.10) is uniquely fixed by the central charge of the underlying low-energy effective theory describing transverse fluctuations of the confining string in four spacetime dimensions [115], while we neglect possible higher-order (in 1/r) terms [116].
Alternative methods to set the scale are discussed in the recent refs. [117,118] and in the works mentioned therein.
The phase structure of the lattice theory is revealed by the expectation value of the plaquette-averaged over the lattice hypervolume and over the six independent (µ, ν) planes-at T = 0, that we denote as U p 0 : similarly to what happens in SU(N ) gauge theories, as β is increased from zero to large values, U p 0 interpolates between a strongcoupling regime, dominated by lattice discretization effects, and a weak-coupling regime, analytically connected to the continuum limit. The two regions are separated by a rapid crossover (or, possibly, a first-order transition) taking place at β c 9.45 [82,83]. This bulk transition is unphysical, and, in order to extract physical results for the thermodynamics of the theory, all of our simulations are performed in the region corresponding to "weak" lattice couplings, β > β c , which is connected to the regime of continuum physics.
In the region of weak couplings, the finite-temperature deconfinement transition is probed by studying the distribution of values and the Monte Carlo history of the bare Polyakov loop (after thermalization): in the confining phase, the distribution of P is peaked near zero, whereas a well-defined peak at a finite value of P develops, when the system is above the critical deconfinement temperature T c . In the vicinity of T c , both the distribution of P values ρ(P ) (with two local maxima, separated by a region in which ρ(P ) gets suppressed when the physical volume of the lattice is increased) and the typical Monte Carlo histories of P (featuring tunneling events, which become more and more infrequent when the lattice volume is increased, between the most typical values of P ) give strong indication that the deconfining transition is of first order, in agreement with previous studies [82,83].
The equilibrium thermodynamics quantities of interest in the present work are the pressure p, the energy density per unit volume , the trace of the energy-momentum tensor ∆ (which has the meaning of a trace anomaly, being related to the breaking of conformal invariance of the classical theory by quantum effects), and the entropy density per unit volume s: they can be obtained from the finite-temperature partition function Z via (2.14) Introducing also the free energy density the pressure can be readily computed using the p = −f identity, which holds in the thermodynamic limit. Using the standard "integral method" of ref. [119], p (or, more precisely, the difference between the pressure at finite and at zero temperature) can thus be obtained as where U p T denotes the plaquette expectation value at the temperature T , and β 0 corresponds to a point sufficiently deep in the confining phase, i.e. to a temperature at which the difference between p and its zero-temperature value is negligible. The integral in the rightmost term of eq. (2.16) is computed numerically, by carrying out simulations at a set of (finely spaced) β values within the desired integration range, and performing the numerical integration according to the trapezoid rule. Although more sophisticated methods, like those described in ref. [120, appendix A], could allow us to reduce the systematic uncertainties related to the numerical integration, it turns out that this would have hardly any impact on the total error budget of our results. The lattice determination of the trace of the energy-momentum tensor ∆ (in units of T 4 ) is even more straightforward, as but it requires an accurate determination of the relation between β and a. We carried out the latter non-perturbatively, using the σa 2 values extracted from the confining potential, as discussed above. 4 4 Note that eq. (2.17) reveals one technical challenge in this computation: on the one hand, as we mentioned above, the simulations have to be carried out at values of the coupling in the region analytically connected to the continuum limit, i.e. β > βc. In practice, this corresponds to relatively fine lattice spacings, or Nt 5. On the other hand, the N 4 t factor appearing on the right-hand side of eq. (2.17) shows that the physical signal is encoded in a difference of average plaquette values Up 0 and Up T (both of which remain finite), that becomes increasingly small when the continuum limit is approached. In practice, the computational costs due to this technical aspect severely restrict the range of Nt values which can be used, and, as a consequence, the lever arm to control the continuum limit.
Finally, and s can be readily computed as linear combinations of p and ∆, using eq. (2.13) and the thermodynamic identity It is worth noting that alternative methods to determine the equation of state have been recently proposed in ref. [121] and in refs. [122][123][124][125].

Numerical results
The first set of numerical results that we present in this section are those aimed to the non-perturbative determination of the scale, i.e. of the relation between the parameter β and the corresponding lattice spacing a. As discussed above, we derive this relation by extracting the string tension in lattice units from the area-law decay of large Wilson loops at zero temperature.
To achieve high precision, this numerical computation is carried out using the multilevel algorithm [114], which yields exponential enhancement of the signal-to-noise ratio for long loops. Fig. 1 shows results for the opposite of the logarithm of Wilson loops of different widths r/a (symbols of different colors) as a function of the loop length in lattice units L/a. The plot displays the results from our simulations at β = 10.4. The comparison of results obtained from a naïve, brute-force computation (empty circles) and with our implementation of the multilevel algorithm (filled squares) clearly shows that the latter are in complete agreement with the former for short loops, and that the multilevel algorithm outperforms the brute-force approach for large values of L, where the numerical values obtained with the latter are affected by dramatic loss of relative precision.
Our results for the string tension in lattice units, as a function of β, are reported in table 1, which also shows the sizes of the lattices on which the corresponding simulations were carried out (all these zero-temperature simulations were performed on hypercubic lattices of size L 4 s ), the range of r/a values used in the fit (including both extrema), as well as the results for the constant term in the Cornell potential in lattice units and the reduced χ 2 . Note that these are two-parameter fits to eq. (2.10), while, as mentioned in section 2, the 1/r term is fixed to be the Lüscher term [115]. In principle, for data at small r/a one could use an improved definition of the lattice distance [126,127], however this type of correction becomes rapidly negligible at large distances, and hence should not change significantly our estimates of σa 2 , which are dominated by infrared physics.
The values of σa 2 thus computed non-perturbatively can be interpolated by a fit to a suitable functional form, in order to get an expression for a as a function of β in the region of interest. In principle, this can be done in various ways (see, for example, refs. [127][128][129]), which, in particular, can include slightly different parametrizations of the discretization effects. One of the simplest possibilities is to fit our data for the logarithm of σa 2 to a polynomial of degree n par − 1 in (β − β 0 ), where β 0 is a value within the range of simulated data. Choosing β 0 = 10.2, a parabolic fit with three parameters, however, yields a large χ 2 red 6.4. Different choices of β 0 and/or of n par ≥ 3 give interpolating functions that are only marginally different (within the uncertainties of the fitted parameters) and do not r / a = 2, multilevel r / a = 3, multilevel r / a = 4, multilevel r / a = 5, multilevel r / a = 6, multilevel r / a = 7, multilevel r / a = 8, multilevel r / a = 2, naïve r / a = 3, naïve r / a = 4, naïve r / a = 5, naïve r / a = 6, naïve r / a = 7, naïve r / a = 8, naïve β = 10.4  bring χ 2 red down to values close to 1. For example, using a cubic, rather than quadratic, polynomial, the fitted curve changes slightly, and the χ 2 red (using all points) changes from 6.4 to 5.2. These unsatisfactory results hint at discretization effects-an indication confirmed by that fact that, excluding the data corresponding to the coarsest lattice spacing, at β = 9.6, from the fit, the χ 2 red goes down to 2.2-and call for a modelling of our lattice results via a functional form that could (at least partially) account for lattice cutoff systematics. Therefore, following ref. [128], we chose to interpolate the values for the string tension in lattice units computed non-perturbatively by a fit to with β 0 = 9.9. This yields c 1 = 0.139 (25), c 2 = 0.450 (99) and c 3 = −0.42 (10), with χ 2 red = 0.73. Among the systematic uncertainties affecting this scale setting (which we include in our final results) are, for example, those related to the possibility of adding a c 4 f 4 (β) term in the denominator, or modifying the functional form for f (for instance, multiplying it by a polynomial in β): while not necessarily better-motivated from a theoretical point of view, 5 these alternative parametrizations do not lead to significant changes in our final results, hence we keep eq. (3.1), with the parameters listed above, as our final determination of the scale. The results of this fit are shown in fig. 2, where, in addition to our results, we also show those obtained in ref. [81], which are essentially compatible with our interpolating curve (with its uncertainty).
From the results of our three-parameter fit to eq. (3.1), the lattice spacing a (and, hence, the physical scale of our simulations and the temperature) can be determined at any β value within the range of interpolation-provided one defines a physical value for the string tension σ: to make contact with real-world QCD, one can, for example, set σ = (440 MeV) 2 . In addition, the derivative of ln ( √ σa) with respect to β is also obtained as and the inverse of this quantity is later used in the computation of the trace of the energymomentum tensor ∆, according to eq. (2.17). Next, we proceed to simulations at finite temperature, which we carried out on lattices of sizes N 3 s × N t (in units of the lattice spacing), where the shortest size N t defines the tem- As we already pointed out, by virtue of thermal screening, an aspect ratio of the order of 4 (or larger) for the "temporal" cross-section of the system is known to provide a sufficient suppression of finite-volume effects in SU(N ) gauge theories at the temperatures under consideration [108], while sizable corrections to the thermodynamic quantities are expected to appear at much higher temperatures [107]. Some tests on lattices of different spatial volume confirm that this is the case for G 2 Yang-Mills theory, too, and did not give us any evidence of significant finite-volume corrections. The parameters of our simulations are summarized in table 2. In order to compute the pressure with respect to its value at a temperature close to zero, according to the method described in sect. 2, for each set of finite-temperature simulations we also carried out Monte Carlo simulations on lattices of sizes (aN s ) 4 , at the same values of the lattice spacing.
The first task consists in identifying, for each value of N t , the critical coupling corresponding to the transition from the confining to the deconfined phase: by varying β, the lattice spacing a can be tuned to 1/(N t T c ). As mentioned in sect. 2, the transition from one (blue histograms), at β = 9.765 (red histograms) and at β = 9.77 (green histograms), reveals the transition from a confining phase at low temperature, in which ρ(P ) has a peak close to zero, to a deconfining one at high temperature, in which ρ has a global maximum for a finite value of P .
phase to the other can be identified by monitoring how the distribution ρ(P ) of values of the spatially averaged, bare Polyakov loop P varies with β. The confining phase is characterized by a distribution peaked at zero, while in the deconfined phase ρ(P ) has a maximum at a finite value of P , and the transition (or crossover) region can be identified as the one in which ρ(P ) takes a double-peak structure, with approximately equal maxima. An example of such behavior is shown in fig. 3, where we plotted the distribution of values of Polyakov loops from lattices with N t = 6 and N s = 32, at three different β values, namely 9.76, 9.765 and 9.77 (corresponding to three different values of the lattice spacing, and, hence, of the temperature). More precisely, in any finite-volume lattice this identifies a pseudo-critical coupling: as usual, the existence of a phase transition is only possible for an infinite number of degrees of freedom, namely in the thermodynamic, infinite-volume, limit. Thus, the actual critical point corresponding to the thermodynamic phase transition is obtained by extrapolation of the pseudo-critical couplings to the infinite-volume limit.
A more accurate way to determine the location and nature of the transition is based on the study of the Binder cumulant [130], defined as which is especially useful in computationally demanding problems (see ref. [131, appendix] for an example). For a second-order phase transition, the values of B interpolate between two different limits at "small" and "large" β (i.e. at low and at high temperature, respectively). As the system volume is increased, B(β) tends to become a function with a sharper and sharper increase in the region corresponding to the critical β. The critical coupling in the thermodynamic limit can thus be estimated from the crossing of these curves. On the other hand, for a first-order phase transition B(β) develops a deep minimum near the transition point [132] (see also ref. [133] for a discussion).
We studied B for different values of the simulations parameters: one example is shown in figure 4, which refers to simulations on lattices with N t = 8 sites in the Euclidean time direction, at different values of β and for different spatial volumes. Note that the values of B vary rapidly within a narrow β-interval, allowing one to get a rather precise indication of the (pseudo-)critical point.
For our present purposes, however, the main qualitative features of the thermal deconfinement transition in G 2 Yang-Mills theory are already revealed by how the Monte Carlo history of the spatially averaged Polyakov loop and the ρ(P ) distribution vary with β and the lattice volume. For β values equal to (or larger than) the critical one, the former exhibits tunneling events between different vacua, which become increasingly rare when the lattice volume grows. This suggest that the passage from the confining to the deconfined phase is a transition of first order. Accordingly, the peaks in the ρ distribution at criticality tend to become separated by an interval of P values, whose probability density is exponentially suppressed when the lattice volume increases.
As an example, in fig. 5 we show the distribution of P values obtained from simulations at the critical point for fixed N t = 8 and fixed lattice spacing, for different spatial volumes (12a) 3 and (16a) 3 .
Our observation of a first-order deconfining transition confirms the results of earlier lattice studies of this model [82,83].
Having set the scale and determined the critical coupling for different values of N t , we proceed to the computation of equilibrium thermodynamic quantities at different lattice spacings, and to the discussion of their extrapolation to the continuum limit. As pointed out in sect. 2, the static observables of interest in this work are related to each other by elementary thermodynamic identities. Since our numerical determination of the equation of state is based on the integral method introduced in ref. [119], the quantity which is computed most directly is the trace of the energy-momentum tensor ∆: as shown by eq. (2.17), it is just given by the difference between the expectation values of the plaquette at zero and at finite temperature, up to a β-dependent factor. The results from our simulations (at different values of N t ) for the dimensionless ∆/T 4 ratio are shown in fig. 6, as a function of the temperature (in units of the critical temperature). Note that the data from simulation ensembles at different N t are close to each other, indicating that discretization effects are under good control.
To get results in the continuum limit, we first interpolated our results from each N t data set with splines, and then tried to carry out an extrapolation to the a → 0 limit, by fitting the values interpolated with the splines (at a sufficiently large number of temperatures) as a function of 1/N 2 t , including a costant and a linear term. The first of these two steps is done as follows: we split the data sets in n int intervals defined by n int −1 internal knots, and computed n int +3 (or n int +2) basis splines (B-splines): they define a function basis, such that every spline can be written as a linear combination of those. The systematics uncertainties involved in the procedure are related to the choice of the number of knots, and to the spline degree (quadratic or cubic); these uncertainties can be estimated by comparing the χ 2 values obtained for different choices, and turn out not to be large.
As for the second step (the pointwise extrapolation to the continuum limit), however, we observed that it leads to results which are mostly compatible with the curve obtained from the interpolation of the N t = 6 data set, except for a few (limited) regions, in which the extrapolation is affected by somewhat larger errorbars-an effect likely due to statistical fluctuations in the ensemble obtained from the finest lattice, which tend to drive the continuum extrapolation. Since the latter effects are obviously unphysical, for the sake of clarity of presentation we decided to consider the curve obtained from interpolation of our N t = 6 data (with the associated uncertainties) as an estimate of the continuum limit. This curve corresponds to the brown band plotted in fig. 6.
Our results for the trace anomaly reveal two very interesting features: 1. When expressed per gluon degree of freedom, i.e. dividing by 2 × d a (where 2 is the number of transverse polarizations for a massless spin-1 particle in 3 + 1 spacetime dimensions, and d a is dimension of the gluon representation, i.e. d a = 14 for G 2 , and d a = N 2 − 1 for SU(N ) gauge group), the results for ∆/T 4 agree with those obtained in SU(N ) Yang-Mills theories.
2. In the deconfined phase (at temperatures up to a few times the deconfinement temperature), ∆ is nearly perfectly proportional to T 2 .
This is clearly exhibited in fig. 7, where our lattice results for ∆/(2d a T 4 ) at N t = 6 are plotted against (T c /T ) 2 , together with analogous results for the SU(3) and SU(4) theories (at N t = 5) from ref. [134]. 6 The collapse of data obtained in theories with different gauge groups is manifest, as is the linear dependence on 1/T 2 in the temperature range shown (implying that ∆ is approximately proportional to T 2 ).
These features were already observed in SU(N ) gauge theories, both in four [134][135][136][137][138] and in three [139,140] spacetime dimensions (the latter provide an interesting theoretical laboratory: see, e.g., ref. [141] and references therein).  The errorbars shown in the plot account for the statistical uncertainties in the computation of the average plaquettes in eq. (2.17), as well as for statistical and systematic uncertainties related to the non-perturbative determination of the scale and of the critical coupling for each β. The brown band denotes the interpolation of our results from the ensembles corresponding to N t = 6. As discussed in the text, such curve turns out to be essentially compatible with the results obtained from an attempt to carry out the continuum extrapolation (up to small deviations in the latter, which are likely due to statistical effects). Thus, we present the brown curve as an approximate estimate of the continuum limit.
Integrating the plaquette differences used to evaluate ∆/T 4 , the pressure (in units of T 4 ) is then computed according to eq. (2.16) for each N t . In principle, one could then extrapolate the corresponding results to the continuum limit; however, like for the trace anomaly, it turns out that, at the level of precision of our lattice data, this leads to results which are essentially compatible with those from our N t = 6 ensemble (within uncertainties, including those related to the extrapolation systematics). Therefore, in fig. 8 we show the results for p/T 4 obtained by numerical integration of the curve interpolating the N t = 6 data (solid red curve): this curve can be taken as an approximate estimate of the continuum limit (up to an uncertainty defined by the band within the dashed red curves). As one can see, at the highest temperatures probed in this work the pressure is growing very slowly (due to the logarithmic running of the coupling with the typical energy scale of the thermal ensemble, which is of the order of T ) and tending towards its value in the Stefan-Boltzmann   [134] at N t = 5). All G 2 results displayed here were obtained from simulations at N t = 6. Note that in this figure the data are plotted against (T c /T ) 2 , in order to reveal the approximately perfect proportionality between ∆ and T 2 in the deconfined phase (at least in the temperature range investigated in this work, up to a few times T c ): with this choice of axes, this feature manifests itself in the linear behavior observed in the plot. 4) so that at temperatures T 3T c the system is still relatively far from a gas of free gluons. Fig. 8 also shows our results for the energy density in units of the fourth power of the temperature ( /T 4 , solid blue curve) and for the entropy density in units of the cube of the temperature (s/T 3 , solid green curve), respectively determined according to eq. (2.13) and to eq. (2.18). Like for the pressure, the uncertainties affecting these two quantities are denoted by the bands enclosed by the dashed curves.

Discussion
The features of this exceptional group (in particular: the fact that it has a trivial center) make the G 2 Yang-Mills model very interesting for a comparison with gauge theories based on classical simple Lie groups. As we discussed, previous works already showed that at zero temperature this model bears several qualitative similarities with QCD: the physical spectrum does not contain colored states, and the potential associated with a pair of static color sources is linearly rising at intermediate distances-before flattening out at very large distances, due to dynamical string-breaking. However, a difference with respect to QCD (in which the color charge is screened by creation of dynamical quark-antiquark pairs, which are absent in pure Yang-Mills theory) is that in G 2 Yang-Mills theory screening is due to gluons. At finite temperature, there is numerical evidence that this theory has a first-order deconfinement phase transition (at which the average Polyakov loop modulus jumps from small to finite values), even though this transition is not associated with the breaking of center symmetry [82].
Our lattice results confirm the first-order nature of the deconfinement transition in G 2 Yang-Mills theory at finite temperature. This is in agreement with analytical studies available in the literature. In particular, a semiclassical study of the confinement/deconfinement mechanism in different Yang-Mills theories was presented in ref. [85] (a related study, discussing the inclusion of fundamental fermionic matter, is reported in ref. [91], while a generalization to all simple Lie groups has been recently presented in ref. [94]). Generalizing a previous study for the SU(2) case [144], the authors of ref. [85] showed how the properties of the high-temperature phase of a generic Yang-Mills theory can be understood, by studying its N = 1 supersymmetric counterpart on R 3 × S 1 , and by continuously connecting the supersymmetric model to the pure Yang-Mills theory, via soft supersymmetry-breaking induced by a finite gluino mass. This analytical study is possible, by virtue of the fact that, when the S 1 compactification length L becomes small, the theory can be reliably investigated by means of semi-classical methods. 8 In particular, analyzing the effective potential describing the eigenvalues of the Polyakov line, it turns out that: • The phase transition is driven by the competition between terms of perturbative origin [145][146][147][148], Bogomol'nyi-Prasad-Sommerfield monopole-instantons, and Kaluza-Klein monopole-instantons (which tend to make the Polyakov line eigenvalues collapse, namely to break center symmetry) and neutral bions (that stabilize the center). 9 • This confining/deconfining mechanism is common to all non-Abelian theories, irrespective of the underlying gauge group. The order of the deconfining transition, however, does depend on the gauge group: for the SU(2) case, the mechanism predicts the existence of a second-order transition, whereas for SU(N ≥ 3) and for G 2 the transition is a discontinuous one.
In addition, the results of our lattice simulations also show that the equilibrium thermodynamic observables in this theory are quantitatively very similar to those determined in previous studies of the SU(N ) equation of state [134,136,137]. In fact, rescaling the various thermodynamic quantities by the number of gluon degrees of freedom, we found that the observables per gluon component in the deconfined phase of G 2 Yang-Mills theory are essentially the same as in SU(N ) theories. This is consistent with the observation (based on an analysis of the gluon propagator in Landau gauge) that confinement and deconfinement should not be qualitatively dependent on the gauge group [87]. The same independence from the gauge group has also been observed in the numerical study of Polyakov loops in different representations in SU(N ) gauge theories [109]: the quantitative similarities with results in the SU(3) theory [152][153][154] are suggestive of common dynamical features.
In particular, our data show that, in the deconfined phase of G 2 Yang-Mills theory, the trace of the stress-energy tensor ∆ is nearly exactly proportional to T 2 for temperatures up to a few times the critical deconfinement temperature T c . This peculiar behavior was first observed in SU(3) Yang-Mills theory [135] (see also ref. [138] for a more recent, highprecision study) and discussed in refs. [155][156][157][158][159]. Later, it was also observed in numerical simulations of gauge theories with SU(N > 3) [134,136,137]. A dependence on the square of the temperature is hard to explain in perturbative terms (unless it accidentally results from cancellations between terms related to different powers of the coupling). Actually, at those, relatively low, temperatures, most likely the gluon plasma is not weakly coupled and non-perturbative effects probably play a non-negligible rôle [160,161]. While one could 8 Note that, since periodic boundary conditions are imposed along the compactified direction for all fields (including fermionic ones), the transition in the supersymmetric theory is a quantum-rather than a thermal -one. 9 For further details about these topological objects, see also refs. [149][150][151] and references therein.
argue that this numerical evidence in a relatively limited temperature range may not be too compelling, it is interesting to note that lattice results reveal the same T 2 -dependence also in 2 + 1 spacetime dimensions [139,140] (for a discussion, see also ref. [162] and references therein). Note that, in the latter case, due to the dimensionful nature of the gauge coupling g, the relation between g 2 and the temperature is linear, rather than logarithmic. Another interesting recent analytical work addressing the thermal properties of G 2 Yang-Mills theory was reported in ref. [92] (see also ref. [95]): following an idea discussed in refs. [90,155,163,164], in this article the thermal behavior of the theory near T c is assumed to depend on a condensate for the Polyakov line eigenvalues, and the effective action due to quantum fluctuations in the presence of this condensate is computed at two loops. The somewhat surprising result is that the two-loop contribution to the effective action is proportional to the one at one loop: this holds both for SU(N ) and for G 2 gauge groups. In order to derive quantitative predictions for the pressure and for the Polyakov loop as a function of the temperature, however, non-perturbative contributions should be included, as discussed in ref. [90]. More precisely, the effective description of the deconfined phase of Yang-Mills theories presented in ref. [90] is based on the assumption that, at a given temperature, the system can be modelled by configurations characterized by a constant (i.e. uniform in space) Polyakov line, and that the partition function can be written in terms of an effective potential for the Polyakov line eigenvalues. By gauge invariance, the Polyakov line can be taken to be a diagonal matrix without loss of generality. For SU(N ) gauge groups, the N eigenvalues of this matrix are complex numbers of modulus 1. Since the determinant of the matrix equals 1, the eigenvalues' phases are constrained to sum up to 0 mod 2π. It is convenient to define rescaled phases (to be denoted as q), that take values in the real interval between −1 and 1; one can then write an effective potential V(q, T ), which includes different types of contributions (of perturbative and non-perturbative nature). For the G 2 gauge group, the construction exploits the fact that G 2 is a subgroup of SO (7), which, in turn, is a subgroup of SU (7). Starting from SU (7), these conditions reduce the number of independent components of q down to 2, the rank of G 2 . Thus, the effective potential can take the form where the first term on the right-hand side simply gives the free energy density for the free gluon gas (according to eq. (3.4) and to the p = −f identity), the next term is the leading perturbative contribution, which can be expressed in terms of Bernoulli polynomials, while the terms within the square brackets are expected to mimic effects relevant close to the deconfinement transition (note the presence of the coefficient proportional to T 2 c , related to non-perturbative physics, in front of the square brackets): see ref. [90] for the precise definitions and for a thorough discussion.
The effective potential in eq. (4.1) depends on the unknown coefficients c 1 , c 2 , c 3 and d 2 , which, in principle, could be fixed using our data. We carried out a preliminary study in this direction, finding that (with certain, mild assumptions) it is indeed possible to fix values for c 1 , c 2 , c 3 and d 2 yielding ∆/T 4 values compatible with our lattice data. However, the quality of such determination is not very good, because the parameters appear to be cross-correlated and/or poorly constrained. Without imposing additional restrictions, an accurate determination of these parameters would require data of extremely high precision (and a genuine, very accurate continuum extrapolation).
After fixing the parameters of the effective potential defined in eq. (4.1), it would be interesting to test, whether the model correctly predicts the behavior of other observables computed on the lattice near the deconfining transition: this is a task that we leave for the future.

Conclusions
The present lattice study of the G 2 Yang-Mills model at finite temperature confirms that this gauge theory has a finite-temperature deconfining phase transition. In agreement with earlier lattice studies [82,83], we found that the latter is of first order, as predicted using semi-classical methods applicable to all simple gauge groups [85]. In particular, the nature of the deconfinement transition, determined by the form of the effective potential experienced by the Polyakov loop eigenvalues, results from the competition of different topological objects (and perturbative effects [145]): neutral bions (which generate repulsion among the eigenvalues) and magnetic bions, as well as monopole-instantons (which generate attraction among eigenvalues, like the perturbative terms).
The study of the equation of state that we carried out also shows that the equilibrium thermal properties of G 2 gauge theory are qualitatively and quantitatively very similar to those of all SU(N ) theories (up to a trivial proportionality to the number of gluon degrees of freedom), and are compatible with the predictions from recent analytical studies, like the one reported in refs. [90,92]. Recently, analogous conclusions have also been reached for supersymmetric theories [165,166], using an approach inspired by ref. [167]. These results corroborate the idea of universality in the thermal behavior of confining gauge theories. To summarize with a pun, one could say that the exceptional thermodynamics in the title of the present paper "is not so exceptional, after all".
Our findings are also interesting to understand the rôle that different topological excitations play in confinement, and give indications about the non-perturbative dynamics relevant at temperatures close to deconfinement, where truncated weak-coupling expansions are no longer quantitatively accurate.
This work could be generalized along various directions. The temperature range that we investigated could be extended, possibly in combination with the technical refinement of using an improved version of the gauge action, as was done for SU(3) Yang-Mills theory in ref. [138]. It is worth remarking that the multilevel algorithm used to set the scale in the present work has been recently generalized to improved actions [168]. With sufficient computational power, it would be interesting to compare the behavior of the thermodynamic quantities in the confined phase with a gas of free glueballs, possibly modeling the spectrum of excited states in terms of a vibrating bosonic string. This type of comparison was successfully carried out in ref. [169] for SU(3) Yang-Mills theory in 3 + 1 dimensions and in ref. [170] for SU(N ) theories in 2 + 1 dimensions.
One could also extend the investigation of the theory, by looking at other observables in the deconfined phase, in particular going beyond those characterizing the equilibrium properties of the QCD plasma. While the present study addresses a model which is interesting as a theoretical dōjō, but which is not physically realized in nature, ultimately our aim is to achieve a deeper understanding of phenomenologically relevant aspects of strong interactions at finite temperature. In particular, transport properties describing the real-time evolution of conserved charge densities in the QGP are of the utmost relevance for experimentalists and theorists alike. The lattice investigation of these quantities, however, is particularly challenging (see ref. [53] for a detailed discussion), and until recently has been mostly limited to the SU(3) theory. Given the aspects of universality that seem to emerge from the present study and from previous works, it would be interesting to investigate, whether also the spectral functions related to different transport coefficients in G 2 Yang-Mills are qualitatively and quantitatively similar to those extracted in the SU(3) theory-albeit this may prove computationally very challenging.
Another possible generalization would be to investigate the equation of state in a Yang-Mills theory based on another exceptional gauge group. The "most natural" candidate would be the one based on F 4 : this group has rank 4 and dimension 52 and, like G 2 , its center is trivial. The smallest non-trivial irreducible representation of this group, however, is 26dimensional, making lattice simulations of this Yang-Mills theory much more demanding from a computational point of view. We are not aware of any previous lattice studies of F 4 gauge theory.

Acknowledgments
This work is supported by the Spanish MINECO (grant FPA2012-31686 and "Centro de Excelencia Severo Ochoa" programme grant SEV-2012-0249). We thank A. Dumitru, C. Korthals Altes and R. Pisarski for helpful discussions. The simulations were performed on the PAX cluster at DESY, Zeuthen.
A General properties of the G 2 group and of its algebra In this appendix we summarize some basic facts about the G 2 group and its algebra. Our discussion mostly follows ref. [74], although (where appropriate) we also provide some additional technical details-in particular as it concerns the representation theory. G 2 is the smallest of the five exceptional simple Lie groups, with dim G 2 = 14. It is a subgroup of SO(7) and coincides with the automorphism group of the octonions. Equivalently, it can be defined as the subgroup of GL(7) preserving the canonical differential 3-form (given by the canonical bilinear form taking the cross product of two vectors as its second argument). G 2 has an SU(3) subgroup, and G 2 /SU(3) is isomorphic to the six-dimensional sphere S 6 [110]. This allows one to decompose a generic G 2 element as the product of a matrix associated with an element of S 6 , times an SU(3) matrix: for an explicit realization, see ref. [82, appendix A]. Another subgroup of G 2 is SO(4) [171].
The Lie algebra of the G 2 generators has dimension 14 and rank 2: its Cartan matrix is so that the Π-system includes a long and a short root, at a relative angle 5π/6 (a unique property among all simple Lie algebras). There exist two fundamental representations, which are seven-and fourteen-dimensional, respectively. The weight diagram of the 7dimensional fundamental representation is given by the vertices of a regular hexagon, plus its center. The representation of dimension 14 is the adjoint representation: its weight diagram is given by the vertices of a hexagram, with the addition of two points at its center.
Tensor products of irreducible representations are not, in general, irreducible. However, they can be decomposed into sums of irreducible representations. For G 2 , the most straightforward algorithm to compute the decomposition of tensor products of irreducible representations is the one based on girdles (see ref. [173] and references therein), which can be briefly summarized as follows.
In the real vector space of dimension equal to the rank n of the group (R 2 in this case) with coordinates x 1 , . . . , x n , consider t distinct points P (i) = P The girdle ξ(λ 1 , . . . λ n ) of a representation of a group is then a particular set of points, with certain well-defined multiplicities: for G 2 , the girdles are irregular dodecagons, that are symmetric with respect to both the x 1 and x 2 axes, and whose vertices in the first (λ 1 , λ 2 ) Table 3. The smallest irreducible representations of the G 2 group, sorted by increasing dimension d, up to 10 5 . 1 denotes the trivial representation, while 7 and 14 denote the two fundamental representations (14 being the adjoint representation). Non-equivalent irreducible representations of the same dimension are distinguished by a prime sign (we conventionally choose to use the prime sign for the representation with the largest value of λ 1 ). Note that, in contrast to the claim of ref.
The G 2 Casimir operators are discussed in ref. [174]; in particular, the functionally independent ones are those of degree 2 and 6. Their eigenvalues can be found in ref. [175, section 5].
The non-perturbative aspects of a non-Abelian gauge theory are related to the topological objects that gauge field configurations can sustain. For the G 2 group, the homotopy groups are listed in table 5-see ref.
[176]-, where Z 1 denotes the trivial group. G 2 is connected, with a trivial fundamental group; its second homotopy group is trivial, too, while the third is the group of integers, hence G 2 gauge theory admits "instantons". n π n (G 2 ) n π n (G 2 ) n π n (G 2 ) 0 Table 5. Lowest homotopy groups of the G 2 group, from ref. [176].
Finally, using the properties of exact sequences (and the additivity of homotopy groups of direct products of groups), it is also trivial to show that π 2 (G 2 /[U(1) × U(1)]) is Z ⊕ Z: this implies that, like for SU(3) gauge theory, when the global G 2 group gets broken to its Cartan subgroup, two types of 't Hooft-Polyakov monopoles appear.