Matrix Models for Deconfinement and Their Perturbative Corrections

Matrix models for the deconfining phase transition in $SU(N)$ gauge theories have been developed in recent years. With a few parameters, these models are able to reproduce the lattice results of the thermodynamic quantities in the semi-quark gluon plasma(QGP) region. They are also used to compute the behavior of the 't Hooft loop and study the exceptional group $G(2)$. In this paper, we review the basic ideas of the construction of these models and propose a new form of the non-ideal corrections in the matrix model. In the semi-QGP region, our new model is in good agreement with the lattice simulations as the previous ones, while in higher temperature region, it reproduces the upward trend of the rescaled trace anomaly as found in lattice which, however, can not be obtained from the previous models. In addition, we discuss the perturbative corrections to the thermal effective potential which could be used to systematically improve the matrix models at high temperatures. In particular, we provide, for the first time, an analytical proof of the relation between the one- and two-loop effective potential: two-loop correction is proportional to the one-loop result, independent of the eigenvalues of the Polyakov loop. This is a very general result, we prove it for all classic groups, including $SU(N)$, $SO(2N+1)$, $SO(2N)$ and $Sp(2N)$.


INTRODUCTION
Over the past thirty years, one of the most important aims in high energy nuclear physics is to explore a new form of matter -the Quark Gluon Plasma (QGP), which has been predicted by the lattice Quantum Chromodynamics (QCD) at finite temperature T . It is believed that under the extremely hot and dense condition created in the heavyion experiments, hadrons dissolve into a gas of almost free quarks and gluons. The understanding of the transition to the QGP phase at high temperature is expected to progress significantly as a result of the Relativistic Heavy Ion Collider at Brookhaven and the Large Hadron Collider at CERN.
In this paper, we will focus on the studies of pure SU (N ) gauge theories in the deconfined phase. At the asymptotic high temperature, one could calculate the thermodynamics from first principles. However, in the semi-QGP region [1,2,3,4], the computation is very challenging. A useful technique is to perform the Monte-Carlo simulations on the lattice which gives insight into parts of the theory inaccessible by other means.
Lattice results for the thermodynamics of pure SU (N ) gauge theories show a rich variety of unexpected behavior near the critical temperature T c . In particular, we are interested in the rescaled conformal anomaly, which is defined as where e(T ) and p(T ) are the energy density and pressure, respectively. The lattice simulations show that the rescaled conformal anomaly is roughly a constant in the temperature region from 1.2T c to 4T c [5,6,7,8]. This demonstrates that∆(T ) is not dominated by a constant "bag pressure". In order to mimic this behavior, one can add a non-ideal term ∼ T 2 in the pressure, as did in the previous matrix models [9,10,11]. Roughly speaking, the dominate contribution to the rescaled conformal anomaly comes from the term ∼ T 2 above 1.2T c . However, what is the fundamental reason for the appearance of the term ∼ T 2 in the correction is still an open question. In principle, one can also consider some other possible forms for the non-ideal term in the pressure. In addition, the latest Wuppertal-Budapest lattice data [8] for pure gauge theory shows that the rescaled conformal anomaly exhibits an upward trend, starting at temperatures above approximately 4T c 1 . Of course, such an behavior can not be described by a term ∼ T 2 in the pressure. The current matrix models are reliable in the whole semi-QGP region, i.e., from T c to about 4T c .
For very high temperature region, in principle, one can use the perturbative approach to compute the thermodynamics. The next-to-next-to-leading order (NNLO) results from hard thermal loop (HTL) resummed perturbation theory agree with the lattice data of the rescaled conformal anomaly at temperatures on the order of 8T c [12]. One can expect that higher order corrections may further improve the agreement and the perturbation prediction works down to even lower temperatures where the matrix model works. However, the calculations would be rather complicated.
Although the commonly used strategy to describe the behavior of∆(T ) in the semi-QGP region is the inclusion of the non-ideal term ∼ T 2 in the pressure, we will propose another possible solution in this paper. As an alternative, our new proposed model is able to describe the thermodynamics in the semi-QGP region as well as the present matrix models. Furthermore, it also predicts the upward trend starting at temperatures about 4T c as found in the lattice simulation. Our results show that the discrepancy of the rescaled conformal anomaly between the new model and lattice data appears around the temperature of 8T c where the perturbation theory becomes reliable.
There is a common feature in the matrix models discussed in Refs. [9,10,11]. All of them use the one-loop thermal effective potential as the ideal contributions [13,14,15]. As a result, the behavior of the matrix model at high temperatures can be systematically improved by adding the perturbative loop corrections. The effective potential is simply the traditional path integral over the gauge fields subject to a constraint [16]. This constraint is that the integration is done while preserving the value of the Polyakov loop at some fixed value. Such a procedure generates a probability distribution for the eigenvalues determined by the fixed value.
The one-loop effective potential has its minima when the value of Polyakov loop equals to ±1 and the pressure calculated from the minimum of the potential equals the known perturbative pressure calculated at vanishing background field. One important goal to compute the loop corrections is to study how the eigenvalue distribution will be affected.
At present, only the two-loop perturbation corrections have been done [17]. The result shows that the eigenvalue distribution doesn't change and the two-loop effective potential is simply a multiplicative and background independent renormalization of the one-loop result [18], where C 2 (A) is quadratic Casimir invariant in the adjoint representation. The above relation was first found for the SU (N ) groups along straight paths from the origin to the degenerate Z(N ) minima [19,20]. These paths run along the edges of the SU (N ) Weyl chamber. However, the minimum of the potential does not exactly follow a straight path as a function of temperature. We should also study if this simple relation still holds inside the Weyl chamber. Aside from SU (N ), it is also important to perform the perturbative two-loop calculations for all other classical gauge groups, including the exceptional group G(2) [11,18].
In Ref. [18], we show that (2) is a very general result for the classic groups including G(2), which holds not only along the edge of the Weyl chamber but also inside. However, by using the MATHEMATICA program, only some specific cases (up to N = 5) of (2) have been verified. For large N, such a computation becomes rather time-consuming. An analytical proof of the simple relation between one-and two-loop effective potential is still needed. In the second part of this paper, we will give an analytical proof of (2) based on the relations of the Bernoulli polynomials. This proof will finally complete the calculation of the two-loop perturbation corrections to the thermal effect potential which has been studied for about twenty years.
The rest of the paper is organized as the following: in Sec. 2, we review the basic idea of the matrix model which was first proposed by Meisinger, Miller and Ogilvie in Ref. [9]. We also briefly discuss the improvements of this model. These improvements make the matrix models more reliable for an quantitative comparison to the lattice data. For more details, we refer the readers to Refs. [9,10,11]. In Sec. 3, we propose a new non-ideal correction in the pressure. Although the form looks very different from those in Refs. [10,11], our new model also predicts the thermodynamics in the semi-QGP region very well. In addition, we also study the upward trend of the rescaled trace anomaly above 4T c with the new model. In Sec. 4, we discuss the two-loop perturbative corrections to the thermal effective potential and provide an analytical proof of the relation given in (2). We prove it for all the classic groups, including SU (N ), SO(2N ), SO(2N + 1) and Sp(2N ). The last section contains the conclusions and outlook.

Review of the Matrix Models
The original matrix model was first proposed by Meisinger, Miller and Ogilvie. Considering the 1-loop effective potential in a constant classic background with a dispersion relation ω k = k 2 + M 2 for the gauge bosons, the model is then given by the first two terms in the high temperature expansion of such an effective potential [9] The classical background field is taken to be diagonal by a gauge rotation and we have q ab is defined as q ab = q a − q b , modulo 1. In this background field the Wilson line is defined as and the Polyakov loop in the fundamental representation is In the confined phase, the potential favors ℓ = 0. While at asymptotic high temperature, ℓ = 1. The lattice simulations show that the value of the Polyakov loop changes near the critical temperature T c . Without the additional mass term in the potential model, the effective potential has its minima when q a = 0(or a Z(N ) equivalent state) which indicates ℓ = 1 at all temperatures. In fact, even including the two-loop perturbative corrections, we still have q a = 0 as the vacuum. In order to model the transition to deconfinement, as shown in (3), we need to include the terms containing the mass scale M which can be treated as the non-ideal corrections to the perturbative contributions.
When the effective potential attains at its minimum, we can identify V as the free energy which is also equal to the minus pressure. By requiring the effective potential is stationary with respect to the eigenvalues, one can determine the values of q a 's at a given temperature. The only parameter M in the potential model in (3) is fixed by requiring the phase transition happens at the critical temperature T c . 2 As a result, there is no free parameter in this model. Therefore, it is also called the 0−parameter matrix model.
The predicted thermodynamics can be found in Ref. [9]. However, for a quantitative agreement between the model and lattice data, it needs to be improved further. In addition, such a model will have a negative pressure near T c which is not physical.
In order to get a better fit to the lattice data, in Ref. [10], an improved matrix model was proposed. The ideal contribution is unchanged as compared to (3). Using the notations in Ref. [10], we denote it as V pt , The corresponding non-ideal corrections take the following form Under the conditions that c 2 = 0 and c 3 = c 1 , the above potential is restored to the one given by (3) with c 1 = 15M 2 /(4π 2 T 2 c ). In the improved model, the ratio between c 1 and c 3 is kept unfixed which is equivalent to adding a term that is independent of the background field and simply proportional to ∼ T 2 . All the parameters c 1 , c 2 and c 3 are temperature independent constants. In order to avoid the pressure being negative near the critical temperature, its value at T c is forced to be zero in the improved model. As a result, there is only one free parameter left which can be determined by fitting the lattice data of the pressure.
Furthermore, to get the correct latent heat at the critical temperature, a "bag" constant term was introduced in the matrix model [11]. Namely, we can replace the parameter c 3 in (8) with c 3 (∞) + (c 3 (1) − c 3 (∞))/t 2 where the variable t is equal to T /T c by definition. The improved models 3 are in good agreement with the lattice data and have been used in other studies [21,22,23,24,25,26].

A New Approach to Improve the Matrix Model
The original model in (3) is very simple which contains no free parameter and can be used to describe the thermodynamics for only a qualitative purpose. Therefore, one should consider the improvements for a quantitative description of the lattice data. As we discussed above, the models with a new term ∼ T 2 q 2 ab (1−|q ab |) 2 significantly improve the results. In this section, we propose another different approach to improve the original matrix model which also leads to a perfect fit to the lattice data in the semi-QGP region.
In Refs. [10,11], the study of the thermodynamics is focused on a temperature region from T c to 3T c . However, as mentioned in the introduction, the rescaled trace anomaly has an upward trend starting at some higher temperature. The previous matrix models fail to reproduce this behavior. With our new model, we will also study the rescaled trace anomaly at higher temperatures.
As before, the ideal contribution in the new model is still given by (7). However, instead of the inclusion of the term ∼ T 2 q 2 ab (1 − |q ab |) 2 in the non-ideal contributions, we consider the temperature dependence of the mass scale M in the original model. Although the above mentioned models treat M as a constant, its temperature dependence was already considered in the quasi-particle models [27,28,29]. Of course, the exact form of the T -dependence near the critical temperature can not be calculated from first principles. In this work, we assume the following form for the mass scale 4 where g(T, c ′ 2 ) is the running coupling which will be considered at the one-loop level in the following, Here, c ′ 2 is another constant parameter that needs to be fixed. With this assumption, the non-ideal terms in the original matrix model in (3) becomes Instead of ∼ T 2 , the temperature dependence ∼ g 2 (T, c ′ 2 )T 3 appears in the above equation. It turns out to be responsible for the explanation of the upward trend of the rescaled trace anomaly above 4T c . Similar as the ∼ T 2 terms in the previous matrix models, these non-ideal contributions can be considered as the high temperature corrections to the ideal contributions which are at the order of ∼ T 4 .
As in Ref. [11], to get the correct behavior at the critical temperature, we should also include two terms that are independent on the background field. Finally, the non-ideal terms in our new model can be expressed as Here, for convenience, we define c ′′ 1 ≡ 4π 2 c ′ 1 /15. The four parameters c ′ i with i = 1, 2, 3, 4 are temperature independent constants. By comparing with (8), we find that (12) can be obtained with the following replacements, The procedure to fix all the parameters follows the same idea as in the previous work [10,11]. To proceed further, we parameterize the eigenvalues of the Polyakov loop q a 's under the straight line ansatz: which assumes the eigenvalues have constant spacing and automatically satisfies the constraint q 1 + q 2 + · · · + q N = 1 for SU (N ) gauge group. In the above equation, 0 ≤ s ≤ 1. The confining vacuum corresponds to s = 1 while the perturbative vacuum is at s = 0. One can also check that with this parametrization, the value of Polykov loop is real. In fact, we can always consider a real-valued Polyakov loop under a global Z(N ) rotation.
Notice that for N > 3, the exact solution for the q a 's doesn't satisfy this straight line ansatz. However, the deviation from the straight line is very small [11]. For simplicity, we adopt the ansatz (14) in this section.
Then we can perform the sums in the potential and redefine the potential as where Here, r ≡ 1 − s. The background field r can be obtained by the variational approach and the solution corresponds to the minimum in the deconfined phase is given by In fact, the results in (16) and (18) can be directly read off from Ref. [11] by using the replacements given in (13). Impose the condition W(0, 1) = W(r 0 (1), 1) = 0 at the critical temperature, one can get the following relations among these parameters An additional constraint on these parameters can be obtained by considering the latent heat which is the jump in the energy density at T c . By construction, both the pressure and the energy density vanish at the confined phase, as a result, the latent heat equals to the rescaled energy density e(T c )/T 4 c and so the trace anomaly (e(T c ) − 3p(T c ))/T 4 c . Based on the matrix model, the straightforward calculation of the latent heat leads to another constraint on the parameters c ′ 2 and c ′ 4 . It reads  Table 1: The values of the four parameters in the matrix model and the corresponding "bag" constants for SU (N ) gauge theories for N =3,4 and 6.
where L(N ) is the latent heat renormalized by the number of perturbative gluons N 2 −1.
According the the lattice result [7], we have L(3) = 0.209, L(4) = 0.287 and L(6) = 0.342. Now there is only one free parameter left which will be fixed by fitting the lattice data of the pressure [7]. In order to make an direct comparison with the previous model, we also consider the SU (N ) gauge theories for N = 3, 4 and 6. The values of the four parameters and the corresponding "bag" constants B = π 2 (N 2 − 1)c ′ 4 T 4 c /45 are summarized in table 1. 5 The "bag" constant B increases with N and the numerical values have a good agreement with those from the matrix model in Ref. [11].  (3) thermodynamics obtained from the lattice simulation and matrix models. We consider the dimensionless pressure p/T 4 , energy density e/(3T 4 ), as well as one third the rescaled trace anomaly defined in (1). All quantities are also scaled by 1/8. Left: the matrix model is taken from Ref. [11]. Right: the matrix model is taken from (15).
The comparisons of the SU (N ) thermodynamics obtained from lattice simulation and matrix models for N = 3, 4 and 6 are shown in Figs.(1), (2) and (3), respectively. Thermodynamic quantities considered here are the dimensionless pressure p/T 4 , energy density e/(3T 4 ), as well as one third the rescaled trace anomaly defined in (1). Two kinds of the matrix models are included. The left parts of these figures show the results predicted by the matrix model in Ref. [11], while the right parts show the results from the model given in (15).
Although having very different forms, both matrix models appear to reproduce the lattice data reasonably well in the semi-QGP region. In the new model, the complication of the inclusion of a T -dependent mass scale is compensated by the absence of the term  (1). All quantities are also scaled by 1/15. Left: the matrix model is taken from Ref. [11]. Right: the matrix model is taken from (15).  Figure 3: Comparison of the SU (6) thermodynamics obtained from the lattice simulation and matrix models. We consider the dimensionless pressure p/T 4 , energy density e/(3T 4 ), as well as one third the rescaled trace anomaly defined in (1). All quantities are also scaled by 1/35. Left: the matrix model is taken from Ref. [11]. Right: the matrix model is taken from (15).
Taking a close look into the figures, we find that for the two-parameter matrix model, deviations from the lattice data become visible at higher temperatures, while for the new model, the agreement is rather good in the entire semi-QGP region. Despite the structures of the non-ideal terms in these models, such a difference is probably induced by the way how one fixes the free parameter. Instead of fitting the lattice data of the pressure as we did in this paper, in Ref. [11], the authors fixed this parameter by considering the behavior of the rescaled trace anomaly.
As already mentioned in the introdudction,∆(T ) is almost constant in the temperature region from 1.2T c to 4T c where the background fields are approximately zero 6 , therefore, the contribution to∆(T ) is dominated by the term ∼ T 2 in the potential 7 . Using the value of∆(T ) from the lattice simulation, the parameter related to the term ∼ T 2 , which is denoted as c 3 (∞) in Ref. [11], can be determined. This determination is universal for any N since the numerical result of∆(T ) is basically independent on N as indicated by the lattice.
However, the above analysis can not be used in the model considered in this paper. As in the previous models, we also find a rapid fall-off of the background fields near T c and their values almost vanish above 1.2T c . However, beside the ∼ T 2 term, the term ∼ T 3 g 2 (T, c ′ 2 ) also contributes to the rescaled trace anomaly. The same behavior of ∆(T ) is also found in the new model because the contribution to∆(T ) coming from the term ∼ T 3 g 2 (T, c ′ 2 ) is also approximately a constant from 1.2T c to 4T c . However, as the temperature increases, such a contribution starts to increase with T . The same trend is found in the lattice simulation [8]. On the other hand, with previous matrix models, the change of∆(T ) is negligible above 1.2T c .
In Fig.(4), we show the thermodynamic quantities for SU (3) up to 8T c . The corresponding lattice data can be found in Ref. [8]. With our new model, the upward trend of the resacled trace anomaly above 4T c is observed which, however, disappears when using the previous matrix models. Quantitatively, a slight discrepancy exists because the values of the parameters are taken from table 1 and we don't attempt to refix them. The lattice data of the resacled trace anomaly is available for even higher temperatures where, however, the model description fails. As compared to the lattice simulation, a slower increasing with temperature is found from our model. We comment that the temperature dependence of the mass scale in (9) is proposed only at the phenomenological level and our consideration is concentrated in the non-perturbative region. Above 8T c , the HTL perturbation theory works very well for the SU (N ) gauge theory.
It is also interesting to see the temperature dependence of the gluon mass M (T ) as given in (9). If we parameterize the mass M (T ) as M 2 (T ) = g 2 ef f T 2 N/6 according to the leading order perturbative result, the behavior of the effective coupling g 2 is given in Fig.(5). On the right hand side of this figure, we plot the ratio of M 2 (T ) to T 2 which is the effective coupling scaled by a factor N/6. It turns out that the effective coupling increases as the temperature approaches to T c and the ratio of M 2 (T ) to T 2 is essentially independent on N . Some visible N -dependence  Lastly, we mention that due to the rapid fall-off of the background fields near the critical temperature, the Polyakov loop predicted by the two-paramter matrix model differs from 1 only in a very narrow region, from T c to 1.2T c . However, the renormalized Polyakov loop from the lattice varies over the entire semi-QGP. Although the increasing rate of the loop near T c found in the new model is relatively smaller as compared to the two-paramter matrix model, the agreement with the lattice results is still not satisfactory.

Two-loop Perturbation Corrections to the Effective Potential
In this section, we discuss the perturbation corrections to the effective potential. The basic properties of these loop corrections are of particular interest as we already mentioned in the introduction. Here, we will focus on the studies of the simple relation between one-and two-loop effective potential as given in (2). Notice that the explicit result of the effective potential up to two-loop is known since long, however, there is nothing in the way one performs the computation that suggests such a simple relation. Previous work can be found in [17,18,19,20], however, an analytical proof is not available so far.
Using the explicit form of the effective potential, we will give an analytical proof of (2) in the following. For completeness, we will prove it for all the classic groups, including SU (N ), SO(2N ), SO(2N + 1) and Sp(2N ).
We rewrite the effective potential up to two-loop order as In We also divide the two-loop effective potential into two parts. Γ (2) = Γ (2) f + Γ (2) i . Γ (2) f denotes the contributions from usual two-loop free-energy diagrams, while Γ (2) i corresponds to the contribution from the insertion diagram [17,18]. In Γ (2) f , the indices a,b and c run over both diagonal and off-diagonal indices 8 . In Γ (2) i , each structure constant contains the diagonal indices d, while b and c denote off-diagonal indices. In addition, we have E −a ≡ (E a ) † for the off-diagonal generators which defines the "minus" indices.
In (2), the Casimir invariant is related to the rank of the group d(r) through For each classic group, we have that C 2 (A) = N − 1 for SO(2N ), N − 1 2 for SO(2N + 1), N + 1 for Sp(2N ), and C 2 (A) = N for the SU (N ) groups.
In order to continue our discussion, some generalities on the classical Lie algebras should be mentioned 9 . For any semi-simple Lie algebra, the commutation relations in the Cartan basis are given by The diagonal generators in the fundamental representation are the components of H which are the orthonormal matrices spanning the Cartan subalgebra. The off-diagonal generators are the orthonormal E α , labelled by the roots α. They are vectors in Cartan space. A root labeled by an off-diagonal index a has d components which are the structure constants f d,a,−a . They can be determined by the commutation relation There is another kind of structure constants f α,β,−α−β which connect off-diagonal generators. These generators are normalized as The proof of (2) can be achieved by the following two steps. Firstly, by using the corresponding commutation relations of the group generators, we are able to get some simplified expressions for the two-loop effective potential which contain only the Bernoulli polynomials while the structure constants appear in (22) and (23) have been get rid of. Secondly, based on a set of relations of the Bernoulli polynomials, the simplified two-loop results are expressed in terms of B 4 (x). Then a direct comparison between the one-and two-loop effective potential proves the renormalization relation given in (2).
In fact, the insertion diagram involves sums over diagonal indices d which can be performed quite easily as they correspond to inner products between the corresponding roots. This has already been done in [18]. For later use, the simplified results for Γ (2) i for each group are list below.
For SU (N ), we have For the orthogonal groups, the result is the following and Γ (2) i (SO(2N + 1)) = Γ (2) i (SO(2N )) + 2g 2 i,j ( B 1 (q i + q j ) + B 1 (q i − q j )) Finally, for Sp(2N ),the simplified insertion diagram reads Notice that in these results, the constraints on the indices i = j and i = l apply.

Proof for SU(N)
Starting from (22) and (23), we are able to calculate the two-loop perturbative correction to the effective potential. First of all, we need to know the structure constants. They can be obtained from the commutation relations of the generators in Cartan basis. For SU (N ), there are N (N −1) off-diagonal generators E ij ≡ λ ij with i, j = 1, · · · , N and i = j. The explicit forms are given by (λ ij ) kl = δ ki δ jl / √ 2. In addition, we have N − 1 traceless diagonal generators H d ≡ λ d with d = 1, · · · , N − 1, The commutators between diagonal generators are obviously zero. The nonvanishing commutators we need are 10 Here, λ d ii is the ith diagonal component of λ d . From (33) we find that the roots α ij = ( λ ii − λ jj ). As mentioned before, the dth component of α ij is the structure constant f d,ij,ji . In addition, (34) indicates that the non-vanishing structure constant is For SU (N ), the arguments in the Bernoulli functions have two different cases: for a diagonal index a, q a = 0; for an off-diagonal index a = ij, q a = q i − q j .
Firstly, we consider (22) in the case where only one of the indices a, b and c is the diagonal index 11 . By using the property of the roots α ij = 1 √ 2 ( e i − e j ) with an orthonormal basis { e i } spanning an N -dimensional space, we can get rid of the structure constants and write this kind of contribution as The other contribution from (22) then becomes which corresponds to the case where all the indices a, b and c are the off-diagonal indices.
Using (35), (37) can be simplified as Notice that i = j = k means i = j, i = k and j = k.
For later use, we split Γ (2) i given in (28) into two parts as The next step is to rewrite the simplified two-loop effective potential in (36), (38), (39) and (40) in terms of B 4 (x). Firstly, we consider Γ (2) f | I and Γ (2) i | I . These two terms are easy to deal with because there is only one argument q i − q j in the Bernoulli polynomials.
Using the definition of the Bernoulli polynomials, we can show the following two identities Using (41), we can express Γ (2) f | I as Γ (2) Furthermore, by using (42), the sum of Γ (2) f | I and Γ (2) i | I becomes It is actually the final result we need for Γ (2) | I . However, for Γ (2) f | II and Γ (2) i | II , we have two different arguments q i − q j and q i − q k in the Bernoulli polynomials simultaneously. To complete this proof, it is very important to use the following non-trivial identity (46) The proof of (45) is given in Appendix B. After this identity is proved, it is straightforward to get the following result In the third and second lines of (47), we have used the following two equations, respectively.
i =j =k These equations can be obtained by using the fact that i = j = k indicates the three indices i, j and k are equivalent which can be interchanged under the summation. Combine (44) and (47) and use the expression of the one-loop effective potential we finally prove the simple relation between one-and two-loop effective potential for SU (N ) as given in (2).

Proof for SO(2N) and SO(2N + 1)
For both groups, there are N (2N − 2) off-diagonal generators E ηi.η ′ j , where In the above equation, i, j = 1, · · · , N and i > j. We define the indices i with an associated sign η. Similarly, j is defined with η ′ . The signs η or η ′ are independently ±1.
For SO(2N + 1) there are 2N additional off-diagonal generators For either of the groups the N -dimensional Cartan subalgebra is spanned by mutually commuting and orthogonal generators H d , with The structure constants can be obtained from the commutation relations 12 From (54) and (55), the roots can be expressed as From (56) and (57), the non-vanishing structure constants are given by 12 For SO(2N ) and SO (2N + 1), if a typical off-diagonal index b is denoted as In (57), with our notation, E ρk.ηi should be understood as −E ηi.ρk if i > k. Similarly for E η ′ j.ρ ′ l .
For these two groups, the arguments in the Bernoulli functions have three different cases: for a diagonal index a, q a = 0; for an off-diagonal index a = ηi, q a = ηq i and for an off-diagonal index a = ηi.η ′ j, q a = ηq i + η ′ q j .
We start with SO(2N ) and split the two-loop effective potential into two parts. For Γ (2) f , we have Γ (2) In order to get (63), we should mention that there are six terms which contribute to the final result and they correspond to six possible ways to get a non-vanishing structure constant. If we symbolically write the structure constant in (61) as f AB,CD,EF , then the six possible ways are the following 13 : In fact, every term has exactly the same form as that in (63) expect for the constraint in the summation. The constraint for an individual term is just one possible way to order the values of i, j and k, for example i > j > k. As a result, the total contribution is given by (63) with the constraint i = j = k.
Similarly as SU (N ), Γ f | I corresponds to the case where only one index is diagonal in (22) while Γ (2) f | II is the part with all the indices being off-diagonal in (22). To get these two equations, we have used (59) and (61) to get rid of the structure constants. To keep the notations compact, we define B n (q i ± q j ) ≡ B n (q i + q j ) + B n (q i − q j ) with n = 1, 2, 3, 4.
For the Γ In the second line of (65), we use the fact that For the sum of Γ (2) i | I and Γ (2) f | I , by using (41) and (42), it is easy to show that Γ (2) i | I + Γ (2) Before we discuss Γ (2) i | II and Γ (2) f | II , we need to prove the following identity (69) Based on the identity given in (45), we can easily prove (68) by using the following relation Then the sum of Γ (2) i | II and Γ (2) f | II can be obtained by summing over the indices i, j and k in (68) with the condition i = j = k.
Adding up (67) and (71) and using the expression of the one-loop effective potential we complete the proof of (2) for SO(2N ). Once SO(2N ) is done, the proof for SO(2N +1) is easy. We only need to consider the extra terms in the effective potential which arise due to the 2N off-diagonal generators E ηi . Using (58) and (60), we can get Combing (73) and (75) and using (41) and (42), we can show Γ (2) In (68), by setting q k = 0, the sum of Γ (2) f | II and Γ (2) i | II then has the following form The one-loop effective potential for SO(2N + 1) is given by Adding up (67), (71), (77) and (78), we prove (2) for SO(2N + 1).

Proof for Sp(2N)
The diagonal generators of Sp(2N ) is Here, σ i are the Pauli matrices with i = 1, 2, 3. The N − 1 matrices λ d are the same as for SU (N ), and we need the additional λ N = 1 The corresponding off-diagonal generators E ij are In addition, we have additional N (N + 1) off-diagonal generators denoted as E ηij ; the first index η is a sign index. They are defined by where σ η = 1 2 (σ 1 + iησ 2 ).
The structure constants can be obtained from the commutation relations 14 Like for SU (N ), we can write the roots for Sp(2N ) in terms of the orthonormal basis { e i }, Here, the first roots are associated with the generators E ηij when i > j and the second roots are associated with the generators E ij . The roots α ηi (which can be also written as α ηii ) come from E ηii . For Sp(2N ), we found that the arguments of the Bernoulli functions are the following: for a diagonal index, q d = 0; for an off-diagonal index, we have q ij = q i − q j and q ηij = η(q i + q j ).
Using the results given in (91), the contribution with one diagonal index in the structure constant in (22) is given by For other contributions from (22), we have In the above equation, the first line is obtained by using (86)  Here, the subscript indicates the constraint on the indices for each structure constant. It is easy to check that each contribution related to the individual structure constant as given in the above equation has the same form and sum of the six terms corresponds to the constraint i = j = k in (93). Similarly, we can get the following non-vanishing structure constants from Eqs. (88) and (89) Using this explicit result, the last line of (93) is obtained. The contribution from the insertion diagram given in (96) has to be rewritten in a proper way in order to combine with the corresponding terms in Γ (2) f . As before, we split it into two parts Similarly as what we did in (65), terms which have zero contribution have been dropped in Γ (2) i . At this point, we can easily show that the combination of Γ (2) f | I and Γ (2) i | I can be expressed in terms of B 4 (x) as For the contributions from Γ (2) f | II and Γ (2) i | II , we can rewrite the result in a proper way in order to make use of (68), The first summation in the above equation is the same as SO(2N ) and the result is given in (71). In order to perform the second summation, we can use (68) in the case q i = q k which gives Here, we have used (41) and (42) and some vanishing terms are dropped according to (111). Then sum over indices i and j with i = j in (99), we find the result for the second summation in (98) reads Notice that the third line in (99) vanishes under the summation. Now we can write down the final result for (98) (101) Together with (97) and the one-loop effective potential for Sp(2N ) we prove (2) for Sp(2N ).

Conclusions and Outlook
We reviewed the basic ideas of the construction of the matrix model which were first proposed by Meisinger, Miller and Ogilvie. The non-ideal contributions in this model have been further improved for a quantitative fit to the lattice simulations on the thermodynamics for SU (N ) gauge theories. Previous work included the Bernoulli Polynomial B 4 (x) in the non-ideal contributions for such a purpose which reproduced the lattice result in the semi-QGP region very well. In this work, we consider the temperature dependence of the mass scale appears in the original matrix model, and the corresponding new model also does a good job in the semi-QGP region, therefore, can be treated as an alternative to the previous models. On the other hand, with the exact temperature dependence given in (9), our new model is able to get the upward trend of the rescaled trace anomaly as found by the recent lattice simulations. Starting at about 4T c , this quantity increases with temperature and such a behavior can not be obtained by using previous matrix models. Up to 8T c , the new model is in good agreement with lattice results. Beyond this temperature, the HTL perturbation theory works well.
In addition, we discussed the perturbative corrections to the one-loop thermal effective potential which is used as the ideal contributions for all the matrix models. In particular, there is a simple relation between the one-and two-loop effective potential as shown in (2). An analytical proof of this relation is first given in this paper. We show it is quite general and holds for all the classic groups.
One could include the two-loop contribution in the matrix models which is expected to improve the high temperature behavior. If we use the same non-ideal form as given in (12) and assume the two-loop matrix model is obtained from (15) by including an overall factor (1 − 5N g 2 /(16π 2 )), an optimal fit shows that the values of the parameters are essentially unchanged as compared to table 1 while the coupling g 2 is extremely small. The same can be found when using the two-parameter model in Ref. [11]. It turns out in our approach, the two-loop corrections are very small and negligible.
On the other hand, for a realistic value of g 2 , one has to refit these parameters in the model and our results show a faster increase of the rescaled trace anomaly above 8T c which indicates a better agreement with the lattice data as compared to the one-loop model. However, the behavior in the semi-QGP region is not satisfactory. It suggests that in order to reproduce the thermodynamics in both semi-QGP and perturbative region, the proper form of the non-ideal contributions that can accommodate the twoloop corrections should be considered. Furthermore, the two-loop corrections do not change the distribution of the eigenvalues of the Polyakov loop, one has to add the nonideal terms to generate a phase transition. Therefore, the calculation at three-loop order is also important which is expected to answer if or not the distribution of eigenvalues can be modified. We postpone these studies to our future work.
In d = 4 dimensions and for k = 0, 1, we have the following four Bernoulli polynomials: with The above expressions are defined on the interval 0 ≤ x ≤ 1 and they are periodic functions of x, with period 1. For arbitrary values of x, the argument of the above Bernoulli polynomials should be understood as x − [x] with [x] the largest integer less than or equal to x, which is nothing but the modulo function. If −1 ≤ x ≤ 1 we can drop the modulo functions and the Bernoulli polynomials reduce to where ǫ(x) is the sign function.

B Proof of Equation (45)
In this appendix, we prove the identity given in (45). To make the proof more clear, we can use the periodicity of the Bernoulli polynomials to impose a constraint on the values of the background field q i . Without losing any generality, we can assume 0 ≤ q i < 1 thanks to the following properties of the modulo operation Here, n i is an integer which satisfies 0 ≤ q i + n i < 1. n j is defined similarly. With our constraint of the values of the background field, the argument of the Bernoulli polynomials is in the interval (−1, 1). In principle, we can directly use (108) to prove (45). However, there is a problem related to B 1 (x) because it has discontinuities at integer x. The value of B 1 (x) at x = 0 depends on the way how x approaches zero, from above or from below 15 . We have B 1 (0 + ) = − 1 2 and B 1 (0 − ) = 1 2 . Due to the discontinuity, q i = q j should be understood as q i = q j ±ǫ. Here, ǫ is an infinitely small (positive) number and the + sign corresponds to q i approaches q j from above, while − sign is for approach from below. This problem of discontinuity seems to make the discussion very complicated if some of the background fields are equal. Fortunately, we find that, in order to prove (45), there is no need to specify a way how q i approaches q j , although the value of F(q i , q j , q k ) depends on this specification 16 .
On the other hand, with the condition q i = q j = q k , it is obvious that we will not encounter the problem of discontinuity as we discussed above. we can directly use (108) to prove (45) by assuming a specific ordering of these background fields, for example, q i < q j < q k . 17 Let's consider the second case where q i = q j = q k . Although B 1 (q i − q j ) is not determined until one specifies a way how q i approaches q j , the product B 1 (q i −q j ) B 3 (q i − q k ) is zero independent on the specification. This is because B 3 (0) = 0. Actually, it is easy to check that both the left and the right sides of (45) are T 4 48 in this case. For the same reason, we don't need a specification on how q i approaches q j in (42).
The last case corresponds to only two of the three background fields are equal. For example, we consider q i = q j . The following discussion simply applies for q i = q k and q k = q j .
Explicitly, we have 15 With our constraint on the background field, the argument x of the Bernoulli polynomials satisfies −1 < x < 1. So we only need to consider the discontinuity at x = 0 for B1(x). 16 Of course, (45) does not depend on this specification as we will see later. 17 The only reason to use a specific ordering of these background fields is to make the sign functions explicit, then the proof becomes straightforward. There are six possible ways to order the values of qi, qj and q k . Of course, the result does not depend on the orderings.
Using the identity we find that the first term in F(q i , q j , q k ) and F(q j , q k , q i ) are cancelled by each other. Therefore, the left side of (45) becomes F(q i , q j , q k ) + F(q j , q k , q i ) + F(q k , q i , q j ) = 4 B 1 (q k − q i ) B 3 (q k − q i ) + 2 B 2 (0) Furthermore, the above result can be rewritten in terms of B 4 (q k − q i ) with the help of (41) and (42), Summing up the above discussions, we have proved (45). This identity is true for any value of the background fields and independent on the way how q i approaches q j in case these two fields are equal.