Dual simulation of a Polyakov loop model at finite baryon density: correlations and screening masses

Computations of screening masses in finite-temperature QCD at finite density are plagued by the sign problem and have been performed so far with an imaginary chemical potential. Here, we use a dual formulation of a Polyakov-loop model which allows the determination of screening masses at real baryon chemical potential. This is a second paper in a series devoted to a detailed study of dual Polyakov-loop models at finite density. While the first paper was mainly devoted to establishing the phase diagram of the model, here we compute correlation functions of the Polyakov loops and the second-moment correlation length at non-zero chemical potential. This enables us to evaluate numerically the screening masses from correlations of the real and imaginary parts of the Polyakov loops. We also compute these masses in the mean-field approximation and compare with numerical results. In addition, we provide a quantitative improvement of the general phase diagram presented in the first paper.


Introduction
Understanding the properties of strongly interacting matter at finite temperatures and densities is exceedingly important in many areas of research, including cosmology, astrophysics of compact objects and phenomenology of heavy-ion collisions.Numerical investigations on a Euclidean space-time lattice, based on Monte-Carlo methods, are hindered by the notorious sign problem: in presence of a non-zero baryon chemical potential the Boltzmann weight turns complex and cannot be used to sample thermal ensemble configurations.Many approaches have been designed in the last decades to circumvent or mitigate this problem, such as Taylor expansion around zero baryon chemical potential, reweighting to small baryon chemical potential, simulations at imaginary potential, complex Langevin simulations and others (see, e.g., [1][2][3][4]).
Recent times have witnessed the development of radically new approaches that aim to fully solve the sign problem by expressing the original partition function and observables in terms of different, usually integer-valued, degrees of freedom, so to make the resulting Boltzmann weight positive definite.These formulations are conventionally called "dual", even though sometimes different names are used, such as "flux line representation" [5].
While several strategies have been developed to get dual formulations for non-Abelian lattice models with fermions [6][7][8], a dual formulation with a positive Boltzmann weight is not yet available for full QCD.In some limiting regimes, however, it has been possible to construct formulations either positive or such that the sign problem is less severe.This is the case of the strong-coupling limit of QCD, where the SU (N ) lattice gauge theory (LGT) can be mapped onto a monomer-dimer and closed baryon loop model [9] (see also Ref. [10] and references therein).Other cases include many effective Polyakov-loop models which can be derived from the full lattice QCD in some specific limits.Examples of such dual models with positive Boltzmann weight are known in the wider context of the SU (N ) [11][12][13][14] and U (N ) groups [7].
In this paper we focus on a 3-dimensional effective Polyakov-loop model, describing the (3 + 1)-dimensional SU (N ) LGT with one flavor of staggered fermions at finite baryon density.It is defined on a 3-dimensional hypercubic lattice Λ = L 3 , with L the linear extension and a unit lattice spacing; ⃗ x ≡ x = (x 1 , x 2 , x 3 ), x i ∈ [0, L − 1] denote the sites of the lattice, l = (⃗ x, ν) is the lattice link in the direction ν; e ν is a unit vector in the direction ν and N t is the lattice size in the temporal direction of the underlying (3+1)-dimensional LGT; periodic boundary conditions are imposed in all directions.The general form of the partition function of the model reads [15,16] Z Λ (β, m, μ; 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.TrU denotes the fundamental character of SU (N ).Integration in (1) is performed with respect to the Haar measure on SU (N ).The effective coupling constant β is a complicated function of the original gauge coupling (its precise form is not important here).The constants A( m) and h ± are given by The pure-gauge part of the effective action is invariant only under global discrete transformations U (x) → ZU (x), 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 follows from (1).Therefore, the model cannot be studied by direct Monte-Carlo methods if µ is non-zero.The model defined in Eq. ( 1) has been investigated in the large-N limit of U (N ) at non-zero µ in [17]; it was found that the dependence on the chemical potential drops out from the free energy.Moreover, a third-order phase transition has been reported in this limit.The large-N limit of the SU (N ) case was studied in [18,19].In Refs.[20,21] a mean-field approximation was used to study the SU (3) model.
Dual formulations with positive Boltzmann weight of the model defined in Eq. ( 1) have been constructed in Refs.[7,11], thus making Monte-Carlo simulations feasible [12,13,22].Mean-field and Monte-Carlo studies turned out to be in quantitative agreement for the energy density and expectation values of the Polyakov loops.In some regions of the parameters β, h and µ the phase diagram of the model has been revealed, but long-distance correlations and, hence the screening masses, have not been computed so far (see, however, [23]).
In a previous paper [24], we have studied an alternative dual form of the partition function (1), originally derived in [14] (its explicit expression is given below).By a finite-size scaling of the magnetization susceptibility, we have studied in great details the phase diagram of the model, identifying three regions in the parameter space (β, h, µ) according to the type of the critical behavior: first or second order phase transition, or crossover.We have computed expectation values of Polyakov loop, baryon density and quark condensate, finding good agreement with meanfield predictions both at zero and non-zero µ.In general, we have found that the overall qualitative picture of the phase diagram and the behavior of all observables fully agree with Refs.[13,22].We have also observed that, in general, all the observables considered in that work are sensitive to the chemical potential, exhibiting a less steep variation across transition when the coupling β (which corresponds to the temperature in the underlying QCD theory) is increased.Qualitatively, the effect of increasing µ has shown to play the same role of the reduction of the quark mass.
In the present paper we exploit the same dual form of the effective Polyakovloop model adopted in [24] to study several correlations of the Polyakov loops and to extract screening (electric and magnetic) masses at finite chemical potential.We consider also the determination of the second-moment correlation length.To the best of our knowledge, this kind of analysis has never been done so far in the framework of dual formulations of SU (N ) effective Polyakov-loop models.It is important not only to size the effect of a non-zero baryon chemical potential on screening effects in the deconfined phase, but also for the investigation of the elusive oscillating phase at finite density [25][26][27], which is ultimately connected to the complex spectrum of the theory.
For a general review of the screening masses in lattice QCD at finite temperature we refer to [28].Simulations at imaginary chemical potential with two flavors of Wilson fermions and with (2 + 1) flavor of staggered fermions have been performed in Ref. [29] and Ref. [30], correspondingly.
We describe now the dual form of the partition function (1).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 (1) 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 ( 3) is valid for all SU (N ) groups and in any dimension.Clearly, all factors entering the Boltzmann weight of (3) are positive.Hence, this representation is suitable for numerical simulations.The Kronecker delta-function in expression (6) 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 (3) has been already tested by us on the example of a 2-dimensional SU (3) model, where we studied correlation functions and three-quark potential [23].
In the following Sections we study the dual representation (3) 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 [31] Q 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 3.
The paper is organized as follows.In the next Section we present a mean-field study of the correlations for SU (3) Polyakov-loop model and calculate the screening masses analytically.In Section 3 we give the definitions of several Polyakovloop correlators both in the standard and in the dual formulations; we define also in this context the second-moment correlation length.Then, we present our numerical Monte-Carlo results for Polyakov-loop correlation functions, screening masses and second-moment correlation length.In Section 4 we summarize our main results and outline the future work.In Appendix we present improved phase diagram of the Polyakov-loop model.

Analytic study of the correlations
First of all, we would like to derive some analytic predictions for the behavior of the correlations and screening masses.Here we calculate the partition and correlation functions of SU (3) model following the approach of Refs.[18,19] developed for the large-N limit.All results presented in this and next Sections are given in terms of dimensionless quantities m = m/T and µ = μ/T .
Consider the following change of variables: The partition function can then be rewritten as The full action in new variables takes the form where S g is the gauge part of the action J(x) is the Jacobian of the transformation [32] and Q(x) represents the quark contribution The integration in ( 10) is performed over the domain where J(x) > 0. Firstly, we look for constant, translation invariant solutions of the saddle-point equations 18dβρ + 1 2J Secondly, we expand the full action around the saddle points.If ρ 0 , ω 0 are the saddle points, we make the substitution ρ(x) → ρ 0 + 1 3 ρ(x), ω(x) → ω 0 + 1 3 ω(x) and expand in powers of fluctuations.The action becomes Here, S cl is the classical action and coefficients b i are functions of ρ 0 , ω 0 .These coefficients are very lengthy and cumbersome, so we choose not to give here their precise expressions.Important is that b 1 and b 2 are always real and non-vanishing, while b 3 is purely imaginary and non-zero only when the chemical potential is non-vanishing.Finally, we calculate the correlation functions by integrating over quadratic fluctuations.For the arbitrary normalized correlations of the fundamental characters of the Polyakov loops one obtains We have used the following notation for the Green function: Constants C 1,2 define the screening masses and are given below.This general result allows one to compute any invariant or non-invariant observable by choosing the appropriate values of the sources η(x) and η(x).E.g., the magnetization M = ⟨ρ(x)e iω(x) ⟩ and its complex conjugate become where G 0 (m) is the zero-distance Green function.Since we have two components of the mean Polyakov loop ( ⟨Tr U ⟩ and Tr U † , or equivalently ⟨Re Tr U ⟩ and Im Tr U † ), we have four components of the Polyakov loop correlation function, which can be gathered in the correlation matrix [30] Γ The off-diagonal terms vanish if µ = 0, and the coefficients in the exponential decay of the connected parts in diagonal terms define the magnetic m M = m 1 and electric m E = m 2 screening masses.For µ > 0 the electric and magnetic sector mix, so each correlation matrix element should be a sum of two terms -one decaying with m M , and second with m E .Eigenvalues of the correlation matrix ( 25) turn out to be Notice that, if m 1 ≤ m 2 , one obtains in the limit of large separation R: ) Then, for example, the connected part of the Polyakov loop correlation becomes The masses m 1 and m 2 play the role of the magnetic and electric screening masses and are given by These expressions for screening masses are the main result of the mean-field analysis.Let us make a few comments on these expressions.
• Using explicit expressions for the coefficients b i , one proves that m 2 ≥ 2m 1 , i.e., the electric mass is at least two times larger than the magnetic mass.
Monte-Carlo simulations support this conclusion.
• If µ is non-vanishing then b 3 ̸ = 0 and masses may, in principle, become complex and m 1 = m * 2 .In this case the correlation would decay exponentially and exhibit oscillatory behavior as follows from Eq. (31).Such behavior is typical for a liquid phase which exists in Z(3) spin model in the presence of external complex magnetic field [27] and in the large-N limit of the SU (N ) Polyakov-loop model [19].However, for the SU (3) model studied here we were not able to find a set of parameters for which the product C 1 C 2 is negative and masses are complex.Either the oscillating region is very narrow or it is not realized in this model.
• The behavior versus β of the two screening masses predicted by analytic expressions (32), (33) is shown in Fig. 1.This can be compared with results of Monte-Carlo simulations presented below.One finds a reasonable qualitative, and sometime even quantitative, agreement.
One can see that both masses grow when going away from the transition point, with the higher mass having very rapid growth after entering the ordered phase.Both masses have their minimum at the transition point, with lower mass getting closer to zero at the transition.3 Monte-Carlo results for correlations and screening masses In this Section we present results of simulations and describe the behavior of the screening masses extracted from the exponential decay of the Polyakov-loop correlations.These masses will be compared with the masses computed from the second-moment correlation length.Details of the lattice setup and update are described in Ref. [24].In this study we used lattices with sizes L = 20, 24 and sometimes L = 32.The phase diagram described in Ref. [24] was obtained from simulations mainly on the lattice L = 20.The larger lattices we used here gave us the possibility to quantitatively improve the phase diagram of the model.Such improved diagram is briefly described in the Appendix.

Two-point correlation functions
To see the impact of non-zero chemical potential on the correlation function behavior, we calculated the values of the two-site correlation functions for several values of parameters.We considered six kinds of correlation functions: The correlation matrix (25) takes the form In the dual formulation the correlation functions can be written as These formulas work for R > 0. For R = 0 both shifts to the n, p variables happen at the same point, so only one ratio remains.The correlations Γ rr , Γ ri and Γ ii can be obtained as linear combinations of Γ nn , Γ nc and Γ cc .
The expressions (36) become unusable when h = 0, and can have 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 large-distance correlation function constant even in the disordered phase.Due to that, we introduce the connected correlation functions: One expects an exponential decay for these connected correlations at least in the disordered phase.Samples of correlation function behavior at different points of the phase diagram are shown in Fig. 2. One can see that, indeed, both in disordered and ordered phases, the correlations decay exponentially.An interesting property is that, while the mass gap, corresponding to the slope of the plots, remains the same for Γ nn , Γ nc , Γ cc , Γ rr and Γ ri correlation functions, it is much larger for the Γ ii correlation function.Also the mass gap for the Γ ii correlation remains more or less constant, and in particular does not vanish in the vicinity of the phase transition.At µ = 0 the off-diagonal terms Γ ri (R) are zero, and one can define magnetic and electric correlation masses as the exponential decay rates for, correspondingly, real-real and imaginary-imaginary connected correlations, At non-zero µ the correlation matrix has non-zero off-diagonal elements, and each element has contributions from two masses: While for small µ values the coefficients for the contribution of the smaller mass to Γ ii , and of the larger mass to Γ rr are small, they make the extraction of the masses (especially the larger one) difficult.Diagonalizing the correlation matrix Γ(R), one can get contributions that depend purely on m 1 or m 2 like in µ = 0 case, which we fit to the exponential decay behavior to obtain the two masses.The masses extracted from fitting the diagonalized correlation matrix are shown in Fig. 3.One can see that the smaller mass behavior defines the transition: it rapidly drops to a small, but non-zero value, at the first-order transition point; it drops to zero (or, more precisely to a small value that goes to zero with L → ∞) at the second-order transition point; it has a smooth minimum in the case of crossover.The higher mass remains constant and of the order of one in the disordered phase, and starts to grow rapidly around the transition, thus meaning that the corresponding correlation drops very fast in the ordered phase and vanishes at distances of a couple of lattice spacings, making therefore very difficult to estimate the mass in that phase.Such behavior agrees well with the mean-field prediction m 2 ≥ 2m 1 .

Second moment correlation length
Since in many cases the distances on which the correlation length can be determined reliably does not exceed five lattice spacings, the mass gap extracted from fitting the exponential decay of the correlation functions has a large uncertainty.
An alternative to this approach is to extract from data second-moment correlation length of the (diagonalized) correlation functions Γ 1 , Γ 2 : where the summation is taken over all vectors R = (R 1 , R 2 , R 3 ) on the lattice.
Assuming Γ j (R) = A e −mR R , and taking the thermodynamic limit L → ∞, we can replace summation by the integration, and trigonometric functions by their Taylor expansions.Then the integration gives ξ 2,j = 1 m + O(1/L).Existence of larger masses m k in the correlation function spectrum introduces a finite correction that becomes larger when the ratio m k /m becomes closer to 1. Thus, in general ξ 2,j ̸ = 1 m even in thermodynamic and continuum limit.We introduce a "second-moment mass gap" m (ξ 2 ) j = 1/ξ 2,j as an approximation to the actual mass gap m.To make this approximation closer to the actual mass gap, and to attempt extracting m 2 , we have to perform the diagonalization of the correlation matrix and take the second moment of the diagonalized correlation functions, thus removing corrections from m 2 to m 1 .Since χ j and F j are linear in Γ j , we apply diagonalization to the χ and F matrices to obtain χ 1,2 and F 1,2 -the moments of the correlation functions corresponding to the lower and higher masses.
The "second-moment mass gaps" extracted in this way are compared with the mass gaps extracted from the correlation function fits in Fig. 3.One can see that the lower masses extracted from the correlation function fits and from the second-moment correlation length are in good agreement, especially close to the transition.The higher masses get very large errors in the ordered phase for µ = 0, and the diagonalization needed for non-zero µ results in large errors for the higher masses at any β.The behavior of the higher mass in the disordered phase at µ = 0 is the same as for the mass extracted from the correlation function fit, up to a shift, that can be explained by a corrections from higher masses.: Dependence on β of the lower and higher screening masses extracted from the diagonalized correlation function fits (blue and violet markers), and from the diagonalized second-moment correlation lengths (yellow and green markers) around the transition point for h = 0.01, µ = 0, L = 16 (left, first order phase transition), h = 0.01, µ = 0.9635, L = 24 (center, second order phase transition), and h = 0.01, µ = 2.0, L = 16 (right, crossover).The data points where the error estimate is above 1/3 of the mass value are removed for plot clarity.

Summary
In this paper we studied the dual formulation of SU (3) lattice gauge theory with one flavor of static staggered fermions.In this approximation the original theory can be first mapped onto an effective Polyakov-loop model.The latter can be further mapped onto a dual formulation with a positive weight.The main idea of this study is to reveal the behavior of the screening masses and the second moment correlation length at finite baryon density.We used a simple mean-field approximation and Monte-Carlo simulations to accomplish such study.Let us briefly summarize our main results.
• Using the larger lattices and better statistics we have quantitatively improved the phase diagram presented in Ref. [24].No qualitative changes have been found.
• To get an idea of how the screening masses behave, we used the mean-field approximation to calculate the correlation matrix of the Polyakov loops.Two screening masses have been extracted analytically from the correlations of the real and imaginary parts of the Polyakov loops.In all cases their behavior agrees with expectations.Especially good agreement with Monte-Carlo data is found in the vicinity of the second-order phase transition.
• We have computed numerically all possible correlations of the Polyakov loops.The screening masses and the second-moment correlation length have been obtained from the fitting of the connected correlations to the exponential decay.In all three cases (first-, second-order phase transitions and crossover) one observes a reasonable agreement with the expectations and with the mean-field results for the lower magnetic mass.The second, higher mass grows very rapidly after the transitions, preventing us from obtaining good fits.The general trend is that when µ is increased, the screening masses exhibit a less steep variation across transition when the coupling β (which corresponds to the temperature in the underlying QCD theory) is increased.
• Both the mean-field and numerical results do not show the existence of the complex mass spectrum in the theory which was found in two-dimensional QCD with the static quarks [25,26], in Z(3) spin model with the complex magnetic field [27] and in the large-N limit of the :Polyakov-loop model similar to the one studied here [19].
The two main approximations used to obtain the dual formulation with a positive Boltzmann weight are the strong coupling approximation for the Wilson action and calculation of the quark determinant in the static approximation.In Ref. [14] we have presented a dual formulation valid at all values of the temporal gauge coupling constant.This formulation preserves the positivity of the Boltzmann weight.Therefore, one of the possible direction for the future work is an extension of the present simulations to the dual representation derived in [14].Another direction is to include the leading corrections to the static quark determinant.It is not yet clear at the moment if such corrections preserve the positivity of the dual weight.This problem is currently under investigation.
We can see that, within uncertainties, the values of γ/ν are spread in a range between 3, which implies a first-order transition, and 0, which holds for crossover, passing through the second order 3-dimensional Ising value, γ/ν = 1.9638 (8) [33].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 [34]: thereafter, for lower quark masses, the crossover region is met.In the simulations of our effective Polyakovloop 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 to make this assignment, extending and modifying the three possible options (first order, second order and crossover) as seen in Table 1.This makes no sense in the thermodynamic limit, but can be helpful in the present context.In Table 1 we introduced a color code, to identify which of these regions a given parameter pair (h, µ) falls in.
In Fig. 4(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 1, allowing us to sketch a tentative phase diagram in the right panel of the same figure.
We have not performed simulations in the absence of an 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 the following values β g = 0.2741(2), a = −0.50(2).The value of β g agrees very well with the value quoted in the literature, β g = 0.2741 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. 4 in the heavy-dense limit, h → 0, µ → ∞.From the data we have, one cannot make unambiguous conclusions about its behavior.Nevertheless, This shows that the line of second order phase transition might persist in the heavy-dense limit of QCD.

Figure 2 :
Figure 2: Behavior of the correlation functions on lattice with L = 20 for different values of parameters β, h and µ.Solid lines represent the best fit to the function (38) with periodic contribution (R → L−R).Plots on the left are in the disordered phase, in the middle -near the phase transition (crossover) point, on the rightin the ordered phase.

Figure 3
Figure3: Dependence on β of the lower and higher screening masses extracted from the diagonalized correlation function fits (blue and violet markers), and from the diagonalized second-moment correlation lengths (yellow and green markers) around the transition point for h = 0.01, µ = 0, L = 16 (left, first order phase transition), h = 0.01, µ = 0.9635, L = 24 (center, second order phase transition), and h = 0.01, µ = 2.0, L = 16 (right, crossover).The data points where the error estimate is above 1/3 of the mass value are removed for plot clarity.

Figure 4 :
Figure 4: (Left) Assignment of each parameter pair (h, µ) to a transition region according to the color code of Table 1.(Right) Estimated phase diagram.

Table 2 :
Summary of the determinations of β pc and γ/ν for different values of h and µ.