Thermodynamics of SU(N) Yang-Mills theories in 2+1 dimensions II - The deconfined phase

We present a non-perturbative study of the equation of state in the deconfined phase of Yang-Mills theories in D=2+1 dimensions. We introduce a holographic model, based on the improved holographic QCD model, from which we derive a non-trivial relation between the order of the deconfinement phase transition and the behavior of the trace of the energy-momentum tensor as a function of the temperature T. We compare the theoretical predictions of this holographic model with a new set of high-precision numerical results from lattice simulations of SU(N) theories with N=2, 3, 4, 5 and 6 colors. The latter reveal that, similarly to the D=3+1 case, the bulk equilibrium thermodynamic quantities (pressure, trace of the energy-momentum tensor, energy density and entropy density) exhibit nearly perfect proportionality to the number of gluons, and can be successfully compared with the holographic predictions in a broad range of temperatures. Finally, we also show that, again similarly to the D=3+1 case, the trace of the energy-momentum tensor appears to be proportional to T^2 in a wide temperature range, starting from approximately 1.2 T_c, where T_c denotes the critical deconfinement temperature.


Introduction and motivation
The phase diagram of strongly interacting matter has been studied in an extensive experimental program since the 1980's. During the last decade, the main heavy ion collision facilities have provided convincing evidence for the existence of a new state of matter, which is qualitatively different from usual hadronic matter, and appears to behave as a nearly ideal fluid [1]. While these results confirm the intuitive theoretical expectation that, at high temperatures or densities, asymptotic freedom leads to deconfinement, i.e. to the liberation of colored particles from hadrons [2], they also reveal that the deconfined plasma is quite different from a gas of nearly free quarks and gluons, and should be rather described as a strongly coupled fluid. This makes the theoretical description of the equation of state of this system at temperatures close to deconfinement particularly challenging: weak-coupling computations [3] have to be pushed to high orders, but the convergence of perturbative expansions in thermal gauge theories is generally poor, and the evaluation of high-order terms is complicated by the appearence of severe infrared (IR) divergences [4]. The latter reveal the mathematically non-trivial structure of perturbative expansions in finite-temperature QCD (with terms which are non-analytical in α s ), and are related to the existence of an ultra-soft, chromomagnetic energy scale, which retains an intrinsically non-perturbative nature, and to long-wavelength modes that are strongly coupled at all temperatures.
As a consequence, numerical computations on the lattice are the main tool to derive the predictions of QCD at the temperatures probed in experiments, and in the last few years, various collaborations have presented results for the QCD equation of state at vanishing chemical potential, obtained from simulations including dynamical quarks at or close to the physical point [5]. 1 At the same time, high-precision lattice results have also been obtained for various equilibrium thermodynamic properties in SU(N ) Yang-Mills theories with a large number of colors N , which is a particularly interesting limit, for several reasons. 2 First of all, the large-N limit at fixed 't Hooft coupling λ = g 2 N and fixed number of flavors N f [9] provides a natural interpretation for some non-trivial features of QCD (such as, for instance, the OZI rule [10]), and leads to a topological classification of Feynman diagrams, in which the dominant contributions come from planar graphs. This is suggestive of an analogy with similar expansions in closed string theory [11]. Moreover, for N → ∞ one expects that all correlation functions of gauge-invariant operators factorize, and that the functional integral describing a large-N gauge theory should be dominated by a single "master" gauge field [12]; the translational invariance properties of the latter are related to the ideas of large-N volume independence [13].
As it concerns the phase diagram of QCD-like theories, it is interesting to note that the 1 By contrast, the progress in lattice simulations of the QCD equation of state at finite net baryon density has been slower, due to the existence of a severe sign problem [6]. As a consequence, in this case one often obtains useful insight from the numerical study of appropriate effective models, see, e.g., ref. [7].
2 Some of the surprising mathematical properties arising in the large-N limit of a generic quantum theory are related to the fact that, generally, the large-N limit can be interpreted as a sort of "classical limit". The meaning of this statement is made precise in ref. [8], with the definition of an appropriate basis of coherent states and a classical Hamiltonian. The construction, however, can be carried out explicitly, and leads to an exact solution, only for certain particularly simple models.
large-N limit leads to interesting implications at finite density, including, in particular, a possible "quarkyonic phase" [14]. Moreover, this limit is also important for applications of the conjectured correspondence between gauge and string theories [15] to study the strongly interacting plasma [16], since these computations are done in the infinite-N limit.
For these reasons, it is interesting to understand how much the thermal properties of strongly interacting gauge theories depend on the number of colors. Recent lattice simulations of large-N Yang-Mills theories at finite temperature [17,18] have revealed that, in the deconfined phase, the bulk thermodynamic observables are essentially independent of N (except for a trivial proportionality to the number of gluons), and that SU(3) [19] is close to the large-N limit. 3 This result is particularly interesting from a theoretical point of view, since it provides support to analytical studies of the QCD plasma based on the approximation of an infinite number of colors, and, furthermore, it can shed light onto the effective degrees of freedom relevant for the plasma near deconfinement, and/or rule out possible effective models.
A different perspective on the hot QCD plasma is based on the study of non-Abelian gauge theories in a lower-dimensional spacetime. While the case of D = 1 + 1 dimensions is essentially trivial [21], in D = 2 + 1 these theories are characterized by rich dynamics [22], and share many qualitative features with their D = 3 + 1 analogues. In particular, at low energies they are characterized by a spectrum of color-singlet states with a finite mass-gap and by linear confinement, and they exhibit a deconfining transition at a finite temperature T c . Moreover, these theories can also be studied using techniques inspired by the AdS/CFT correspondence: comparing the results obtained from first-principle lattice computations with those derived from the gauge/string duality can provide a useful test-bed for the application of holographic methods to study strongly coupled systems in D = 2 + 1 dimensions, such as those relevant for condensed matter systems at criticality [23]. Finally, one further motivation to look at the Yang-Mills equation of state in 2 + 1 dimensions stems from the observation that, in D = 3 + 1, the trace of the energy-momentum tensor in the deconfined phase appears to be proportional to T 2 over a broad range of temperatures [18,19,24]. This behavior seems to be at odds with the expectation from perturbative computations, which would rather predict a logarithmic dependence on the temperature. In order to understand the physical origin of this characteristic behavior in the physical case of D = 3 + 1, it is instructive to investigate whether the same phenomenon also occurs for a generic SU(N ) Yang-Mills theory in the lower-dimensional setup, given that the theories in 3 + 1 and in 2 + 1 dimensions have both some similarities and some obvious qualitative differences.
For these reasons, in this work we present a systematic study of finite-temperature SU(N ) Yang-Mills theories in D = 2+1 dimensions. Having already discussed the confining phase of these theories in a previous work [25], in the present article we focus on the equilibrium thermodynamic properties in the deconfined phase. In particular, we investigate the strongly coupled regime (close to the deconfinement temperature T c ), where physical quantities cannot be reliably caculated via weak-coupling expansions, and compare the theoretical results obtained from two different non-perturbative approaches: holographic computations based on the gauge/string correspondence, and numerical simulations in the lattice regularization.
The structure of this paper is as follows. First, in section 2 we construct a holographic model, which is expected to describe the deconfined finite-temperature phase of strongly coupled SU(N ) gauge theories in D = 2 + 1. Then, in section 3, we review the basic properties of SU(N ) Yang-Mills theories in 2 + 1 spacetime dimensions, and introduce their regularization on a Euclidean lattice. Next, in section 4 we present a set of high-precision results from our numerical lattice simulations of these theories, for different values of N ranging from 2 to 6. After discussing the extrapolation to the thermodynamic and continuum limits, we investigate the dependence of the pressure (p), of the trace of the energy-momentum tensor (or interaction measure, denoted by ∆), and of the energy ( ) and entropy (s) densities on the number of colors, and compare the prediction for ∆ to the holographic model. Section 5 includes a discussion of our findings and their implications. The appendix A reports the details of the computation of the lattice Stefan-Boltzmann limit in D = 2 + 1 and D = 3 + 1 dimensions.

A holographic model 2.1 Generalities
The gauge-gravity correspondence ("holography" for short) [15] successfully reproduces most of the salient features of large-N gauge theories. Although the first examples of holography involved supersymmetric and conformal quantum field theories, the correspondence was soon generalized to more realistic examples, including theories with linear confinement and no supersymmetry [26], hence in the same class as QCD. However, these models (sometimes referred to as models built in a "top-down approach"), that stem from D-brane constructions in type IIA or IIB string theory, generally have an infinite number of undesired scalar operators in their spectrum, arising from the Kaluza-Klein modes on the internal extra-dimensions of the ten-dimensional parent theory.
In the meantime, an alternative "bottom-up" holographic approach has been developed [27]: it uses minimal ingredients to model the desired features of confining gauge theories on the gravity side, in a more direct and "economic" fashion. Our approach here consists of an advanced version of the bottom-up construction, that is known as "improved holographic QCD" (IHQCD) [28,29]see ref. [30] for a review, and ref. [31] for similar constructions in the literature. In what follows, we first explain and review the setup of IHQCD.
Holography generally associates the energy dependence of the field theory with a radial direction r perpendicular to the D dimensions of the Minkowski spacetime in which the gauge theory is defined. Therefore the most economic "bottom-up" approach to a D-dimensional QFT involves a (D + 1)-dimensional gravitational background. In addition, the holographic correspondence associates a bulk field to each operator that is relevant or marginal in the IR. 4 For pure Yang-Mills theory, there are two such marginal operators: the energy-momentum tensor T µν and the gluon operator Tr F 2 /Λ 4−D where Λ is the dynamically generated energy scale of the theory (for D ≤ 4). The bulk fields that are dual to these operators are the metric g µν (r) and the dilaton field φ(r). The dependence of these fields on the radial coordinate r corresponds to the renormalization group scale dependence of the corresponding operators in the field theory. In particular, the profile of φ(r) encodes how the (dimensionless) coupling constant g 2 eff = g 2 /Λ 4−D runs with the energy. 5 On the other hand, in order to make φ(r) run with r in a dynamical gravitational setup, one needs to turn on a potential for it. Therefore, the minimal general relativity action of the IHQCD is: (where the ellipsis denotes some boundary counter-terms that should be introduced to render the variational problem on geometries with a boundary well-defined; we will not need the explicit form of these terms here). The coefficient ξ is an unspecified normalization constant 6 that will not play an important rôle in what follows. In particular, it can be absorbed into φ by a redefinition. Note that the action is proportional to N 2 and to a positive power of M P , which denotes a "reduced" Planck mass, 7 thus, in the large-N limit of the gauge theory, gravitational interactions are suppressed. For convenience, we keep the normalization of the scalar kinetic term unspecified, except that we assume ξ > 0. We also assume that the scalar potential has a single AdS minimum, that corresponds to the UV limit of the dual field theory: where is the AdS length scale.

Vacuum solution:
The solution to the action in eq. (1) that corresponds to the vacuum of the dual field theory is of the form: Here the scale factor of the metric b 0 (r) is of the AdS form b 0 (r) = /r only if the dilaton potential is constant: V = D(D − 1)/ 2 . More generally, when V is a non-trivial function of φ, the function b 0 attains the AdS form only in the UV, i.e. for r → 0, while it deviates from AdS in the IR limit, i.e. for r → ∞. Its profile can be determined by solving the Einstein's equations, given the potential V(φ). The IR part of the geometry characterizes the confinement properties in the vacuum of the dual field theory. In ref. [29], the various confining asymptotics were classified. In particular, we shall be interested in the confining geometries (with gapped and discrete spectrum) of the form: where the integration constant Λ corresponds to the dynamically generated energy scale, and the ellipsis denotes sub-leading terms that are typically logarithmic in r [29]. The parameter α is directly related to the large-φ asymptotics of the dilaton potential in eq. (1). Using Einstein's equations, it is straightforward to show that eq. (4) follows from a potential of the form: where ξ is the normalization factor appearing in front of the kinetic term in eq. (1). Therefore, α is a parameter of the action, rather than of the particular solution.
The glueball spectrum of the field theory can be obtained by solving the fluctuation equations of the dilaton and the metric, which can generally be expressed as a Schrödinger equation of the form −ψ (r) + V S (r)ψ(r) = m 2 ψ(r), where ψ denotes a generic glueball wave-function that corresponds to normalizable fluctuations of the bulk fields. 8 Here the Schrödinger potential V S has the following asymptotics: Therefore we observe that the Schrödinger potential is bounded both in the ultraviolet (UV) and in the IR limits, hence the glueball spectrum is gapped and discrete, if and only if α > 1 [29]. We will take α > 1 in the rest of our discussion. In particular, in the limit of large mass, one can use the WKB approximation to write down an approximate expression for the glueball spectra: where C is a constant that depends on the particular model and on the type of glueball.
Thermodynamics: The temperature is introduced by Wick-rotating the time direction t → iτ and compactifying the Euclidean time: τ ∼ τ + 1/T .
This geometry corresponds to the thermal ensemble in the confined phase, that we call the thermal gas (TG) solution. When evaluated on the TG solution, the action in eq. (1) yields the free energy of the dual theory. Neglecting the action counter-terms-which are denoted by the ellipsis in eq. (1)-, the result turns out to be divergent, due to the infinite volume of the asymptotic AdS space. Here we shall adopt a particular choice of renormalization, which corresponds to tuning the counter-terms, so that the on-shell thermal gas action vanishes. 9 The deconfined phase of the field theory is described by another solution with the same UV asymptotics, the black-hole (BH) background: The scale factor b(r) and the dilaton φ(r) are generally different from their counterparts in the thermal gas, eq. (8). The blackness function f (r) can be solved in terms of b(r) using Einstein's equations as: It is a monotonically decreasing function starting as f = 1 at r = 0 and vanishing at r = r h . The latter corresponds to the event horizon of the black-hole, where the 00-component of the BH metric in eq. (9) vanishes. The temperature of the deconfined gluonic ensemble is given by the Hawking temperature of the BH: The location of the horizon r h determines the temperature of the system. The entropy density of the ensemble is given by the Bekenstein-Hawking entropy: where we defined the entropy density dividing the total entropy by the total spatial volume and by the square of the number of colors. 10 In order to obtain s as a function of T , one has to invert the variables using eq. (11). The free energy of the deconfined phase is given by the value of the action in eq. (1), evaluated on the BH background. As we have fixed the counter-terms for the gravity action in eq. (1) by the requirement that they cancel the on-shell value of the TG action A(T G), the free energy is given by F = A(BH) − A(T G), hence the terms denoted by the ellipsis in eq. (1) can be ignored, because they cancel in the difference.
A practical way to determine the thermodynamics of the gluon plasma is as follows. The general relativistic system satisfies the first law of thermodynamics (see ref. [32] and references therein), therefore one can obtain the free energy density directly by integrating the entropy, Here the functions s(r h ) and T (r h ) are defined in eq. (12) and in eq. (11). The reason for the upper bound of integration is as follows. As the location of the horizon r h tends to infinity, the size of the BH becomes smaller and smaller and for r h → ∞ the BH and TG solutions coincide [32]. Therefore, the difference between the actions of the two solutions, hence the free energy F , should vanish. One can also notice this by comparing the metric functions b(r) and b 0 (r) of the BH and of the TG solutions in the limit r h → ∞, and observe that the difference b(r) − b 0 (r) vanishes exponentially in r h [32,33], in the entire region 0 ≤ r < r h . Once the free energy is determined as in eq. (13), one can change the variable to T instead of r h , using eq. (11), and obtain the free energy F (T ). Given F (T ), all other thermodynamic variables follow by standard thermodynamic identities. In particular, one can demonstrate [32] the existence of a Hawking-Page phase transition at a finite temperature T c , if the parameter α in eq. (4) satisfies the condition: α ≥ 1 [34]. The Hawking-Page transition corresponds to the confinement-deconfinement transition in the dual field theory.
We are particularly interested in the behavior of the interaction measure ∆ as a function of T . We define the normalized interaction measure by: where we used the standard thermodynamic relation: ∆ = (DF + ST )/V D−1 . The result can be immediately obtained, once the black-hole geometry is found, using the formulas given above. This is carried out explicitly in the next subsection.

Model construction
In this subsection, we shall discuss a general holographic construction that describes a normalized interaction measure∆ which decays with the temperature as 1/T in D = 2 + 1 dimensions. A more ambitious aim would be to construct a holographic model that fits all lattice data presented in section 4. This would certainly be possible by engineering the dilaton potential and fixing the parameters of the model, which are the parameters in the dilaton potential, the integration constants of the equations of motion for b(r) and φ(r), and M P . This aim was successfully achieved in the case of D = 3 + 1 [18,34,35]. Here, however, we restrict our investigation to the general behavior of the interaction measure and leave a more detailed holographic construction to future work. For this purpose, one can use the following simple, semi-analytic construction. Instead of starting from a given dilaton potential and obtaining the metric functions b(r) and f (r) by solving the Einstein's equations, one can simply make an Ansatz for the scale factor b(r), obtain f (r) from eq. (10), and derive the thermodynamic functions using the formulas presented in the previous subsection. Although this method is not exact but only approximate, 11 it turns out to provide a good approximation to thermodynamics in a large range of temperatures, as has been discussed in ref. [36]. The general conclusion, that we derive in the following, will be independent of the details of this approximation.
At this point, the question is: What Ansatz should one take for the scale factor b(r)? To answer this question, we first recall that b(r) tends to its TG analogue b 0 (r) both in the UV and in the IR limits, therefore eq. (4) leads to: where the ellipsis denotes subleading terms. A simple Ansatz that satisfies both limits is: 12 As discussed after eq. (4), the parameter α controls the properties of the theory in the IR. In particular, α should be larger than 1 for a confining theory. Interestingly, α also controls the nature of the confinement-deconfinement phase transition at T c : The transition tends to become a continuous one (in particular, a second-order one) as α → 1 [33,37]. Here we show that the parameter α also controls the decay of the interaction measure in the range of temperatures between T c and some intermediate value T i T c . The interaction measure is obtained from eq. (14) by numerical integration. 13 Fig. 1 shows the interaction measure ∆ (normalized dividing by T 3 N 2 ) as a function of the temperature T (in units of the deconfinement temperature T c ), for different values of α, in the case of D = 2 + 1 spacetime dimensions. We observe that the curve becomes steeper with increasing values of α; in section 4, we show a comparison of the curve corresponding to α = 3/2 to the numerical results obtained from lattice simulations, see fig. 3. One technicality in our present calculation is the value of T c . In general it is a number of the same order as Λ: with c 0 being some constant depending on the model. In order to determine c 0 (for a given Λ), we search for a value r c , such that for r h = r c the free energy difference vanishes: f (r c ) = 0. By definition, this point corresponds to the phase transition. The transition temperature T c is then obtained by evaluating eq. (11) at r h = r c .

General conclusions from holography
One can understand the dependence of the slope of the reduced interaction measure on the parameter α semi-analytically as follows. From eq. (14), one can write: Trace of the energy-momentum tensor from holography Figure 1: Comparison of the interaction measure ∆ (normalized by T 3 and by N 2 ) that follows from the holographic construction, for different choices of the parameter α. The solid (orange) curve is obtained for α = 3/2, the dashed (brown) curve corresponds to α = 2, and finally the dotted (maroon) curve is the result for α = 3.
where, again, we used the standard thermodynamic relation ∆ = (DF + ST )/V D−1 and we rewrote the free energy F using the first law as (where the choice for the lower bound of the integral comes from the requirement that F (T c ) = 0 at the transition). Note that, for extremely high temperatures, the expression on the right-hand side of eq. (19) vanishes, by asymptotic conformality of the theory. Technically, for very large T the leading-order term in a 1/T expansion of the integral cancels the first term in eq. (19). The contribution that we seek for is the sub-leading term of the integral in 1/T . For large enough temperatures T /Λ 1, one can use the asymptotically conformal result to relate T to r h . This follows from the AdS black-hole expression where the blackness function in eq. (9) becomes: f = 1 − (r h /r) D . As the background turns into the asymptotically AdS BH for large T , substituting this expression in eq. (11) one finds: Using this expression in eq. (17) and in eq. (12), one immediately obtains an approximate expression from eq. (19) as: where we defined the constants: while c 0 is defined in eq. (18). Evaluating the integral in eq. (22) for large T yields: Clearly, this expression is valid only for α = D. The case α = D should be treated separately, and one finds:∆ The conclusion of this semi-analytic calculation is that, as T increases, the normalized interaction measure falls off as a function of T /T c , and the precise shape of the fall-off is a power-law determined by the parameter α. In particular, this power is not directly related to the spacetime dimensionality D, but rather to the nature of the deconfinement transition, which is determined by the value of the exponent α in eq. (4). This is in agreement with known results from lattice computations. In particular: • In D = 3 + 1 dimensions, for pure Yang-Mills theory one expects linear confinement with an asymptotically linear glueball spectrum. Then, eq. (7) determines α = 2 and from eq. (24) one expects∆ ∼ 1/T 2 for some range of temperatures above T c . This has indeed been observed in lattice simulations [18,19,24].
• In D = 2 + 1 dimensions, the confinement-deconfinement transition tends towards being continuous. In particular, for N = 2 and 3 it is known to be a second-order transition, and for N = 4 it may be continuous or a weakly first-order one; a general statement that one can make is that Yang-Mills theories in 2 + 1 dimensions are more inclined to exhibit continuous or weakly first-order transitions than in D = 3 + 1 [38][39][40]. On the other hand, as we discussed above, the nature of the transition in holography is determined by the exponent α. For D = 2 + 1 one expects 14 it to be smaller than 2. Furthermore, in order to have confinement at zero temperature, α should also be larger than 1, see eq. (4). Therefore we expect 1 < α < 2 for the D = 2 + 1 theories under study. Indeed, as discussed in sec. 4, we found that α = 3/2 gives a very good fit to the lattice results.
One should be cautious with the various approximations that we made in this section. First of all, we adopted a semi-analytic approach in determining the holographic background, by setting the scale factor of the metric via eq. (17). In the well-studied case of 3 + 1 dimensions, this approximation works quite well, as shown in ref. [36]. Guided by these results, we made an educated guess for the scale factor of the metric. Another approximation is to treat the relation between T and r h as in the case of AdS, eq. (21). Strictly speaking, this is only valid for large T , where r h is small enough, so that the asymptotically AdS region sets in. However, numerical studies show that this is also a good approximation in a large range of temperatures, ranging all the way from the limit of infinite T down to near T c [35,36]. The final approximation is the large-N limit. Although we believe that the first two approximations can be justified at least in a finite range of temperatures, the validity of the latter cannot be assessed by simple arguments, but should be checked through case-by-case studies. In the case of 3 + 1 dimensions, the SU(N ) Yang-Mills lattice data show that, indeed, the equilibrium thermodynamic quantities per gluon are essentially independent of N [18], and agree well with holographic calculations [35,41]. We hope that this general conclusion also holds for theories in 2 + 1 dimensions. This should be checked by a thorough study of holographic models in D = 2 + 1, that we plan to pursue in the future. In the following sections, after introducing the setup to study SU(N ) Yang-Mills theories in D = 2 + 1 dimensions on the lattice, we compare the lattice results for the equation of state with the predictions of the holographic model that we discussed in this section.

Yang-Mills theories in 2 + 1 dimensions
The continuum formulation of Yang-Mills theories with SU(N ) gauge group in D = 2+1 spacetime dimensions can be defined via the Euclidean functional integral: where g 2 0 (which has the dimensions of an energy) is the bare square gauge coupling, F αβ (x) denotes the non-Abelian field strength tensor, and the functional integration is done over the non-Abelian gauge field A µ (x), taking values in the adjoint representation of the algebra of the gauge group. Like in D = 3 + 1, SU(N ) Yang-Mills theories are also asymptotically free in D = 2 + 1 dimensions; since g 2 0 is dimensionful, perturbative computations for processes at a momentum scale k can be organized as series in powers of the ratio g 2 0 /k [42].
To make the definition in eq. (26) mathematically well-defined at the non-perturbative level, one can introduce a gauge-invariant lattice regularization. The common choice is a regularization on a cubic, isotropic lattice Λ of spacing a, which allows one to trade the continuum field A µ (x) for an at most countable (and finite, if the spacetime is truncated to a finite volume) set of matrices U µ (x), which are defined on the oriented links joining nearest-neighbor lattice sites. The U µ (x) matrices represent parallel transporters on the lattice links, and take values in the adjoint representation of the gauge group. Their dynamics is defined by: where dU α (x) denotes the Haar measure for each U α (x), and S E L is the Wilson gauge action [43]: which tends to S E in the naïve (i.e. tree-level) continuum limit a → 0, with corrections O(a 2 ), and which, being defined in terms of the trace of the plaquette variable: is exactly gauge-invariant at all values of the lattice spacing a.
The expectation value of a physical observable O on the lattice is defined by: This quantity is a ratio of high-, but finite-dimensional, finite, ordinary group integrals, and can be estimated numerically by importance sampling over an ensemble of configurations of link matrices. Previous lattice computations [38][39][40] have shown that, similarly to the D = 3 + 1 case, also D = 2 + 1 non-Abelian gauge theories exhibit linear confinement at low energy (i.e. the interquark potential V (r) grows as σr at large distances r) and a gapped, discrete spectrum of color-singlet glueball states. Furthermore, they also undergo a deconfining phase transition at a finite temperature T c , associated with the spontaneous breakdown of a global center symmetry.
In order to "set the scale" (i.e., to determine the value of the spacing a as a function of β) we used the accurate non-perturbative results reported in ref. [38], which lead to the following expression for the temperature (in units of T c ) as a function of β: .
The statistical and systematic uncertainties in this scale determination are set by those on the value of the critical temperature over the square root of the zero-temperature string tension σ in the continuum limit and in the large-N limit from ref. [38]. The latter reports T c / √ σ = 0.9026 (23) in the N → ∞ limit of this ratio, with a finite-N correction term proportional to N −2 (and valid down to N = 2), whose coefficient is 0.880 (43). So our uncertainty on the temperature scale can be estimated to be of the order of 1%, and does not have a real impact on our analysis (hence, for the sake of clarity, in our plots we do not show the errorbars on the temperature). Another potential source of systematic uncertainties is given by the choice of the physical observable to set the scale: since all numerical simulations are performed at finite values of the spacing a, the determination of the physical scale using different observables can be affected by different discretization artifacts; however, the quantitative effect of the induced systematic uncertainty is small, O(a 2 ). For a comparison with alternative non-perturbative definitions of the scale, see, e.g., refs. [39,40,44]. Finally, note that, in principle, one could also use a perturbative definition of the scale; in particular, for the D = 3 + 1 case high-order perturbative computations in the lattice scheme are available in the literature [45]. However, since in the present work we are interested in a temperature regime where non-perturbative effects are expected to be non-negligible, our determination of the scale is completely non-perturbative.
In our simulations, we generated the gauge configurations using code implementing a 3 + 1 combination of local overrelaxation and heat-bath updates on SU(2) subgroups [46]; for part of our simulations, we also used the Chroma suite [47]. In the following, we denote the cardinality of our configuration ensembles as n conf . Having defined the lattice action via eq. (28), the only parameters that fix the physical setup of our simulations are the number of colors N , the Wilson gauge action parameter β, and the sizes of the lattice along the space-like and time-like directions, which can be expressed in units of the lattice spacing a as aN s and aN t , respectively. As usual in a Euclidean QFT setup (assuming periodic boundary conditions for the bosonic fields), the latter quantity is related to the physical temperature via aN t = T −1 .
The parameters of the simulations that we performed for this work are summarized in table 1; we chose N s N t , which guarantees a good approximation of the thermodynamic limit [48]. In fact, finite-volume effects in the T > T c phase are known to be strongly suppressed, due to the screening phenomenon in the deconfined plasma. Note that, as pointed out in ref. [48] (for the D = 3 + 1 case), the thermodynamics of a gas of free gluons is sensitive to finite-volume corrections, which depend on the product of the linear spatial size of the lattice times the temperature. In the limit of very high temperatures, such corrections lead to quantifiable corrections to ordinary thermodynamic relations. However, at the relatively moderate temperatures probed in the present lattice simulations, the numerical evidence from all previous studies (both in D = 3 + 1 and in D = 2 + 1 dimensions) indicates that screening makes finite-volume corrections essentially negligible for simulations on lattices with N s /N t ≥ 4. In particular, the accurate numerical study of SU(3) thermodynamics in D = 2 + 1 dimensions presented in ref. [39] provided convincing evidence for the strong suppression of finite-volume effects. On the other hand, deviations from the thermodynamic limit in the confining phase are exponentially suppressed by the finiteness of the mass gap: if m 0 denotes the mass of the lightest glueball, then, typically, lattices of linear size aN s ≥ 4/m 0 are such, that systematic effects due to the volume finiteness play an essentially negligible rôle in the lattice computation error budget (which is dominated by finite-cutoff effects, and by statistical uncertainties due to the finite cardinality of the sampled configuration ensemble). The results at zero temperature are obtained from simulations on cubic lattices of volume (aN s ) 3 . In addition to the simulations listed in table 1, for the SU(2) and SU(4) gauge groups we also analyzed the configurations corresponding to the N t = 6 and N t = 8 ensembles (and their respective T = 0 counterparts) taken from ref. [25].
The equation of state of Yang-Mills theories in D = 2 + 1 dimensions can be easily obtained from elementary thermodynamic identities. Let Z(T, V ) denote the partition function for an isotropic system of two-dimensional "volume" V at temperature T ; in the thermodynamic limit V → ∞, the pressure p is related to the free energy F = −T ln Z by pV = −F , and to the trace of the energy-momentum tensor ∆ = T µ µ via: As discussed above, deviations from the thermodynamic-limit relation pV = −F due to the finiteness of the lattice volume can be neglected at the temperatures investigated in this work. Finally, the energy and entropy densities (denoted as and s, respectively) can be obtained from = ∆ + 2p and sT = ∆ + 3p. Our determination of the equation of state on the lattice is done according to the "integral method" [49]: the trace of the energy-momentum tensor is extracted from differences of U 2 T , the expectation value of the average trace of the plaquette at a temperature T : so that the pressure is obtained by integration over β: starting from a lower integration extremum β 0 corresponding to a temperature sufficiently deep in the confined phase. We performed the numerical evaluation of the integral in eq. (34) comparing the trapezoid rule with the method described by eq. (A.4) in ref. [50], which is characterized by systematic errors O(n −4 β ). Since our scan in β values is very fine, the systematic error related to the choice of the numerical integration method has a negligible rôle in the error budget.

Numerical results
In this section, we present our numerical results for the basic equilibrium thermodynamic properties in D = 2 + 1 SU(N ) Yang-Mills theories with N = 2, 3, 4, 5 and 6 colors, and compare them to the predictions from the holographic model introduced in section 2. By virtue of asymptotic freedom, one expects that in the high-temperature limit the thermodynamics of these theories reduces to that of a gas of non-interacting gluons, whose equation of state in the continuum reads:  Table 1: Parameters of the new lattice simulations performed for this work: N denotes number of colors, N t and N s are, respectively, the lattice sizes along the time-like and space-like directions (in units of the lattice spacing). n β denotes the number of β-values (i.e. of temperatures) that were simulated, in each β min ≤ β ≤ β max interval; the T = 0 and finite-T statistics at each β-value are shown in the last two columns. For N > 2, all T = 0 simulations were performed on lattices of size (aN s ) 3 . Our analysis also includes part of the data from the simulations reported in ref. [25].
where ζ(3) 1.20205690316 . . . is Apéry's constant. On the lattice, eq. (35) is affected by cutoff corrections: where the correction factorR I (N t ) can be either estimated numerically or evaluated analytically, order by order in an expansion in powers of N −2 t [51]. For the Wilson action and the integral method that we used in this work, the latter computation yields: (see the appendix A for details). However, it is important to stress that the cutoff artifacts encoded by this correction factor are suppressed at temperatures close to T c , and hence we do not rescale our numerical results byR I (N t ). Since the right-hand side of eq. (35) is proportional to the number of gluon degrees of freedom (one transverse polarization state for each of the N 2 − 1 color d.o.f.), in the deconfined phase it is natural to normalize the dimensionless ratios p/T 3 , ∆/T 3 , /T 3 and s/T 2 obtained in gauge theories with a different number of colors, by dividing them by N 2 − 1. Note, however, that, while this is expected to make the results from different groups collapse onto the same curve in the limit of very high temperatures, in principle there is no obvious reason to expect the same to be true also at moderate temperatures close to T c , where the deconfined plasma is far from being weakly coupled (and the thermodynamics could perhaps be dominated by different degrees of freedom, with unknown scaling properties with N ). Figure 2 shows our results for the pressure per gluon, in units of T 3 , as a function of T /T c . The plot shows the results that we obtained for the various gauge groups, from simulations on lattices with N t = 6 sites in the compactified Euclidean time direction, and for the corresponding spacelike volumes listed in table 1. As we shall discuss in the following, simulating finite-temperature lattices at this value of N t turns out to be an optimal choice, given that it allows one to reach very precise numerical results for the plaquette differences appearing in eq. (34), while keeping the systematic effects due to the finiteness of the lattice cutoff under control. The first, striking result manifest from fig. 2 is the nearly perfect scaling of p with N 2 − 1: in the deconfined phase, the numerical results of the pressure per gluon collapse on the same curve, for all the gauge groups that we simulated (denoted by symbols of different colors). This is analogous to what happens in 3 + 1 dimensions [17,18], and can be clearly contrasted to the behavior in the confining phase, where, on the contrary, all bulk thermodynamic quantities scale proportionally to O(N 0 ), i.e. are independent of N , in the large-N limit [25]. While in principle it is reasonable to expect that the equation of state should be qualitatively similar in all of these theories, the remarkable quantitative agreement among different gauge groups that our results reveal is completely nontrivial. Generally, for SU(N ) gauge theories without quark fields, the leading-order finite-N corrections with respect to the large-N limit are expected to be proportional to N −2 , and hence could amount to relative deviations of the order of 10% for SU (3), or even larger for SU (2). On the contrary, our high-precision lattice results do not reveal any statistically significant evidence 15 of such dependence on N . Figure 2 also reveals that the approach to the continuum Stefan-Boltzmann limit (which is about 0.1913 . . . , out of the vertical axis range) is relatively slow: at temperatures slightly above 7 T c , the pressure is still approximately 15% off from the Stefan-Boltzmann value, indicating that the plasma is still far from an ideal gas of free massless gluons. A qualitatively similar feature is also observed for Yang-Mills theories in 3 + 1 dimensions [18,19] (for which, however, it is important to observe that the energy scale dependence of the physical coupling is different). Finally, figure 2 also shows the prediction (solid orange curve) from the holographic model discussed in section 2, for α = 3/2. Since the holographic prediction for p/[T 3 (N 2 − 1)] is obtained by integration of ∆/[T 4 (N 2 − 1)] over T , based on eq. (32), we fixed the integration constant by imposing consistency of the holographic model with the lattice data at high temperatures. 16 15 The only differences among results corresponding to different gauge groups can be interpreted as statistical fluctuations, and/or in terms of the small systematic uncertainty related to the scale setting, as discussed in section 3 (for the sake of clarity, the horizontal errorbars associated with the accuracy limits on our temperature determination are not displayed in the figures). 16 Note that the approximations involved in the construction of the simple holographic model considered here lead to some deviations from the lattice data for temperatures close to the transition region, and, in particular, to an unphysical non-vanishing value for p/[T 3 (N 2 − 1)] for T → T + c . Since the pressure is a continuous function of the temperature, this would imply a non-vanishing value for p/[T 3 (N 2 − 1)] also for T → T SU (2) SU (3) SU (4) SU (5) SU (6) holographic model, for α = 3 / 2 Trace of the energy-momentum tensor Figure 3: Same as in fig. 2, but for the trace of the energy-momentum tensor per gluon, in units of T 3 . The solid orange curve is the corresponding prediction from the holographic model discussed in section 2, for α = 3/2. Note that, in principle, the holographic model could be refined, to match the lattice data also at temperatures lower than 2.5 T c , through an appropriate choice of the dilaton potential. We postpone a more detailed discussion about this issue to future work. Figure 3 shows our results for the temperature dependence of the trace of the energy-momentum tensor ∆, normalized in units of T 3 and per gluon. The data shown in this plot are the same that we used to evaluate the pressure in fig. 2, hence in the deconfined phase they show the same, approximately perfect, proportionality to N 2 − 1. For this observable, however, one also clearly sees that, in the confining phase (in which the number of physical states is independent of N -except for the special case N = 2: see ref. [25] for a discussion), the results corresponding to different gauge groups do not follow this proportionality law. This effect was not clearly visible in fig. 2, because the pressure and the interaction measure are related to each other by eq. (32), and the signal for p at T < T c is much smaller than at T > T c .
Since p is obtained by numerical integration according to eq. (34), whereas ∆ is directly related to the plaquette expectation values, see eq. (33), it is most natural to compare the lattice results and the predictions of our holographic model for the interaction measure. This is shown by the orange line in fig. 3, which corresponds to the prediction for α = 3/2, as discussed in section 2. As one can see, our holographic model accurately captures the non-trivial temperature dependence of ∆ for all temperatures T 2.5 T c . This good quantitative agreement breaks down in the region of temperatures closer to T c , where the holographic curve falls below the lattice results, and the agreement is only qualitative. While in principle the holographic prediction could be adjusted to fit the lattice data over an even broader temperature range (by including more terms in the dilaton potential), we emphasize that, even for the simple setup discussed here, the model already gives a quantitatively correct description for the high-temperature fall-off of ∆/[T 3 (N 2 − 1)] and, as discussed in section 2, it relates it to the nature of the deconfinement phase transition.
Our lattice results for the other two equilibrium thermodynamic quantities (the energy and the entropy density) are shown in fig. 4 and in fig. 5, respectively. Being linear combinations of the pressure and the interaction measure, these quantities obviously exhibit the same, accurate proportionality to N 2 − 1 as p and ∆, and reveal the same type of mismatch between lattice data and the holographic prediction for temperatures close to T c .
Another interesting problem that we investigated in our simulations is the following: In D = 3 + 1 dimensions, several authors observed that, in the deconfined phase, the trace of the energy-  momentum tensor appears to be proportional to T 2 over a rather broad temperature range [18,19,24]. Since several different interpretations have been proposed for this phenomenon [24], it is interesting to investigate whether a similar effect also occurs in D = 2+1 dimensions. We address this issue in figure 6, by showing our results for the dimensionless ratio ∆/[T 3 (N 2 − 1)], plotted as a function of T c /T : if, in the temperature range under consideration, the interaction measure is dominated by a contribution proportional to T 2 , this should result in a linear behavior in the plot. This is indeed clearly seen in the figure, hence we confirm that, similarly to the D = 3 + 1 case, also in D = 2 + 1 dimensions there is a large temperature interval, starting from the value where ∆/[T D (N 2 − 1)] has its maximum, in which the interaction measure of Yang-Mills theories exhibits a quadratic dependence on T . This figure also shows that, at least in the temperature range T ≥ 2.5 T c , the holographic model captures this type of temperature dependence very well.
The implications of this result are twofold. On the one hand, our finding can be useful to shed light on the nature of the phenomenon in D = 3 + 1. In particular, due to the qualitative differences of Yang-Mills theories in 2 + 1 versus 3 + 1 dimensions, our result might help to rule out some mechanisms that have been proposed to explain the phenomenon in D = 3 + 1, if they are expected to be at work also in the lower-dimensional case. On the other hand, our holographic model leads quite naturally to a power-law decay of ∆/[T D (N 2 − 1)] with the temperature, and, even more interestingly, it suggests a connection between the order of the deconfining phase transition, and the exponent of such power-law decay. As discussed above, the fact that, in general, the deconfinement transition tends towards being more discontinuous when the spacetime dimensionality increases from 2 + 1 to 3 + 1, can thus be directly related to the change from a 1/T to a 1/T 2 fall-off for ∆/[T D (N 2 − 1)].
Finally, we conclude this section with a discussion of the finite-cutoff effects affecting our lattice results. As we mentioned, most of the results presented in this paper are based on finitetemperature simulations using lattices with N t = 6 sites in the compactified Euclidean time direction. One may wonder, whether the corresponding results are close enough to the continuum limit or not. According to our analytical expansion of theR I (N t ) factor in eq. (37), it turns out that, for this value of N t , the lattice Stefan-Boltzmann limit (evaluated with the integral method and the Wilson action on an isotropic cubic lattice) differs from the value in the continuum by approximately 5%: an effect much larger than the statistical uncertainties and the other systematic errors affecting our data. Thus, in principle one may be tempted to rescale all our lattice results by dividing byR I (N t ). SinceR I does not depend on N , this would not change the fact that the thermodynamic quantities are nearly perfectly proportional to the number of gluons, but would lead to slightly different (smaller) numerical values for p, ∆, and s. However, in the temperature region investigated in this study, this naïve rescaling of the results would not be correct: the physical reason is that the distortion of the Stefan-Boltzmann limit encoded byR I is due to modes near the lattice cutoff, and those are not relevant for the physics at temperatures of the order of T c . For this reason, we chose not to rescale our numerical results byR I , but rather to repeat our simulations (except for the computationally most demanding gauge group SU(6)) at the same temperatures and at the same space-like volumes, on finer lattices, with N t = 8. Given that the leading discretization effects of the Wilson action are O(a 2 ), this corresponds to reducing the lattice artifacts by approximately a factor 2. From eq. (34) (in which the plaquette mean values are always O(1), for any N t ) it is also easy to see that, when N t is increased, the difference appearing in the integrand on the right-hand side is affected by a fast decay of the signal-to-noise ratio. As a consequence, it would become extremely difficult to get sufficiently precise results from much finer lattices. Fortunately, however, the discrepancy between our N t = 6 and N t = 8 results for the equilibrium thermodynamic quantities considered in this work appears to be very small, as figure 7 shows. This holds in the whole range of temperatures that we studied, and for all values of N from 2 to 5. Since our N t = 6 and N t = 8 data sets yield compatible results, we did not attempt a continuum extrapolation of the thermodynamic quantities, and we can safely state that, to the level of precision we reached, our results from the N t = 6 lattices are already compatible with the continuum limit.

Conclusions
In this work, we presented a non-perturbative study of the equilibrium thermodynamic properties in the deconfined phase of (non-supersymmetric) SU(N ) Yang-Mills theories in 2 + 1 dimensions, using holographic computations and numerical simulations based on the lattice regularization. This allowed us to combine the advantages of both tools: the former enables one to gain analytical insight on the dynamical properties of these strongly coupled systems, while the latter (once the thermodynamic and continuum limits are taken) provides numerical results obtained from an ab initio approach, directly based on the microscopic definition of the theories for any number of colors, without any assumption or uncontrolled systematic uncertainty. First, we introduced a holographic bottom-up model, inspired by the IHQCD model [28][29][30], which describes the non-trivial dynamics of these strongly interacting non-Abelian gauge theories in the large-N limit. This model reveals a non-trivial relationship between the order of the deconfinement phase transition, and the dependence of the trace of the energy-momentum tensor ∆ on the temperature. In particular, for non-Abelian gauge theories in 2 + 1 dimensions (which, typically, are characterized by a tendency towards a second-order or a weaker first-order transition than in 3+1 dimensions), at temperatures of the order of T c the model favors a behavior approximately compatible with a 1/T decay for the dimensionless ratio ∆/T 3 .
Then, we defined the non-perturbative regularization of SU(N ) Yang-Mills theories on a (2+1)-dimensional Euclidean lattice, and performed a set of high-precision numerical simulations to study their equation of state at T ≥ T c . We compared the results obtained for different numbers of colors, up to N = 6, and found that the trace of the energy-momentum tensor and the related bulk thermodynamic quantities per gluon are independent of N , reflecting a strikingly accurate scaling of the equation of state in the deconfined phase, over the whole temperature range that we probed (up to about 7.5 T c ). This holds for all the gauge groups that we studied, including SU (2), and-at least for these equilibrium thermodynamic observables-supports the potential quantitative relevance of analytical computations relying on the large-N limit (including, in particular, those based on the gauge/gravity correspondence). We also found that, in all the theories that we simulated, ∆ exhibits a clear, characteristic quadratic dependence on the temperature. Both these findings are analogous to those which have been obtained for non-  Abelian gauge theories in 3 + 1 dimensions [18,19,24]. Finally, we compared the lattice results with the prediction of our holographic model, finding good quantitative agreement, at least for temperatures not too close to T c .
In the future, we plan to extend the present study, by investigating the holographic model in more detail, and by looking at different observables, which could be compared with the results of lattice simulations.

Acknowledgements.
We warmly thank J. Engels A Lattice cutoff corrections to the Stefan-Boltzmann limit in 2 + 1 and 3 + 1 dimensions In this appendix, following the calculation in ref. [51], we derive the correction to the Stefan-Boltzmann limit due to cutoff effects on the lattice, for the Wilson discretization of SU(N ) Yang-Mills theory in D = d + 1 dimensions, for d = 2 and 3. Our goal is to evaluate the first few terms of the correction to the Stefan-Boltzmann limit, in an expansion in powers of N −2 t , where N t denotes the number of lattice points in the Euclidean time direction. We take N s , the number of lattice sites along the space-like directions, to be infinite (corresponding to the thermodynamic limit).
Throughout this appendix, we work in lattice units, i.e., we set the lattice spacing a to unity, and denote the pressure as p, the spatial volume as V , and the temperature as T . Moreover, in the following, we use the ∼ = notation to mean equality of two quantities, up to terms which are negligible to the order of precision of our computation.
Notation used throughout this calculation includes: so that: Elementary identities used in this calculation include: as well as: and the expansion: Furthermore, we also use the following finite-sum formula [52]: ln ω 2 + sin 2 πl N t = 2 ln sinh(N t x/2) 2 Nt−1 sinh(x/2) and the limit: implying: Finally, in the following we use Z 1DOF to denote the partition function for one bosonic, massless degree of freedom on the lattice. We concentrate on the integral method for the lattice determination of the pressure [49], in which p(T ), the pressure at a given temperature T , is defined with respect to its value at T = 0. The thermodynamic definition of the pressure reads: and, for an isotropic system, in the thermodynamic limit it reduces to: The pressure can be written as: Thus, taking into account that in the integral method the pressure is defined w.r.t. to its value at T = 0 (obtained as the N t → ∞ limit): Changing variables to y i = N t sin(p i /2), and expanding p i as: so that: Thus the dimensionless ratio p/T 4 can be written as: The last expression can be readily evaluated, using the formulas listed above and integration by parts. In particular, the first two terms are: (10) and: Next, note that, introducing the following polar parametrization for the y i coordinates:    y 1 = y sin θ cos φ y 2 = y sin θ sin φ y 3 = y cos θ and denoting c = cos θ, one gets: which reproduces the expression given in ref. [51], and extends it to the next order in powers of (π/N t ) 2 .

A.2 D = 2 + 1
An analogous calculation can be done for a system in D = 2 + 1 dimensions, for which eq. (A.1) gets replaced by: while eq. (A.2) gets replaced by: Accordingly, the Stefan-Boltzmann limit for the dimensionless ratio p/T 3 , as evaluated on a finite lattice using the integral method, reads: Similarly to the D = 3 + 1 case, each term appearing in the last expression can be evaluated separately. In particular, one easily finds that: 9) and − N 2 − 1 2N 2 t π 2 Next, introducing the following polar parametrization for the y i coordinates: y 1 = y cos θ y 2 = y sin θ , it is easy to prove that: In contrast to the D = 3 + 1 case, the latter expression cannot be rewritten as a simple power series in (π/N t ) 2 with rational coefficients, because, for odd values of x, the π x /ζ(x) ratio is not just an integer (or a rational) number. However, note that, defining: it is possible to write: S + (9) 1.0020083928 . . .