Thermodynamics and phase diagrams of Polyakov-loop extended chiral models

We study the thermodynamics and phase diagrams of two-flavor quantum chromodynamics using the Polyakov-loop extended quark-meson (PQM) model and the Pisarski-Skokov chiral matrix ($\chi M$) model. At temperatures up to $T\approx2T_c$ and baryon chemical potentials up to $\mu_B=\SI{400}{\mega\electronvolt}$, both models show reasonable agreement with the pressure, energy density, and interaction measure as calculated on the lattice. The Polyakov loop is found to rise significantly faster with temperature in models than on the lattice. In the low-temperature and high baryon density regime, the two models predict different states of matter; The PQM model predicts a confined and chirally restored phase, while the $\chi M$ model predicts a deconfined and chirally restored phase. At finite isospin density and zero baryon density, the onset of pion condensation at $T=0$ is at $\mu_I={1\over2}m_{\pi}$, and the transition is second order at all temperatures. The transition temperature for pion condensation coincides with that of the chiral transition for values of the isospin chemical potential larger than approximately \SI{110}{\mega\electronvolt}. In the $\chi M$ model they also coincide with the transition temperature for deconfinement. The results are in good overall agreement with recent lattice simulations of the $\mu_I$--$T$ phase diagram.


I. INTRODUCTION
The first phase diagram of QCD appeared in the 1970s, and at the time it was thought that it consists of two phases: A hadronic low-temperature phase and a high-temperature phase of deconfined quarks and gluons. Today, the conjectured phase diagram in the µ B -T plane is far more complicated. In particular, it is believed that the deconfined quark phase at high density and low temperature consists of various colorsuperconducting phases, with different patterns of spontaneous symmetry breaking. Some of these phase may even be inhomogeneous, see Refs. [2-4] for reviews. However, only a few exact results are known: due to asymptotic freedom we know that at asymptotically high temperature QCD is in a plasma phase of weakly interacting quark and gluons. Similarly, due to asymptotic freedom and the existence of an attractive interaction via one-gluon exchange, we have a superconducting colorflavor locked phase at asymptotically high densities. A severe problem in the efforts to map out the QCD phase diagram between these asymptotic regions is that one cannot use lattice simulations at finite baryon density due to the so-called sign problem. The fermion determinant is complex and one cannot apply standard Monte-Carlo techniques based on importance sampling. Thus one typically has to resort to low-energy effective models * afolkest@mit.edu † andersen@tf.phys.ntnu.no such as the quark-meson (QM) model, the Nambu-Jona-Lasinio (NJL) model, and their Polyakov-loop extended versions [5].
There are other external parameters that one can introduce in addition to the temperature and the baryon chemical potential. For example, one can add a (strong) magnetic background field B. QCD in strong magnetic fields is relevant in e.g. heavy-ion collisions [6][7][8] and compact stars [9]. One can also use an independent chemical potential µ f for each quark flavor. In twoflavor QCD, this implies that one uses µ u and µ d , or equivalently, µ B in addition to the isospin chemical potential µ I . Isospin asymmetry and the possibility of Bose condensation of charged pions may also be relevant to compact stars. An advantage of QCD in a magnetic field or at finite isospin density (but at zero µ B ) is that there is no sign problem, and one can therefore use standard Monte-Carlo techniques to study the phase diagram of these systems. This opens up for the possibility to confront results from model calculations with those of the first-principle method of lattice QCD.
In this paper, we study the thermodynamics of twoflavor QCD using the Polyakov-loop extended quarkmeson (PQM) model and the Pisarski-Skokov chiral matrix model (χM ) [1] adapted for two flavors in the meanfield approximation. In this approximation, the mesonic fields are treated at tree level while the fermion fields are integrated over in the Gaussian approximation. At one loop, dropping the mesonic fluctuations is equivalent to working in the large-N c limit. Sometimes the no-sea approximation is made, which simply means that one discards the fermionic quantum fluctuations. How-ever, one should keep vacuum fluctuations since there is no apriori reason to omit them, not even in low-energy models of an underlying theory. Secondly, it turns out that the inclusion of quantum fluctuations change the order of a phase transition in some cases; in the two-flavor QM model the chiral transition changes from first order to second in the chiral limit, showing the importance of keeping them. Moreover, in almost all mean-field calculations to date, the parameters of the Lagrangian are determined at tree level. This is inconsistent since the effective potential has been determined in the one-loop large-N c approximation. The parameters should always be determined at the same level of accuracy as the effective potential, otherwise erroneous results may occur. For example, the onset of pion condensation at T = 0 takes place when the isospin chemical potential equals half the pion mass. This exact result is only reproduced if the matching of the parameters is done in a consistent manner.
The paper is organized as follows. In Sec. II, we discuss various aspects of the Polyakov loop and its properties. In Sec. III, the gluonic sector of the PQM and χM models is reviewed, while in Sec. IV their chiral sectors are discussed. Complications related to minimizing the effective potential at nonzero baryon chemical potential is discussed in Sec. V. In Sec. VI, we present the main result of the paper, namely the thermodynamic functions and the phase diagram in the µ-T and µ I -T planes. Our results are compared with lattice results. In Sec. VII, we summarize and conclude. Four appendices are devoted to technical details.

II. CENTER SYMMETRY AND THE POLYAKOV LOOP
Let T a be the generators of SU (N c ) in the fundamental representation. A gauge transformation of the QCD gluon field A µ = A a µ T a is of the form for any Ω(x) in the fundamental representation of SU (N c ). This transformation leaves the gluonic Lagrangian invariant, and is thus a symmetry of the action of the pure gauge theory. However, when studying QCD at finite temperature T = β −1 in the imaginary time formalism, choosing a generic Ω(x, τ ) ruins the periodicity of A µ (x, τ ) in imaginary time τ , as required for the field configurations summed over in the partition function Z. 1 Restricting ourselves to transformations that satisfy avoids the problem, but there is a larger group of symmetries that preserves the imaginary time periodicity of A µ . Consider instead a generic gauge transformation that satisfies for some G(x, τ ) ∈ SU (N c ). Let A µ be the transformed field. We then get If G(x, τ ) is constant in space and imaginary time and commutes with A µ for all (x, τ ), then the gauge field is periodic. Since Ω(x, τ ) is a matrix in the fundamental representation of SU (N c ), which is irreducible, G is proportional to the identity matrix by Schur's lemma. Let G = λI Nc , where I Nc is the N c × N c identity matrix and λ ∈ C. Since we know that G ∈ SU (N c ), we have that λ = λ n is one of the N c -th roots of unity, and all possible matrices G are given by G n = λ n I Nc = e −2πin/Nc I Nc , n = 0, . . . , N c − 1 . (5) Clearly {G n } forms a finite group that is isomorphic to Z Nc , and it is the center group of SU (N c ). We refer to aperiodic gauge transformations, characterized by G n = I Nc , as twisted gauge transformations or center transformations.
In pure gauge theory, the expectation of the Polyakov loop operator is an order parameter for deconfinement [10][11][12]. For QCD with dynamical quarks, it is an approximate order parameter, similar to the chiral condensate. The thermal Wilson line is given by where T τ denotes time ordering. Here A a 4 is the Euclidean temporal gauge field that replaces the Minkowski temporal gauge field in the Euclidean QCD action through the replacement A 0 → iA 4 , with A 4 Hermitian [13].
The Wilson line is not invariant under (periodic) gauge transformations Ω(x, τ ), but transforms as L(x) → Ω(x, β)L(x)Ω † (x, 0). Taking the trace over color indices, however, yields a gauge-invariant operator, which is the definition of the Polyakov loop operator Under twisted gauge transformations, the Polyakov loop operator transforms nontrivially, . Thus, we see that the Polyakov loop is gauge invariant (n = 0), but not center symmetric. Therefore the thermal expectation value of the Polyakov loop operator transforms as Thus if Φ = 0, the center symmetry is spontaneously broken. While Φ is related to the free energy of a heavy quark, the conjugate Polyakov loopΦ is the analogue of Φ for antiquarks, and it is obtained from Φ by replacing T a → −(T a ) * = −(T a ) T , ie. by switching to the complex conjugate representation. This is shown in Appendix D. After using that tr X = tr X T for any matrix X, we findΦ if the fields A a 4 (x, τ ) are real. More generally, if A a 4 (x, τ ) are complex, we find whereT τ denotes anti-time ordering. In full QCD one would not have to worry about complex A a 4 fields, since any particular A a 4 configuration occurring in the path integral should be real. Thus, when averaging over all field configurations in full QCD one has Φ † = Φ . However, the A a 4 background fields we will deal with in the effective models should be thought of as mean fields A a 4 , which in QCD with µ B = 0 are obtained by averaging real field configurations with potentially complex weights e −S E . The distinction between Φ † andΦ matters only when dealing with Φ( A 4 ) andΦ( A 4 ), i.e. the loops evaluated at the mean field A 4 , rather than the expectation values of the loops themselves. In summary:Φ( A 4 ) = Φ † ( A 4 ) even though Φ = Φ † . In the effective models, the former occurs. We will return to this in the following when we discuss minimization of the effective potential at µ B = 0. Using the expression forΦ is anyway always correct. If we use Φ † in place of Φ without an expectation value, it is implied that A a 4 is real.
In the original paper by McLerran and Svetitsky [10], it was argued that e −βF (x1,...,xn,y1,...,yñ) = (11) where F is the color-averaged free energy of a configuration of quarks located at x 1 , . . . x n and antiquarks located at y 1 , . . . yñ. We can thus interpret −T ln Φ(0) as the free energy of a single quark and −T ln Φ † (x) as the free energy of a single antiquark. If Φ = 0, this implies that the free energy of a quark is infinite, or that quarks are confined.
Another way to think of confinement is in terms of the quark propagator. The Polyakov loop is proportional to the expectation value of the traced propagator of a heavy quark analytically continued to imaginary time.
In Appendix D, we show that (12) where V is the volume and with no sum over the color index a. One can take the vanishing of the propagator and therefore the Polyakov loop as a sign of confinement.
In the context of the PQM and χM models, it is convenient to choose a gauge which simplifies the Polyakov loop as much as possible. The Weyl gauge A 4 = 0 would make the Polyakov loop trivial; however this gauge is not compatible with the periodicity requirement of the gauge field in the imaginary time formalism [14,15]. Instead one can choose the so-called static gauge [16], where Furthermore, one can rotate the gauge fields so that A 4 is in the Cartan subalgebra of the Lie algebra of SU (N c ) [13]. In the case of N c = 3, the gauge field in the Polyakov gauge can be written as where λ 3 and λ 8 are the two diagonal Gell-Mann matrices. Defining we can express the background gauge field as The thermal Wilson line can then be written as (17) Taking the trace to obtain the Polyakov loop and its conjugate yields When A 4 is constant in space, and thus also r and q, we see that at the classical level. Thus we conclude that a state with q = 1, r = 0 is a deconfined state, while a state with q = 0, r = 0 is a confined state.

III. GLUONIC SECTOR
In this section, we discuss the gluonic sector of the effective(/grand) potential Ω of the PQM and χM models, which is the main difference between the two models. They are somewhat different, but they both involve a phenomenological pure-glue potential with a few parameters that are determined such that several physical quantities from pure glue lattice simulations are reproduced.

A. PQM model
It is known from lattice simulations that a first-order phase transition, corresponding to gluonic deconfinement, happens at T 0 = 270 MeV in pure SU (3) gauge theory [21]. A first-order transition is what is expected on the basis of universality, as argued by Svetitsky and Yaffe in Refs. [11,12]. In addition to the knowledge of the location of the phase transition, various thermodynamical properties such as the pressure and internal energy as function of temperature have been established [22][23][24]. Finally, one also has simulations of the value of Polyakov loop as function of temperature [25].
With the knowledge of for example T 0 , P (T ) and Φ(T ) from lattice simulations, one can write down a phenomenological potential U glue (Φ,Φ, T ) that reproduces these three quantities. The first requirement necessitates that the form of the effective potential admits a first-order transition in the first place. For the second criterion, we can find the pressure from the effective potential as P = −U glue (Φ(T ),Φ(T ), T ) with Φ andΦ evaluated at the minimum of U glue . Regarding the form of the potential, several things can be said on general grounds. It must be symmetric under center transformations since the gluonic action is center symmetric. Remembering the transformation rule for Φ andΦ under center transformations, we see that the potential can be a function the terms ΦΦ, Φ 3 andΦ 3 only. Additionally, there is no reason for any asymmetry between Φ andΦ in a pure gluonic system, and we thus require that the potential is symmetric under Φ ↔Φ. Finally, we should demand that the minimum of U glue at low temperatures is at Φ =Φ = 0, while at high temperatures it should equal or asymptotically approach Φ =Φ = 1.
Several potentials have been suggested in the literature [26][27][28][29], and some of the more frequently used are compared in Ref. [30]. The number of fit parameters vary from two [27] to seven [28]. One of the models by Ratti, Rößner, Thaler and Weise [29], which is the one we use for the PQM model, takes the form with the temperature-dependent coefficients We take all the parameters except T 0 as given in the original paper, meaning We note that the potential presented above does not take into account the backreaction of quarks onto the gluonic sector. However, from the fact that the running coupling in QCD depends on the number of quark flavors, as evident from the one-loop expression it is natural to let T 0 = T 0 (N f ), so that the behavior in the gluonic sector also depends on N f , since g determines the strength of the interactions between the gauge fields. The authors of Ref. [31] parametrize this dependence as where the constants areT = 1.77 GeV and α 0 = 0.304. This expression is heuristically obtained by assuming that the temperature dependence of g is governed by (25) with Λ = T and that the deconfinement phase transition occurs at a specific coupling, so that we can solve

B. Chiral matrix model
The gluonic part of the chiral matrix model was developed as an effective model for pure SU (3) gauge theory in Refs. [32,33], where the degrees of freedom are r and q only. The potential consists of two terms; one is obtained by integrating out a fluctuating gauge field in a background gauge field A 4 to one loop. The other term models a nonperturbative contribution and is added by hand.
The one-loop perturbative contribution to the effective potential of pure SU (3) Yang-Mills theory reads where The second term in Eq. (28) is the Weiss potential, first calculated by Weiss [14,34] and Gross, Pisarski, and Yaffe [35]. To drive the system to confinement at low temperatures one adds a phenomenological potential, which is chosen to be of the form where with four fit parameters c 1 , c 2 , c 3 and T d . At temperatures below roughly T ∼ T d the ∼ T 2 term will dominate over the ∼ T 4 perturbative term, and it can thus drive the system to confinement with an appropriate choice of the fit parameters. The T 2 behavior is chosen since it has been observed in lattice data that the subleading contribution to the pressure goes as ∼ T 2 [36,37]. The parameters c 1 and c 3 are chosen so that the pressure in the confined phase of the pure gauge theory is zero and so that a phase transition happens at T d . The former is an approximation, but it is reasonable since the pressure of the confined phase in SU (3) gauge theory is very low compared to the deconfined phase, as lattice data show [22]. Furthermore, we choose T d = 270 MeV, which is roughly the deconfinement temperature in SU (3) gauge theory [22]. Then only c 2 remains as a fit parameter. It is determined by fitting the interaction measure (E − 3P )/T 4 predicted by the full gluonic potential, to lattice data, and the result is c 2 = 0.830 [1], which gives In Fig. 1, we show contour plots of the perturbative (left panel) and the nonperturbative (right panel) contributions to the gluonic potential in the chiral matrix model. The perturbative contribution V pt (q, r) has minima at q = r = 0, q = 0 and r = 1 etc., and maxima at q = 1, r = 0 and q = r = 1 2 etc. The potential reflects the center symmetry of SU (3). The nonperturbative potential V nonpt (q, r) behaves the opposite way; it has minima where the perturbative potential has maxima and vice versa. This potential also reflects the center symmetry. The full gluonic potential is then the sum of the two contributions, and the latter has been constructed so that the two terms are competing. At high temperatures V pt dominates while for low temperatures V nonpt dominates.

IV. CHIRAL SECTOR
In this section, we will discuss the chiral sector of the two models. In Ref.
[1], Pisarski and Skokov add a phenomenological quark term to the χM model which is not present in the QM model, but apart from this the chiral sector in the χM model corresponds to the QM model. Adding this phenomenological term to the PQM as well, their chiral sectors are identical. Note however that while Ref.
[1] includes the strange quark, we treat only the two light quark flavors.

A. Quark-meson model
To obtain the two-flavor quark-meson model we couple two N c -plets of fermionic fields via Yukawa interactions to the linear sigma model with an approximate SU (2) L × SU (2) R symmetry. The fields ψ 1 and ψ 2 are taken to represent up and down quarks, respectively. The Lagrangian of the two-flavor quark-meson model in Minkowski space is where ψ is the flavor doublet and π is the isospin triplet (π 1 , π 2 , π 3 ) T , with π ± = (π 1 ± iπ 2 )/ √ 2.σ is a scalar isospin singlet that will attain a vacuum expectation value that corresponds to the chiral condensate. The τ i are the Pauli matrices acting in flavor space and γ 5 ≡ iγ 0 γ 1 γ 2 γ 3 . µ I is the isospin chemical potential and µ = 1 3 µ B the quark chemical potential.
Let us identify the symmetries of the QM model. In the chiral limit, meaning h = 0, the QM Lagrangian has a global SU (N c ) × SU (2) L × SU (2) R × U (1) B symmetry, while at the physical point (h = 0) the symmetry is The U (1) B symmetry gives rise to conservation of baryon number, with the associated baryon chemical potential µ B = 3 2 (µ u + µ d ) = 3µ, where µ u and µ d are the up and down quark chemical potentials, respectively. The isospin chemical potential µ I is given by The full chiral sector of the PQM and χM models is obtained by coupling the quarks to a temporal gauge field in the Euclidean Lagrangian, whereγ 4 = γ 0 andγ i = −iγ i are the Euclidean gamma matrices, δ µν the Euclidean metric, and µ ∈ {1, 2, 3, 4}. We take m 2 < 0 so thatσ attains a vacuum expectation value v, and defineσ = v + σ. Here v is the chiral condensate. To obtain the effective potential U chiral , we work to one-loop order and neglect bosonic fluctuations. As mentioned, the latter approximation is equivalent to taking the large-N c limit. The contribution to the thermodynamic potential from the Lagrangian (35) coupled to the background gauge field then consists of two terms; a vacuum term arising from the tree-level mesonic potential and the fermion determinant, and a thermal piece coming from the same fermion determinant.
For µ I = 0 and with ∆ ≡ gv, we write where we have, to one-loop order after renormalization and consistent parameter fixing, Here m q , m π and m σ are the physical quark, pion and sigma masses at T = 0, respectively, while f π is the pion decay constant. The quantity ∆ is, in addition to the rescaled chiral condensate, the constituent quark mass, and satisfies ∆(T = µ = 0) = m q . The derivation of Eq. (38) can be found in Appendices A and B, and the definition of the functions F (m 2 ) and F (m 2 ) are given in Eqs. (C5) and (C6). Equation (39) is obtained in the same way as one would calculate the free fermion partition function, except one now has a complex effective chemical potentialμ j that differs for each quark color j, withμ with [A 4 ] jj the j-th diagonal element. This derivation requires using the Polyakov gauge, where A 4 is diagonal.
The thermal quark potential U q,T as function of q and r at µ = 0, T = 100 MeV, and ∆ = 300 MeV is shown in Fig. 2. We see that this potential, whose qualitative shape is mostly unchanged for other values of ∆ or T , drives (q, r) towards q = r = 0, which corresponds to deconfinement. Thus, it is expected that the addition of quarks to the gluonic potential lowers the deconfinement temperature.
We finally note that for µ = 0, Eq. (39) can become complex. This will be discussed in Sec. V.

B. A phenomenological quark term
In addition to the (partly) phenomenological gluonic sector, Pisarski and Skokov add to the χM model a phe- nomenological quark term. In the two-flavor case, it is given by where m cur is the current quark mass. This term is added in order to achieve that ∆ → m cur in the hightemperature limit. Let us show how this works: In the high-temperature limit we expect q = r = 0, and we can thus set the Polyakov loop to be Φ =Φ = 1. Let us furthermore assume that µ = 0 and T ∆. We can Thus, to leading order we find As we assume high temperatures, we consider the potential only up to subleading temperature dependence ∼ T 2 . Furthermore, we assume that ∆ is small, which we expect in the high-temperature phase where chiral symmetry is approximately restored. Using this we keep only leading and subleading terms in ∆. Thus, in the high temperature limit we find that the effective potential goes as (44) Minimizing this potential with respect to ∆, we immediately find and we expect ∆ → m cur in the high-temperature limit. When we later in this section investigate the thermodynamics of the PQM and χM models, we will assess the effects of U q,cur on the thermodynamic functions. Due to its ad hoc nature it should preferably affect the thermodynamics minimally while still achieving its purpose of ensuring the quark mass to approach the current quark mass in the approximately chirally restored phase.
There is one more problem we must face. For the case of µ = 0, the quark effective potential is real for any q, r ∈ R since the two terms in Eq. (39) are complex conjugates of each other. However, upon introducing of µ = 0 this breaks down, and the potential becomes complex in general.
The solution suggested in Refs. [1, 38,39] is, when µ = 0, to let the background A 4 field become non-Hermitian by setting q ∈ R and r = iR with R ∈ R. As discussed in Sec. II, this is not as unreasonable as it first seems, since A 4 in the PQM and χM models represents the mean field of a quantum field, A 4 = A qu 4 , and when µ = 0 in full QCD the Euclidean action becomes complex. Because of this, even if each field configuration is real, when carrying out the path integral where we weight each field configuration with e −S E , we might get that A qu 4 is complex. Inserting r = iR into Eqs. (18) and (19), we find which gives that both are real, but with different values. For the RRTW potential, which is a function of Φ andΦ, it is clear that the Polyakov-loop potential becomes real, and thus the full potential is also real for all (q, R) ∈ R 2 . For the χM model this is also the case, since the only potentially complex terms in Eqs. (29) and (32) are and However, the problem is that the full effective potential Ω χM is unbounded as a function of R for low temperatures, and that Ω PQM is unbounded as a function of R for all temperatures. The part of the potential that depends on R in the χM model,Ṽ = U q,T + U χM , is shown in Fig. 3 as function of (q, R). We see that there is no minimum for |R| < 1. This behavior persists for any R. In Refs. [38,39] the authors deal with this by suggesting to choose the physically realized state to be the lowest saddle point of Ω(q, R). This recipe gives Φ =Φ with both being real, thus giving a real effective potential. That we obtain Φ =Φ is desirable, since we do not expect the free energy of single a quark to be equal to that of a single antiquark when µ = 0. However, choosing a saddle point is arbitrary and does not follow from any known principle of thermodynamics. It also is pointed out in Ref. [40] that interface tensions cannot be calculated within this scheme.
An alternative approach, which is the one used in the following, is to keep q, r ∈ R and minimize Re Ω under the constraint Im Ω = 0. If we interpret a complex Ω as signaling an unstable state, this might be reasonable. It turns out that a global minimum of Re Ω can always be found at r = 0, and with r = 0 we always have Im Ω = 0. This means that we can set r = 0 and minimize Ω freely with respect to q only. However, this scheme gives Φ =Φ ∈ R, which is not what we expect from the quark/antiquark free energy interpretation of Φ andΦ. However, the equality of the two can be seen as a result of the fact that we are doing a mean-field treatment of A 4 instead of a mean-field treatment of the actual Polyakov loops. The quantities we are calling Φ andΦ in the PQM model are, in the Polyakov gauge, which are not equivalent to the expectation value of the Polyakov loop quantum operators. Thus, the free energy interpretation should not be taken too seriously.

VI. THERMODYNAMICS
We now have all the ingredients needed to investigate the thermodynamics of the PQM and χM models at one loop in the large-N c limit. For a given temperature and chemical potential we numerically solve with r = 0 and require that we have a global minimum, where the full effective potential Ω for the two models reads We also add the term U q,cur to the PQM model, for the same reason that it is added to the χM model. The parameters P 0,χM and P 0,PQM are constants that we subtract from the effective potential so that the condition is satisfied for each of the two models. This constant will turn out to be small and has a negligible effect on the thermodynamics. However, it makes thermodynamic quantities divided by T 4 better behaved at temperatures close to zero. Once ∆ and q are determined as functions of T and µ, we can determine Ω as a function of T and µ only. We can then calculate the pressure P , quark density n q = N /V , energy density E and interaction measure I = (E − 3P ) as functions of µ and T via the relations To determine the one-loop couplings we use the following values for the masses and the pion decay constant which yields the parameters which are the one-loop values of the running couplings in the MS scheme at the renormalization scale This is the scale that is consistent with σ =0 in the on-shell scheme. The sigma particle is a broad resonance whose mass is usually taken to be in the 400 MeV to 800 MeV range, and the most recent estimated mass range is 400 MeV to 550 MeV [41]. We have chosen a value of 500 MeV, but we vary it to gauge the sensitivity of our results.

A. Order parameters
In Fig. 4, we show the order parameters ∆(T ) ∆(T =0) and Φ(T ). We point out that ∆ does not go to zero at high temperatures, but rather approaches ∆ ≈ m cur , as expected from the discussion in Sec. IV B. We find that the χM model reaches full deconfinement at T ≈ 250MeV, while the PQM model reaches Φ = 1 more slowly, and it is in a "semi-deconfined" state between roughly 200 and 400 MeV. Like in Ref.
[1] we find that the Polyakov loop in both models rises faster than on the lattice, which can be seen from Fig. 6 where the models are compared to lattice data from Refs. [42,46]. One can define the pseudo-critical temperature for example by the temperature at which the order parameter has dropped to half its zero-temperature value. Another definition is the temperature at which the derivative of the order parameter has its peak. In this paper, we will stick to the latter. In Fig. 5, we show the derivatives of the order parameters ∆(T )/∆(T = 0) and Φ(T ). From the figure, we see that the pseudocritical temperatures for the chiral and deconfinement transitions coincide for both models, with the inflection points of ∆ being located at The uncertainty is given by varying σ from 400 MeV to 550 MeV, with the lowest sigma mass corresponding to the lowest T c and vice versa.

B. Pressure, energy density and interaction measure
Let us now turn to the thermodynamic functions. Figure 7 shows the comparison between the Stefan-Boltzmann (SB) normalized pressure as calculated on the lattice [43][44][45] and in the two chiral models, with each data set plotted against T /T c for its respective T c (this also applies to plots in the following). For the (2 + 1)-flavor lattice data we normalize with T c = 155 MeV [43,46], while for the model data the T c s are given by (67). The two-flavor lattice data from Ref. [45] are obtained directly as function of T /T c without knowledge of T c . The pressure in the model data and (2 + 1)-flavor lattice data is normalized with, while the two-flavor lattice data, which are not continuum extrapolated, are normalized with the relevant Stefan-Boltzmann pressure for a discretized space-time (see Ref. [45] for details). The uncertainty bands in the HotQCD data correspond to uncertainty in the continuum extrapolation. The uncertainty bands in the models are obtained by varying the sigma mass within the uncertainty range given in Ref [41], which as mentioned is 400 MeV to 550 MeV. The lowest m σ corresponds to the lowest temperature, and vice versa. Both the PQM and χM models show reasonable agreement with lattice data above T = T c , although with a slightly lower pressure. Below and around T = T c the PQM model appears to have a pressure that is significantly lower than what lattice data show. However, below and around T c we expect mesons to exist and contribute to the pressure, and by neglecting mesonic fluctuations in the model, we have underestimated the pressure. Since the pions have masses of ∼ 140 MeV below T c while the quarks have masses of ∼ 300 MeV, we expect that the mesons would provide a significant contribution to the pressure in this range. For temperatures below T c , the agreement with the χM model is worse, since there is a small but nonzero pressure causing P/P SB to blow up for low temperatures due to the T 4 -dependence of P SB . However, this does not mean that the pressure diverges or that it is large. It only means that a small non-zero pressure exists for T > 0. This pressure is insignificant, as we see if we compare the pressure of the χM and the RRTW models without SB-normalizing, which is done in Fig. 8.
The energy density E and the interaction measure I, both normalized with E , are shown in Fig. 9 (with error bars are obtained in the same way as for the pressure). We find fairly good agreement between the PQM model and two-flavor lattice data up to T ∼ 1.5T c . The peak of the interaction measure in the χM model is shifted to higher values than what is seen in the PQM model and two-flavor lattice data. The χM model also has an interaction measure that is negative for low temperatures and a peak that is too low.
In Fig. 10 we plot the Boltzmann-normalized version of the quantities at µ B = 3µ = 200 MeV and µ B = 400 MeV. The model data are compared with (2 + 1)-flavor lattice data from Ref. [47]. We compare with lattice data where the chemical potentials for the light flavors are µ L = 3µ u = 3µ d = 200 MeV and µ L = 400 MeV, while the strange chemical potential is chosen so that the net strangeness density is zero. Since µ s = µ d = µ u , we use the notation from Ref. [47] and denote 3µ d = 3µ q as µ L instead of µ B in the (2 + 1)-flavor simulation. We compare with lattice data where net strangeness density is zero, since this scenario should resemble the two-flavor situation more than when µ u = µ d = µ s , for which the strangeness is nonzero.
We see that the pressure, which increases as function of temperature at nonzero baryon chemical potential, agrees fairly well with lattice data for both models. For the interaction measure, we find that the PQM model has a significantly higher peak than lattice data at µ B = 400 MeV, while the χM model is in better agreement.

C. Phase diagram
In this subsection we present various phase diagrams of the two models. We first discuss the different phases in the µ-T plane. Then we move on to the phase diagram in the µ I -T plane, where we include the possibility of condensation of charged pions. In calculating the phase diagrams, we have dropped the U q,cur term, since its effect on thermodynamics and critical temperatures is found to be entirely negligible. Figures 11 and 12 display the phase diagrams for the two models, where the pseudocritical temperatures corresponding to the inflection points of ∆ and Φ are indicated, in addition to the temperature where Φ = 1 2 . We see that the chiral and deconfinement phase transitions happen roughly simultaneously also for nonzero chemical potentials. Note however that referring to the inflection point of the Polyakov loop as "deconfinement" in the regime of high chemical potential is mis- leading. It is correct that the chiral symmetry in the models is approximately restored above the orange lines in Figs. 11 and 12, since we can see from Figs. 13 and 14 that ∆ → 0 quickly for temperatures higher than the crossover temperature T c . However, it is not correct to assume that quarks are deconfined everywhere outside the phase boundaries, since the inflection point of the Polyakov loop can be relatively far away from the region where center symmetry is approximately restored (Φ ≈ 1). This is visible from Figs. 15 and 16, which show the value of the Polyakov loop in the µ−T plane. We see in the PQM model that the Polyakov loop is close to the confining value of Φ = 0 also for µ > 300MeV, given low temperatures. Interestingly, we see that we approach deconfinement in the χM model in the high-density limit, which is not the case in the PQM model. This is a major difference between the two models. However, it is hard to assess which behavior best reflects QCD due to the lack of lattice data in that region. We also note that at sufficiently large temperatures, some of the crossovers become first order phase transitions, with the transition from crossover to first order marked by a critical point. In the χM model the critical points of of the two transition lines coincide, while for the PQM model only the line of chiral transition has a critical point. This is another qualitative difference between the two models.

D. Pion condensation
We now move on to discuss the phase diagram in the µ I -T plane, which requires a treatment of Bose condensation of charged pions. For simplicity, we set the baryon chemical potential to zero in the remainder of this section.
In addition to an expectation value ofσ we now allow for a nonzero expectation value of π 2 1 + π 2 2 denoted by π 0 . Introducing ρ = gπ 0 in analogy with ∆, the treelevel potential can be written as The quark energies can be read off from the zeros of the quark determinant and read where we have defined The effective potential at T = 0 then is where the last term is the one-loop contribution. It cannot be evaluated analytically for nonzero ρ. In Ref. [49], it was evaluated by isolating the divergent pieces and writing U vac = V div + V fin . The divergent term V div was then evaluated using dimensional regularization, and the poles in (evaluating integrals in d = 3 − 2 dimensions) were removed by renormalization using the MS scheme in the usual way. The running parameters were finally eliminated in favor of the physical masses and the piondecay constant. The final result is where the finite contribution V fin is which must be evaluated numerically. Eq. (76) reduces to Eq. (38) for ρ = 0; this can be easily seen by noting that V fin = 0 in this case. The medium-dependent part of the one-loop effective potential at µ B = 0 is Note that this term vanishes at T = 0. As discussed previously, we see that two and two terms are complex conjugates of each other, and U q,T is thus real, reflecting that there is no sign problem when µ B = 0.
In Ref. [49], it was shown tha the zero-temperature effective potential (76) exhibits a second-order phase transition at exactly µ c I = 1 2 m π . This was done by expanding U vac in powers of ρ, U vac = α 0 + α 2 ρ 2 + α 4 ρ 4 evaluating it at ∆ = m q (i.e. in the vacuum). The critical chemical potential is defined by α 2 = 0. The transition is second order at µ I = µ c I since α 4 was found to be positive for this value of the isospin chemical potential. Figures 17 and 18 show the solutions ∆(T, µ I ) and ρ(T, µ I ) for the χM model (note the different axis orientations in the two plots). These plots are similar for the PQM model. We clearly see that no pion condensation occurs for low µ I .  [52]. Furthermore, to the right of this line we have µ I > ∆, and the energies for the u andd quarks (72) and (73) are no longer minimized by |p| = 0, but rather by |p| = µ 2 I − ∆ 2 . This can been seen as a signal of a transition to a BCS state [50][51][52]. However, we do not have a thermodynamic phase transition, since the same O(2) symmetry is broken on both sides of the dashed line. In Refs. [53][54][55] the phase diagram in the µ I -T plane is mapped out using lattice methods for 2+1 flavors. The phase diagram in Fig. 19 is in especially good agreement with that obtained from the lattice: The chiral and deconfinement transition lines coincide for small values of the chemical potential and meet the BEC transition line at (µ meet I , T meet ). For chemical potentials larger than µ meet I , the BEC and chiral lines coincide. Finally, the deconfinement line penetrates smoothly into the BEC phase. The athors of Refs. [53][54][55] identify this line inside the O(2)-broken phase as the BEC-BCS transition line. Again, the same O(2) symmetry is broken on either side of that line.

VII. SUMMARY
In this paper, we extended the chiral matrix model of Pisarski and Skokov [1] to finite baryon and isospin chemical potential. For temperatures up to approximately up to 2T c and baryon chemical potentials up to µ B = 400 MeV, this model and the PQM model show reasonable agreement with lattice results for a number of thermodynamic functions. However, the Polyakov loop rises faster with temperature than on the lattice. A significant difference between the models was found in the deconfinement phase diagram. In the χM model the deconfinement transition also goes from a crossover to a first order transition, with the critical point located at the same point as the critical point for the chiral transition. This is not the case in the PQM model, where the deconfinement transition is a crossover for all µ. Furthermore, the chiral matrix model predicts deconfinement in the low-T , large-µ regime, while the PQM model predicts a quarkyonic phase. Thus, the two models predict different phases of matter in the low-temperature, high-density regime, which is the most significant difference between the two models.
Regarding pion condensation at finite temperature, the two models predict essentially the same phase diagram; the only difference is that the deconfinement transition merges with the other lines at large chemical potentials in the χM model, while in the PQM model the deconfinement line penetrates into the BEC/BCS phase. The phase diagram is in overall good agreement with the lattice results of [53][54][55].

VIII. ACKNOWLEDGMENTS
The authors would like to thank R. D. Pisarski for helpful discussions and the CP-PACS collaboration for providing their lattice data.
After renormalization, we find the one-loop potential where the subscript is a reminder that the renormalized parameters are in the MS scheme. In Appendix ??, we show how one can relate the running parameters in the MS scheme to the parameters in the OS scheme and hence the physical masses and the pion-decay constant. Substituting the parameters (B23)-(B26) into Eq. (A6). we obtain Eq. (38).

Appendix B: Parameter fixing
In this Appendix, we discuss the parameter fixing in the quark-meson model using the on-shell scheme. This was first done in Ref. [56]. At tree level, the parameters of the Lagrangian can be expressed in terms of the phsyical masses and the pion decay constant as Beyond tree level, these parameters become running parameters in the MS scheme and the relations (B1)-(B4) no longer hold. The counterterms in this scheme are chosen such that they exactly cancel the ultraviolet divergences coming from the loops. In the on-shell scheme, the counterterms are chosen such that they exactly cancel the loop corrections that appear in the calculations and the parameters therefore still satisfy the above treelevel relations and are not running. Using that the bare parameters in the two renormalization schemes are the same, we can relate the corresponding renormalized parameters.
The first renormalization condition we impose is that σ = 0, i.e. that the loop correction to the one-point function vanishes and that the minimum of the renormalized effective potential coincides with that of the classical mesonic potential. The classical one-point function is denoted by Γ (1) = it = i(h − m 2 π v) and the classical minimum is then given by the equation of motion t = 0. Let δΓ (1) be the one-loop large-N c correction to the one-point function. The renormalization condition σ = 0 is then The first on-shell renormalization condition on the two-point function is that the counterterms exactly cancel the loop corrections that have not been eliminated by the renormalization condition σ = 0. This gives the mass counterterms where the four-dimensional integrals A(m 2 ) and B(m 2 ) have been defined in Eqs. (C2) and (C3). In the onshell scheme one also takes as a renormalization condition that the residue of the propagator at the pole mass equals unity. This implies dΣ σ,π (P 2 ) dP 2 where Z σ,π is the wavefunction renormalization counterterm. One finds Let us now return to the renormalization condition (B5), which reads The relation t = (h − m 2 π )v implies upon variation a relation among the counterterms, In order to find δh OS , we need to compute δv 2 OS . The one-loop correction to the quark-pion vertex is of order N 0 c and so is the one-loop correction to the quark field, implying Z ψ = 1. Consequently, √ Z π Z g 2 g 2 = 1, or δg 2 g 2 + δZ π = 0. A similar argument now applies to m q = gv; since the quark mass correction at one-loop is of order N 0 c , we find δgv + gδv = 0. Combining these relations, we can write Eq. (B12) as . (B13) We finally use Eqs. (B1)-(B2) to find relations among the corresponding counterterms This yields The bare parameters in the Lagrangian are independent of the renormalization scheme. This implies the following relations among the renormalized parameters in the two schemes where we have used that m 2 = m 2 OS etc. The counterterms in the on-shell scheme consist of a pole in plus finite terms. The former is exactly the counterterm in the MS scheme. Moreover, the parameters in the on-shell scheme are expressed in terms of the physical masses and the pion decay constant. We can then express the running parameters in the MS scheme as We note from Eqs. (B25) and (B27) that the product ∆ = gv = g MS v MS , i.e. it does not run with the renormalization scale.
Substituting Eqs. (B23)-(B27) into the effective potential (A6), we obtain Eq. (38). We have emphasized the importance of matching the parameters in the oneloop large-N c approximation for consistency. For example, the onset of pion condensation at T = 0 takes place only at µ I = 1 2 m π if the parameters are determined in the same approximation as the effective potential itself. Moreover, to show the effects of renormalization, we show in Fig. 21 the one-loop effective potential, with cou-plings determined at tree level (dashed orange line) and one loop at large-N c (solid blue line) with a sigma mass of 500 MeV. Due to the term ∝ −∆ 4 log ∆ 2 m 2 q in Eq. (38), the potential will always be unbounded from below for large values of ∆. However, only in the case where the parameters consistently have been determined at the same accuracy of the effective potential, is there a local minimum such that we actually can study the phase transition at finite T and µ. Using tree-level matching for m σ = 500 MeV, leads to a vacuum effective potential that cannot be used.

Appendix C: Integrals
In order to renormalize the PQM and χM models, we need to evaluate some vacuum integrals in four dimensions. These integrals are divergent in the ultraviolet and we regularize them using dimensional regularization in d = 4 − 2 dimension and the MS scheme. We define The integrals needed are where the functions are F (P 2 ) = 2 − 2q arctan 1 s , (C5) and s = 4m 2 q P 2 − 1. We also need some three-dimensional integrals. In this case the integrals are defined as in Eq. (C1) but now with d = 3 − 2 instead of d = 4 − 2 . We also use q instead of Q as integration variable to distinguish the two cases. The integrals are In this appendix, we show the relation between the fermion propagator and the Polyakov loop in the nonrelativistic limit, i.e. for heavy quark masses. We follow Lowell and Weisberger [57] to construct the nonrelativistic limit of the fermion sector in QCD. We first define the operator U as where the sum over latin indices is only over spatial components. We also define a new fermion field Ψ via It can be easily shown that the operator U is unitary. The quark part of the Lagrangian can then be expanded in powers of m −1 as We next define where q andq † are column N c -plets and the upper and lower two-component spinors of Ψ. If we use the Dirac representation of the gamma matrices in which γ 0 = diag(1, 1, −1, −1), the Lagrangian (D3) can be written as L ψ = q † (−m + i∂ t + gA a 0 T a )q +q(m + i∂ t + gA a 0 T a )q † . (D5) In the Dirac representation, the upper and lower components of the Dirac spinors can be interpreted as the particle and antiparticle, and Eq. (D5) shows that the quark and antiquark degrees of freedom decouple in this limit. To get the Lagrangian into the final form, we usẽ qT aq † =q i T a ijq † j = −q † j T a ijq i = (q † ) TT aqT , where we have definedT a = −(T a ) T . A partial integration yields q∂ tq † −(∂ tq )q † = (q † ) T ∂ tq T , and redefiningq † to be a row object andq a column object, ie. (q † ) T →q † and q T →q, the Lagrangian becomes L ψ = q † (−m + i∂ t + gA a 0 T a )q +q † (−m + i∂ t + gA a 0T a )q = q † Dq +q †Dq , where we have defined the operators Finally, the Hamiltonian density is given by The quark Hamiltonian thus is We now want to evaluate the quantity q a (x, 0)|q a (x, −iβ) , which is the zero-temperature Green's function, [G(x, t; x, 0)] aa = q a (x, 0)| e −iHqt |q a (x, 0) , (D11) analytically continued to imaginary time t = −iτ with τ = β for a quark state evolving under H q . Furthermore, A 0 contained in H q is a classical background field. Since L q is quadratic in the quark fields, the propagator is in practice given by a free quantum field theory. The propagator for a quadratic Lagrangian L = q † Dq is the solution to the equation DG(x, t; x , 0) = iδ(x − x )δ(t). (D12) With D as defined in (D7), we find [i∂ t + gA a 0 (x, t)T a − m] G(x, t; x , 0) = iδ(x − x )δ(t) .
(D13) When the delta functions are zero, this is just the Schrödinger equation, which for a time-dependent Hamiltonian H(t) has the well known solution where T is the time ordering operator. With the delta functions included we see by insertion that a solution is given by where θ(t) is the Heaviside step function. This is the retarded propagator, which we have chosen since we work in the non-relativistic limit.
In analytically continuing this formula to imaginary times, we will have that that G(−iτ, x; x , 0) = 0 for imaginary times τ < 0. This is because we should analytically continue to imaginary time before carrying out the path integral implicit in (D11) to get an analogue of (D12) in imaginary time. We find where we have used that δ(x = 0) = V −1 . Thus the fermion propagator analytically continued to imaginary times is proportional to the Polyakov loop. The vanishing of the latter implies the vanishing of the former and is taken as a definition of confinement.