Dual simulation of a Polyakov loop model at finite baryon density: phase diagram and local observables

Many Polyakov loop models can be written in a dual formulation which is free of sign problem even when a non-vanishing baryon chemical potential is introduced in the action. Here, results of numerical simulations of a dual representation of one such effective Polyakov loop model at finite baryon density are presented. We compute various local observables such as energy density, baryon density, quark condensate and describe in details the phase diagram of the model. The regions of the first order phase transition and the crossover, as well as the line of the second order phase transition, are established. We also compute several correlation functions of the Polyakov loops.


Introduction
The properties of strongly interacting matter at finite temperatures and densities remain in the focus of intensive theoretical and numerical studies (see Ref. [1] for a recent review). The full understanding of these properties is still far from satisfactory, especially at finite baryon chemical potential, due to the famous sign problem. Many approaches to solve this problem, partially or completely, have been designed during the last decades and a certain progress has been achieved within such methods as Taylor expansion and reweighting at small baryon chemical potential, simulations at imaginary potential, complex Langevin simulations and some others (see, e.g., the reviews [2][3][4]).
One of the approaches attempting to fully solve the sign problem relies on rewriting the original partition function and important observables in terms of different, usually integer valued, degrees of freedom such that the resulting Boltzmann weight is positive definite. Conventionally, all such formulations are referred to as dual formulations, though sometime a different name can be used (e.g., flux line representation) [5]. There are several routes to construct such dual theory for non-Abelian lattice models with fermions [6][7][8]. While a dual formulation with a positive Boltzmann weight has not yet been constructed for full QCD, positive formulations (or formulations where sign problem appears to be very soft) are already known for few important cases. One such case refers to the strong coupling limit of QCD, where the SU (N ) lattice gauge theory can be mapped onto a monomer-dimer and closed baryon loop model [9] (for a recent development of this direction, see Ref. [10] and references therein). Another important case is represented by many effective Polyakov loop models which can be derived from the full lattice QCD in certain limits. A dual representation with positive Boltzmann weight is known for some SU (N ) [11][12][13] and U (N ) Polyakov loop models [6] (for recent advances, see Ref. [14]). Two of these versions have been studied numerically in [12,13,15]. The emphasis in these simulations was put on establishing the phase diagram of the model in the presence of the baryon chemical potential and on computing local observables which can be obtained by differentiating the dual partition function with respect to some of the parameters entering the action of the theory. Important class of observables not yet computed in dual formulations are the correlations of the Polyakov loops. These correlations can be related to screening (electric and magnetic) masses at finite temperatures. Understanding the properties of such masses would lead to an essential progress in our comprehension of the high temperature QCD phase as a whole. While correlations and related masses have been subject of numerous and intensive calculations at zero chemical potential (see [16] and references therein), it seems to be an extremely difficult problem to compute these masses for the real baryon chemical potential with available simulation methods. So far correlations and screening masses have been computed only at imaginary chemical potential in [17]. A closely related and intriguing problem is the appearance of a hypothetical oscillating phase at finite density [18][19][20]. Such a phase is ultimately connected to the complex spectrum of the theory and requires computations of long-distance correlations with real baryon chemical potential.
Here and in a forthcoming paper we study a somewhat different, but equivalent dual form of the effective Polyakov loop model presented in [14]. This form of the dual representation has been already used by us in [21] for the computation of correlation functions related to the three-quark potential. We believe it is well suited to address the problem of screening masses at finite densities, at least in the framework of the available positive dual formulations. In the present paper we describe the Polyakov loop model we work with, its dual representation and several observables. We compute also some local observables and reveal the phase structure of the model. We shall also present preliminary results for the Polyakov loop correlations. In a companion paper we will give a detailed study of screening masses, based on the computation of correlation functions and of the second moment correlation length at finite density. This will also allow us to draw some conclusions about the existence of an oscillating phase in the model.
We work on a 3-dimensional hypercubic lattice Λ = L 3 , with L the linear extension and a unit lattice spacing. The sites of the lattice are denoted by is the lattice link in the direction ν; e ν is the unit vector in the direction ν and N t is the lattice size in the temporal direction. Periodic boundary conditions are imposed in all directions. Let G be the SU (N ) group and U (x) an element of G, then dU denotes the (reduced) Haar measure on G and TrU the fundamental character of G.
In this paper we shall study an effective 3-dimensional Polyakov loop model which describes a (3 + 1)-dimensional lattice gauge theory with one flavor of staggered fermions. The general form of the partition function of the model is given by Here, B g (β) is the gauge part of the Boltzmann weight and B q (m, µ) is the determinant for static quarks. There are many forms of B g (β) and B q (m, µ) discussed in the literature. In what follows we use the weight that can be obtained on an anisotropic lattice and in the limit of vanishing spatial gauge coupling β s after explicit integration over all spatial gauge fields (see, for instance, Refs. [13,22,23] and references therein): For SU (N ) the effective coupling constant β is related to the temporal coupling where I n (x) is the modified Bessel function and λ i refers to the fundamental representation of SU (N ) and is equal to λ i = δ 1i . The Boltzmann weight of static staggered fermions can be presented as where the determinant is taken over group indices and Below we shall use the dimensionless quantities m = m ph /T and µ = µ ph /T . Then for SU (3) one has A(m) = e 3m .
The resulting Polyakov loop model we work with takes the form In this model the matrices U (x) play the role of Polyakov loops, the only gaugeinvariant operators surviving the integration over spatial gauge fields and over quarks. The integration in (6) is performed with respect to the Haar measure on G. The pure gauge part of the SU (N ) model is invariant under global discrete transformations U (x) → ZU (x), with Z ∈ Z(N ). This is the global Z(N ) symmetry. The quark contribution violates this symmetry explicitly. Another important feature of the Boltzmann weight is that it becomes complex in the presence of a chemical potential, as it follows from (6). Therefore, the model cannot be directly simulated if µ is non-zero.
In the absence of static quarks the Polyakov loop model exhibits a first order phase transition at the critical point β c ≈ 0.274. The global Z(N ) symmetry gets spontaneously broken above β c . At the finite density the model defined in Eq. (6) and its several variations have been studied both numerically via simulation of dual formulations [12,13,15] and analytically via mean-field approximation [24,25]. Mean-field and Monte Carlo study are in quantitative agreement for the expectation values of energy density and Polyakov loop. This allowed to reveal the phase diagram of the model, at least in some regions of the parameters β, h and µ. In this paper we confirm the qualitative picture found in previous study and give further details on the behavior of local observables, including the baryon density and the quark condensate. Also, first results for Polyakov loop correlations will be presented.
The paper is organized as follows. In Section 2 we formulate our dual representation of the model valid for all SU (N ) groups and in all dimensions. We present also results of an analytical study of the model based on strong coupling expansion and mean-field approximation. The phase diagram of the 3-dimensional SU (3) model is studied numerically in details in Section 3, where we discuss also simulation results for some local observables as the baryon density and the quark condensate. In Section 4 we present preliminary results for the Polyakov loop correlation functions. The summary and outline for future work is done in Section 5.

Dual formulation of the Polyakov loop model
In this Section we describe the dual form of the partition function (6). This dual representation will be used in the next Sections for numerical simulations of the model. All details of the derivation can be found in [14]. In the case of one flavor of staggered fermions the partition function (6) can be presented, after an exact integration over Polyakov loops, as where l i , i = 1, ..., 2d are 2d links attached to a site x and The sum over σ runs over all partitions of n + k, and d (σ/1 m ) is the dimension of a skew representation defined by a corresponding skew Young diagram, σ + q N = (σ 1 + q, . . . , σ N + q) (for more details we refer the reader to Ref. [14]). Equation (7) is valid for all SU (N ) groups and in any dimension. Clearly, all factors entering the Boltzmann weight of (7) are positive. Hence, this representation is suitable for numerical simulations. The Kronecker delta-function in expression (10) represents the N -ality constraint on the admissible configurations of the integer-valued variables s(l) and r(l). This constraint can be exactly resolved only in the pure gauge model when h ± = 0. In this case the dual representation (7) has been already tested by us on an example of 2-dimensional SU (3) model, where we studied correlation functions and three-quark potential [21].
In the following Sections we study the dual representation (7) via Monte Carlo simulations for the 3-dimensional SU (3) model. In this case the function R N (n, p; h ± ) takes the form The function Q 3 (n, p) is the result of the group integration and is given by [26] Q N (n, p) = λ min(n,p) where d(λ) is the dimension of the permutation group S r in the representation λ, q = (p − n)/N (when q is not an integer Q N (n, p) = 0). Important is the fact that both local observables and long-distance quantities can be computed with the help of this dual representation. Explicit expressions for the correlation functions of the Polyakov loops will be given in Section 4. Below we list some local observables which we are going to compute here and which can be obtained by taking suitable derivatives of the partition function (7).
• Magnetization and its conjugate • Susceptibility • Energy density • Baryon density • Quark condensate In the last two equations k(x) and l(x) are the summation variables from (10). Before discussing results of numerical simulations we would like to present some results obtained by simple analytical methods. These results can serve as an additional check of numerical data and lead to a better understanding of the whole phase diagram of the model and of the behavior of the different expectation values.

Mean-field solution
Another obvious approach which can be used to get qualitative description of the model (6) is mean-field approximation. Within the mean-field method one obtains an approximate phase diagram of the model. Also, it allows to calculate various local observables, as the free energy density, the baryon density as some others.
We use here one of the simplest mean-field schemes, applied to a similar Polyakov loop model in [24,25]. In this scheme the mean-field approximation reduces to the following replacement: where The partition function gets the form Using the integration methods developed in [14,26] the mean-field partition function is presented as and explicit form of R 3 (r, s) reads where the coefficients C ij are given by The mean-field equations (22), as well as all local observables, can now be computed numerically and compared with the strong coupling results and numerical data. In the pure gauge case, h = 0, one finds a first order phase transition at the critical value β g ≈ 0.2615. This value, as well as the magnetization, matches the corresponding results obtained earlier by the mean-field method [24,25]. Fig. 1 compares our numerical data with the strong coupling expansion and mean-field results for some typical values h = 0.6, µ = 0.5. More mean-field results and comparison with simulations will be given in the next Section.
3 Phase diagram of the model

Lattice setup
To explore the phase structure of the model, we simulate numerically the partition function (7). Important ingredient of such simulations is the way we treat the triality constraint in (10). In the pure gauge theory, h = 0, this constraint is solved exactly in terms of genuine dual variables and the resulting theory can be simulated with the usual Metropolis update [21]. Instead, at non-zero values of h, the function R N (n, p) in (7) is explicitly expanded in series (11). In this formulation every configuration of link variables s(l) and r(l) has non-zero Boltzmann weight, thus allowing us to use again the simple Metropolis update algorithm instead of more complicated worm-like algorithms usually adopted to probe dual model formulations. The values of the function Q 3 (n, p) are computed beforehand and stored in an array, which is then used in simulations. An extra-benefit of the absence of the triality condition is the possibility to calculate a correlation function as the expectation value of a product of one-site observables instead of a product over the path connecting the sources (see Section 4). Let us mention that our dual formulation allows to simulate all SU (N ) models on equal footing: the only difference between different N is encoded in the function Q N (n, p) which, as said, can be computed prior to simulations. Thus, the present approach can be easily extended to any values of N .
For simulations above h = 0.01, we performed 2 · 10 5 thermalization updates, and then made measurements every 100 whole lattice updates, collecting a statistic of 10 5 measurements. To estimate statistical error, a jackknife analysis was performed at different blocking over bins with size varying from 10 to 10 4 .
When we move to smaller values of h, we see that transition probabilities of the one-spin Metropolis update become smaller. To counter this, while keeping the simulation time reasonable, we used the combination of a multihit update with an update skipping procedure: we go through each link variable and attempt n hit Metropolis updates, but before each attempt we calculate the a priori probability that the update is rejected and generate the number n rej of rejects before the first acceptance (from the corresponding geometric distribution); in this way we can skip n rej updates and therefore save on the number of updates left to perform on the current link. Then, if we still need to do some updates, we choose the update for the link using the probabilities calculated in presumption that the update is accepted. Using this technique allows us to vary n hit for different simulation parameters in the range 5 -2000, while keeping the simulation time constant for a fixed acceptance rate.
Below h = 0.01 we performed 2·10 5 thermalization updates, with measurements taken every 50 whole lattice updates, collecting a statistics of 10 6 measurements. The bins size in the jackknife analysis varied from 10 to 10 5 .

Critical behavior
A clear indication of the presence of different phases can be seen from the inspection of the distribution and the scatter plot of the Monte Carlo equilibrium values for the absolute value of the magnetization near a transition value of β: for some choice of the parameters h and µ, we observe two separate peaks and two separate spots, respectively, which is suggestive of a first order transition; for other choices we see instead just one single peak and a single spot, respectively, as expected for a second order transition or a crossover.
Another indication comes from the behavior of the magnetization and its susceptibility versus β near the transition for different values of (h, µ). Figs. 2-4 show that, when h is kept fixed and µ is varied from µ = 0. to µ = 2., the transition softens, the "jump" of the magnetization at the pseudo-critical β becoming less steep and the peak of the susceptibility less pronounced, which suggests that different regimes are being explored. According to the finite-size scaling analysis discussed below, the three regimes represented in Figs. 2-4 correspond to first order, second order and crossover, respectively. To characterize in a quantitative way the phase structure of the model in the (h, µ)-plane, we performed a standard finite-size scaling analysis on the peak value of the magnetization susceptibility χ. Since we had to explore a three-parameter space, we could not afford to perform high-statistic simulations if not on lattices with linear sizes L = 10 and 16. We determined χ for several β values in the transition region and fitted them to a Lorentzian, thus getting the position of the peak, which gives the pseudocritical coupling β pc , and its height. Comparing the dependence of the peak height on the lattice size L with the scaling law we estimated the critical-index ratio γ/ν and collected all our determinations, as many as 202, in Table 1.  2.83 (6) 3.2 0.267843 (5) 2.50 (11) 3.3 0.267156 (7) 2.14 (5) 3    We can see that, within uncertainties, the values of γ/ν are spread in a range between 3., which implies a first order transition, and zero, which holds for crossover, passing through the second order 3-dimensional Ising value, γ/ν = 1.9638(8) [27]. These sparse values of γ/ν are evidently an artifact of the relatively small lattice sizes we could simulate. If we could approach the thermodynamic limit, we would see that values of γ/ν concentrate around the values of 3. (first order), 1.9638 (8) (second order in the 3-dimensional Ising class) and 0. (crossover). This is expected since we know that at µ = 0 in the pure gauge limit of QCD or for heavy enough quark masses there is a whole region of first order deconfinement transitions in the m u,d -m s plane (the famous Columbia plot), delimited by a line of second order critical points in the 3-dimensional Ising class [28]: thereafter, for lower quark masses, the crossover region is met. In the simulations of our effective Polyakov loop model at non-zero density toward the thermodynamic limit we should see the continuation of the line of second order critical points to non-zero values of the chemical potential. For the lattice volumes considered in our study, we are not able to make a clear-cut assignment of each choice of the parameters h and µ to one of the three transition regions. Using the determination of γ/ν, we tried anyhow to make this assignment, extending and modifying the three possible options (first order, second order and crossover) as in Table 2. This makes no sense in the thermodynamic limit, but can be helpful in the present context. In Table 2 we introduced a color code, to help identifying at a glance in which of these regions a given parameter pair (h, µ) falls.
In Fig. 5(left) each parameter pair (h, µ) considered in our simulations is represented by a colored dot in the (h, µ) plane, according to the color code defined in Table 2, allowing us to sketch a tentative phase diagram in the right panel of the same figure. A different visualization of the general phase diagram is presented in Fig. 6, where in a 3d plot γ/ν values are reported in correspondence of each (h, µ) pair. Fig. 7 shows the values of β pc obtained for each considered pair (h, µ). It is interesting to note that the red points (close to the second order transition) lie approximately in a one plane, while all points associated with other phases lie on a curved surface. This is also seen from Table 3, where we collected the values of β pc close to the second order phase transition.
We have not performed simulations in the absence of external field. To find the critical value β g at h = 0, i.e. for the pure gauge theory, and at µ = 0., we performed a simple fit of the form β c (h) = β g + ah, as suggested in [13]. We obtained following values β g = 0.274991, a = −0.5568. The value of β g agrees very well with the value quoted in the literature, β g = 0.274 and reasonably well with the mean-field result β g = 0.2615. Another important question concerns the shape of the critical line shown in the right panel of Fig. 5 in the heavy-dense limit, h → 0, µ → ∞. From the data we have, one cannot make unambiguous conclusions about its behavior. Nevertheless, data are well fitted by the function µ c = −a ln h+c−bh 2 , with a = 0.932, c = −3.1, b = 2367. This shows that the line of second order phase transition might persist in the heavy-dense limit of QCD.

Quark condensate and baryon density
In this subsection we study the behavior of the quark condensate and the baryon density in different phases and compare numerical results with mean-field predictions. In the static approximation for the quark determinant one cannot observe the phenomenon of spontaneous breaking of the chiral symmetry. Indeed, even in the strong coupling region, using Eq. (18), one can easily obtain for the quark condensate Q = 0 for all µ in the massless limit h = 1. The same result in this limit can be obtained within mean-field approach and from numerical simulations which we performed for various values of β and µ. Nevertheless, we think it might be instructive to see the behavior of the condensate in the three regimes corresponding to first and second order transitions and to crossover. A more traditional observable to study in the effective Polyakov loop models is the baryon density. In the dual formulation the baryon density B and the quark condensate Q are given by Eqs. (16) and (17), respectively. We have computed the right-hand sides of these equations both numerically and using the mean-field approach, for the same values of h and µ.
The behavior of baryon density and quark condensate as functions of β depends strongly on the phase of the system. Left panels of Figs. 8 and 9 show the typical γ/ν=1.9638 (8) μ γ/ν h Figure 6: Values of the critical index ratio γ/ν for each considered choice of the pair (h, µ). The color code is as in Table 2. The plane at γ/ν = 1.9638(8) corresponds to the γ/ν value for a second order phase transition in the 3-dimensional Ising class [27] .
behavior of Q and B for h = 0.01 and various values of the chemical potential corresponding to different phases of the model. These phases are characterized as before by the values of γ/ν and can be approximately read off from Fig. 6 or, more precisely from the Table 1. One observes a rapid change of both Q and B at the first order phase transition. This rapid change becomes smoother and smoother when parameters are gradually changed toward the second order line and then to the crossover regime. The right panels of the same Figures compare Monte Carlo data with the mean-field predictions in the three regions. One can conclude that mean-field reproduces numerical simulations with good accuracy.
The quark condensate as a function of β is shown in Fig. 10(left) for vanishing value of the chemical potential and several values of h. The right panel of the same figure demonstrates the approach to the saturation of the baryon density at zero quark mass, h = 1.
As follows from Fig. 1, the baryon density B is a decreasing function of β at sufficiently large h = 0.6, while Fig. 9 shows that for small h values (large mass) B is an increasing function of β. This conclusion is supported by all three methods of calculations used in this paper. Therefore, there should exists some value of the quark mass where the density changes its qualitative behavior. We have not tried to estimate this value.

Large-distance behavior of the correlations
To see the impact of a non-zero chemical potential on the correlation function behavior, we calculated the two-point correlation functions for several values of parameters. We considered six kinds of the correlation functions:  Table 2. The red points have been made bigger than the others to highlight them (see comment in the text).
In the dual formulation the correlation functions can be written as Γ nn (r) = R 3 (n(0) + 1, p(0)) R 3 (n(0), p(0)) R 3 (n(r) + 1, p(r)) R 3 (n(r), p(r)) Γ na (r) = R 3 (n(0) + 1, p(0)) R 3 (n(0), p(0)) R 3 (n(r), p(r) + 1) R 3 (n(r), p(r)) These formulas work for r > 0. For r = 0 both shifts to the n, p variables happen at one point, so only one ratio remains. The correlations Γ rr , Γ ri and Γ ii can be obtained as linear combinations of Γ nn , Γ na and Γ aa . The expressions (30) become unusable when h = 0, and can have a bad convergence properties for very small h, or very large µ values. We have checked by comparing the numerical results with the strong coupling expressions for small β values, that the results can be relied on for h > 0.005 and µ < 3.
Since we work at non-zero h, the average traces can become non-zero, making the correlation function constant even in the disordered phase. Due to that,  Table 3: Values of β pc and γ/ν for h and µ belonging to the region "very close to second order". These values correspond to the bigger red points in Fig. 7.
we introduce the subtracted correlation functions, subtracting the corresponding average trace inside the correlations: For these subtracted correlations we expect an exponential decay, at least in the disordered phase. Samples of correlation function behavior in different regions of the phase diagram are shown in Figs. 11, 12 and 13. One can see that, indeed, both in disordered and ordered phases, the correlations decay exponentially. While the mass gap, corresponding to the slope of the plots, remains the same for Γ nn , Γ na , Γ aa , Γ rr and Γ ri correlation functions, it is much larger for the Γ ii correlation function. Also the  Table 2 is used here and below to differentiate the phases of the system. mass gap for the Γ ii correlation remains more or less constant, and in particular does not vanish in the vicinity of the phase transition. The difference in mass gaps can be explained by noting that at µ = 0 Γ rr and Γ ii correspond to the color-magnetic and color-electric sectors having different mass gaps m M and m E , m M < m E (see [17]). While at non-zero µ these two sectors should mix, so all the correlators we study should decay with m M , it is possible that in our case the mixing is small causing the real large distance mass gap for the color-electric sector to be visible only on distances larger than the ones for which we have reliable results.
In the ordered phase we observe increasing of the correlation function slope (at the same β values) with the increase of µ, implying that the mass gap grows with µ. This means increasing of the screening effects: at finite density non-zero µ pushes system deeper in the deconfined phase. This is in qualitative accordance with the results of [17] obtained with imaginary µ.
We will address these questions in more detail in a future work, which is under preparation. For now, we can note that at least in the disordered phase also the second moment correlation length for the imaginary-imaginary correlations is substantially different from the one for the real-real correlations.

Summary
Revealing the phase diagram of QCD at finite temperature and non-zero baryon chemical potential, as well as the nature of strong interacting matter at high temperatures, remains one of the important challenges of high energy physics. In spite of enormous efforts, many aspects of these problems are still far from unambiguous resolution. The several methods developed to solve the sign problem in finite density QCD have their own advantages and drawbacks. Dual formulations, as the one used here, solve the sign problem completely, but are restricted so far to strong coupled regions in the case of non-Abelian gauge theories. Even in this region one can get valuable information about the phase structure of non-Abelian models and behavior of various observables. The study presented here is in the same spirit of Refs. [12,13,15], though we have used a different dual formulation of the Polyakov loop model, built in Ref. [14]. To corroborate our findings we have also compared in many cases simulation results with the strong coupling expansion of the dual model and with the mean-field analysis. Let us briefly recapitulate our main results.
• The phase diagram of the model was studied in great details. We have classified three regions in the parameter space of the model (β, h, µ) according to the type of the critical behavior: first or second order phase transition, or crossover. The values of the ratio of critical indices γ/ν is different in different regimes.
• As main observables we computed expectation values of the Polyakov loop and its conjugate, the baryon density and the quark condensate.
• Our dual formulation allows us to compute correlation functions of the Polyakov loops. We have presented some preliminary results for such correlations at non-zero chemical potential.
• It is interesting to note that the mean-field results agree very well with numerical simulations both at zero and non-zero µ. Also, the mean-field results are in good qualitative agreement with a similar analysis in Ref. [24].
The overall qualitative picture of the phase diagram and the behavior of all observables fully agree with the picture described in Refs. [13,15].
All observables considered in this work have shown sensitivity to the chemical potential. The general trend is that when µ is increased, they exhibit a less steep variation across transition when the coupling β (which corresponds to the temperature in the underlying QCD theory) is increased. Qualitatively, one can  say that increasing µ plays effectively the same role as a reduction of the quark mass.
The most important direction for the future work is a detailed study of the different correlations of the Polyakov loops and the extraction of screening (electric and magnetic) masses at finite chemical potential. Also, an investigation of the oscillating phase and the related complex masses can be accomplished within our dual formulation. All these problems will be addressed in a companion paper which is currently under preparation.    Figure 11 for h = 0.01, µ = 2.0 (crossover region).