Dual simulation of a Polyakov loop model at finite baryon density: phase diagram and local observables
Oleg Borisenko, Volodymyr Chelnokov, Emanuele Mendicelli, Alessandro Papa
DDual simulation of a Polyakov loop model atfinite baryon density:phase diagram and local observables
O. Borisenko † , V. Chelnokov , , ∗ , E. Mendicelli ‡ , A. Papa , ¶ Bogolyubov Institute for Theoretical Physics,National Academy of Sciences of Ukraine,03143 Kyiv, Ukraine Istituto Nazionale di Fisica Nucleare, Gruppo collegato di Cosenza,I-87036 Arcavacata di Rende, Cosenza, Italy Institut f¨ur Theoretische Physik, Goethe-Universit¨at Frankfurt,Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany Department of Physics and Astronomy, York University,Toronto, ON, M3J 1P3, Canada Dipartimento di Fisica, Universit`a della Calabria,I-87036 Arcavacata di Rende, Cosenza, Italy
Abstract
Many Polyakov loop models can be written in a dual formulation whichis free of sign problem even when a non-vanishing baryon chemical potentialis introduced in the action. Here, results of numerical simulations of a dualrepresentation of one such effective Polyakov loop model at finite baryondensity are presented. We compute various local observables such as energydensity, baryon density, quark condensate and describe in details the phasediagram of the model. The regions of the first order phase transition andthe crossover, as well as the line of the second order phase transition, areestablished. We also compute several correlation functions of the Polyakovloops. e-mail addresses : † [email protected], ∗ [email protected], ‡ [email protected], ¶ alessandro.papa@fis.unical.it1 a r X i v : . [ h e p - l a t ] N ov Introduction
The properties of strongly interacting matter at finite temperatures and densitiesremain 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 fromsatisfactory, especially at finite baryon chemical potential, due to the famous signproblem. Many approaches to solve this problem, partially or completely, havebeen designed during the last decades and a certain progress has been achievedwithin such methods as Taylor expansion and reweighting at small baryon chemicalpotential, simulations at imaginary potential, complex Langevin simulations andsome others (see, e.g. , the reviews [2–4]).One of the approaches attempting to fully solve the sign problem relies onrewriting the original partition function and important observables in terms ofdifferent, usually integer valued, degrees of freedom such that the resulting Boltz-mann weight is positive definite. Conventionally, all such formulations are referredto as dual formulations, though sometime a different name can be used ( e.g. , fluxline representation) [5]. There are several routes to construct such dual theory fornon-Abelian lattice models with fermions [6–8]. While a dual formulation witha positive Boltzmann weight has not yet been constructed for full QCD, positiveformulations (or formulations where sign problem appears to be very soft) arealready known for few important cases. One such case refers to the strong cou-pling limit of QCD, where the SU ( N ) lattice gauge theory can be mapped ontoa monomer-dimer and closed baryon loop model [9] (for a recent development ofthis direction, see Ref. [10] and references therein). Another important case is rep-resented by many effective Polyakov loop models which can be derived from thefull lattice QCD in certain limits. A dual representation with positive Boltzmannweight is known for some SU ( N ) [11–13] and U ( N ) Polyakov loop models [6] (forrecent advances, see Ref. [14]). Two of these versions have been studied numeri-cally in [12, 13, 15]. The emphasis in these simulations was put on establishing thephase diagram of the model in the presence of the baryon chemical potential andon computing local observables which can be obtained by differentiating the dualpartition function with respect to some of the parameters entering the action ofthe theory.Important class of observables not yet computed in dual formulations are thecorrelations of the Polyakov loops. These correlations can be related to screening(electric and magnetic) masses at finite temperatures. Understanding the proper-ties of such masses would lead to an essential progress in our comprehension of thehigh temperature QCD phase as a whole. While correlations and related masseshave been subject of numerous and intensive calculations at zero chemical potential(see [16] and references therein), it seems to be an extremely difficult problem tocompute these masses for the real baryon chemical potential with available simula-tion methods. So far correlations and screening masses have been computed onlyat imaginary chemical potential in [17]. A closely related and intriguing problem isthe appearance of a hypothetical oscillating phase at finite density [18–20]. Such aphase is ultimately connected to the complex spectrum of the theory and requirescomputations of long-distance correlations with real baryon chemical potential.Here and in a forthcoming paper we study a somewhat different, but equivalentdual form of the effective Polyakov loop model presented in [14]. This form of the2ual representation has been already used by us in [21] for the computation ofcorrelation functions related to the three-quark potential. We believe it is wellsuited to address the problem of screening masses at finite densities, at least inthe framework of the available positive dual formulations. In the present paperwe describe the Polyakov loop model we work with, its dual representation andseveral observables. We compute also some local observables and reveal the phasestructure of the model. We shall also present preliminary results for the Polyakovloop correlations. In a companion paper we will give a detailed study of screen-ing masses, based on the computation of correlation functions and of the secondmoment correlation length at finite density. This will also allow us to draw someconclusions about the existence of an oscillating phase in the model.We work on a 3-dimensional hypercubic lattice Λ = L , with L the linearextension and a unit lattice spacing. The sites of the lattice are denoted by (cid:126)x ≡ x = ( x , x , x ), x i ∈ [0 , L − l = ( (cid:126)x, ν ) is the lattice link in the direction ν ; e ν is the unit vector in the direction ν and N t is the lattice size in the temporaldirection. Periodic boundary conditions are imposed in all directions. Let G bethe SU ( N ) group and U ( x ) an element of G , then dU denotes the (reduced) Haarmeasure on G and Tr U the fundamental character of G .In this paper we shall study an effective 3-dimensional Polyakov loop modelwhich describes a (3 + 1)-dimensional lattice gauge theory with one flavor of stag-gered fermions. The general form of the partition function of the model is givenby Z Λ ( β, m, µ ; N, N f ) ≡ Z = (cid:90) (cid:89) x dU ( x ) (cid:89) x,ν B g ( β ) (cid:89) x B q ( m, µ ) . (1)Here, B g ( β ) is the gauge part of the Boltzmann weight and B q ( m, µ ) is the deter-minant for static quarks. There are many forms of B g ( β ) and B q ( m, µ ) discussedin the literature. In what follows we use the weight that can be obtained on ananisotropic lattice and in the limit of vanishing spatial gauge coupling β s afterexplicit integration over all spatial gauge fields (see, for instance, Refs. [13, 22, 23]and references therein): B g ( β ) = exp (cid:2) β ReTr U ( x )Tr U † ( x + e ν ) (cid:3) . (2)For SU ( N ) the effective coupling constant β is related to the temporal coupling β t by β = 2 D F ( β t ) with D F ( β t ) = (cid:18) C F ( β t ) N C ( β t ) (cid:19) N t , C F ( β t ) = ∞ (cid:88) k = −∞ det I λ i − i + j + k ( β t ) (cid:12)(cid:12)(cid:12) ≤ i,j ≤ N , (3)where I n ( x ) is the modified Bessel function and λ i refers to the fundamental rep-resentation of SU ( N ) and is equal to λ i = δ i . The Boltzmann weight of staticstaggered fermions can be presented as B q ( m, µ ) = A ( m ) det [1 + h + U ( x )] det (cid:2) h − U † ( x ) (cid:3) , (4)where the determinant is taken over group indices and A ( m ) = h − N , h ± = he ± µphT , h = e − N t arcsinh m ≈ e − mphT . (5)3elow we shall use the dimensionless quantities m = m ph /T and µ = µ ph /T . Thenfor SU (3) one has A ( m ) = e m .The resulting Polyakov loop model we work with takes the form Z = (cid:90) (cid:89) x dU ( x ) exp (cid:34) β (cid:88) x,ν ReTr U ( x )Tr U † ( x + e ν ) (cid:35) × (cid:89) x A ( m ) det [1 + h + U ( x )] det (cid:2) h − U † ( x ) (cid:3) . (6)In this model the matrices U ( x ) play the role of Polyakov loops, the only gauge-invariant operators surviving the integration over spatial gauge fields and overquarks. 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 discretetransformations U ( x ) → ZU ( x ), with Z ∈ Z ( N ). This is the global Z ( N ) symme-try. The quark contribution violates this symmetry explicitly. Another importantfeature of the Boltzmann weight is that it becomes complex in the presence of achemical potential, as it follows from (6). Therefore, the model cannot be directlysimulated if µ is non-zero.In the absence of static quarks the Polyakov loop model exhibits a first orderphase transition at the critical point β c ≈ . Z ( N ) symmetry getsspontaneously broken above β c . At the finite density the model defined in Eq. (6)and its several variations have been studied both numerically via simulation ofdual formulations [12, 13, 15] and analytically via mean-field approximation [24,25]. Mean-field and Monte Carlo study are in quantitative agreement for theexpectation values of energy density and Polyakov loop. This allowed to revealthe 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 studyand give further details on the behavior of local observables, including the baryondensity and the quark condensate. Also, first results for Polyakov loop correlationswill be presented.The paper is organized as follows. In Section 2 we formulate our dual rep-resentation of the model valid for all SU ( N ) groups and in all dimensions. Wepresent also results of an analytical study of the model based on strong couplingexpansion 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 alsosimulation results for some local observables as the baryon density and the quarkcondensate. In Section 4 we present preliminary results for the Polyakov loop cor-relation functions. The summary and outline for future work is done in Section 5. In this Section we describe the dual form of the partition function (6). This dualrepresentation will be used in the next Sections for numerical simulations of themodel. All details of the derivation can be found in [14]. In the case of one flavorof staggered fermions the partition function (6) can be presented, after an exact4ntegration over Polyakov loops, as Z = ∞ (cid:88) { r ( l ) } = −∞ ∞ (cid:88) { s ( l ) } =0 (cid:89) l (cid:0) β (cid:1) | r ( l ) | +2 s ( l ) ( s ( l ) + | r ( l ) | )! s ( l )! (cid:89) x A ( m ) R N ( n ( x ) , p ( x )) , (7) n ( x ) = d (cid:88) i =1 (cid:18) s ( l i ) + 12 | r ( l i ) | (cid:19) + 12 d (cid:88) ν =1 ( r ν ( x ) − r ν ( x − e ν )) , (8) p ( x ) = d (cid:88) i =1 (cid:18) s ( l i ) + 12 | r ( l i ) | (cid:19) − d (cid:88) ν =1 ( r ν ( x ) − r ν ( x − e ν )) , (9)where l i , i = 1 , ..., d are 2 d links attached to a site x and R N ( n, p ) = ∞ (cid:88) q = −∞ N (cid:88) k,l =0 (cid:88) σ (cid:96) n + k δ n + k,p + l + qN d ( σ/ k ) d ( σ + q N / l ) h k + h l − . (10)The sum over σ runs over all partitions of n + k , and d ( σ/ m ) is the dimen-sion of a skew representation defined by a corresponding skew Young diagram, σ + q N = ( σ + 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 fac-tors entering the Boltzmann weight of (7) are positive. Hence, this representationis suitable for numerical simulations. The Kronecker delta-function in expres-sion (10) represents the N -ality constraint on the admissible configurations of theinteger-valued variables s ( l ) and r ( l ). This constraint can be exactly resolved onlyin the pure gauge model when h ± = 0. In this case the dual representation (7) hasbeen already tested by us on an example of 2-dimensional SU (3) model, where westudied correlation functions and three-quark potential [21].In the following Sections we study the dual representation (7) via MonteCarlo simulations for the 3-dimensional SU (3) model. In this case the function R N ( n, p ; h ± ) takes the form R ( n, p ) = Q ( n + 1 , p ) (cid:0) h + + h − + h + h − + h h − (cid:1) (11)+ Q ( n, p ) (cid:0) h + h − + h h − (cid:1) + Q ( n, p + 1) (cid:0) h − + h + h h − + h h − (cid:1) + Q ( n + 1 , p + 1) (cid:0) h + h − + h h − (cid:1) + Q ( n + 2 , p ) h + h − + Q ( n, p + 2) h h − . The function Q ( n, p ) is the result of the group integration and is given by [26] Q N ( n, p ) = (cid:88) λ (cid:96) min( n,p ) d ( λ ) d ( λ + | q | N ) , (12)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 quantitiescan be computed with the help of this dual representation. Explicit expressions forthe correlation functions of the Polyakov loops will be given in Section 4. Belowwe list some local observables which we are going to compute here and which canbe obtained by taking suitable derivatives of the partition function (7).5 Magnetization and its conjugate M = (cid:104) Tr U ( x ) (cid:105) , M ∗ = (cid:10) Tr U † ( x ) (cid:11) , (13) • Susceptibility χ = L (cid:0)(cid:10) (Tr U ( x )) (cid:11) − (cid:104) Tr U ( x ) (cid:105) (cid:1) , (14) • Energy density E = 13 L ∂ ln Z∂β = 23 βL (cid:88) l (cid:104) s ( l ) + | r ( l ) |(cid:105) , (15) • Baryon density B = 1 L ∂ ln Z∂µ = 1 L (cid:88) x (cid:104) k ( x ) − l ( x ) (cid:105) , (16) • Quark condensate Q = 1 L ∂ ln Z∂m = N − L (cid:88) x (cid:104) k ( x ) + l ( x ) (cid:105) . (17)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 presentsome results obtained by simple analytical methods. These results can serve asan additional check of numerical data and lead to a better understanding of thewhole phase diagram of the model and of the behavior of the different expectationvalues. The formulation (7) allows a straightforward expansion in powers of β . The freeenergy of the 3-dimensional SU (3) model can be written as F = 3 m + ln R (0 ,
0) + ∞ (cid:88) k =1 β k f k . (18)The first three coefficients f k are given by f = 3 p p f = 34 (cid:0) p + p p − p p + 5 p p + 5 p p + 10 p p p (cid:1) f = 18 (cid:0) p p − p p p − (cid:0) p p p + p p p (cid:1) ++ 20 (cid:0) p p + p p (cid:1) + 150 p (cid:0) p p + p p (cid:1) ++ 84 p p (cid:0) p + p p (cid:1) + 60 (cid:0) p p p + p p p (cid:1) ++ 30 p ( p p + p p ) + 15 p ( p p + p p ) ++ 15 p ( p p + p p ) + p p + 3 p p ) , (19)6here p kl = R ( k, l ) /R (0 ,
0) and the coefficients R can be easily calculated fromEqs. (11) and (12), resulting in R (0 ,
0) = 1 + h + h + h + 2 h cosh(3 µ ) ,R (0 ,
1) = h e − µ (cid:0) h (cid:1) + he µ (cid:0) h + h (cid:1) ,R (1 ,
0) = h e µ (cid:0) h (cid:1) + he − µ (cid:0) h + h (cid:1) ,R (0 ,
2) = he − µ (cid:0) h (cid:1) + h e µ (cid:0) h (cid:1) ,R (2 ,
0) = he µ (cid:0) h (cid:1) + h e − µ (cid:0) h (cid:1) ,R (1 ,
1) = 1 + 2 h + 2 h + h + 2 h cosh(3 µ ) ,R (0 ,
3) = R (3 ,
0) = (cid:0) h (cid:1) + 2 h cosh(3 µ ) ,R (1 ,
2) = 2 h e − µ (cid:0) h (cid:1) + he µ (cid:0) h + 2 h (cid:1) ,R (2 ,
1) = 2 h e µ (cid:0) h (cid:1) + he − µ (cid:0) h + 2 h (cid:1) . (20)All local observables listed above can be obtained from the expansion (18). Thestrong coupling expansion converges, presumably in the region β ≤ β c ( h, µ ), where β c ( h, µ ) is the phase transition or crossover point. One expects that in this regionnumerical data agree reasonably well with strong coupling results. Another obvious approach which can be used to get qualitative description of themodel (6) is mean-field approximation. Within the mean-field method one obtainsan approximate phase diagram of the model. Also, it allows to calculate variouslocal 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 Polyakovloop model in [24, 25]. In this scheme the mean-field approximation reduces to thefollowing replacement: (cid:88) x,ν
ReTr U ( x )Tr U † ( x + e ν ) −→ d (cid:88) x (cid:0) ω Tr U ( x ) + u Tr U † ( x ) (cid:1) , (21)where u = (cid:104) Tr U ( x ) (cid:105) = 1 dβ ∂ ln Z mf ( u, ω ) ∂ω , ω = (cid:104) Tr U † ( x ) (cid:105) = 1 dβ ∂ ln Z mf ( u, ω ) ∂u . (22)The partition function gets the form Z = [ Z mf ( u, ω )] L d , (23) Z mf ( u, ω ) = A ( m ) (cid:90) dU exp (cid:2) dβ (cid:0) ω Tr U + u Tr U † (cid:1)(cid:3) × det [1 + h + U ] det (cid:2) h − U † (cid:3) . (24)Using the integration methods developed in [14, 26] the mean-field partition func-tion is presented as Z mf ( u, ω ) = ∞ (cid:88) r,s =0 ( dβω ) r r ! ( dβu ) s s ! R ( r, s ) (25)7nd explicit form of R ( r, s ) reads R ( r, s ) = (cid:88) i,j =0 C ij Q ( r + i, s + j ) + C Q ( r + 2 , s ) + C Q ( r, s + 2) , (26)where the coefficients C ij are given by C = 2 cosh 3 m + 2 cosh 3 µ , C = 2 cosh m , C = e µ , C = e − µ ,C = 2 e − µ cosh 2 m + 2 e µ cosh m , C = 2 e µ cosh 2 m + 2 e − µ cosh m . (27)The mean-field equations (22), as well as all local observables, can now be com-puted numerically and compared with the strong coupling results and numericaldata. In the pure gauge case, h = 0, one finds a first order phase transition at thecritical value β g ≈ . h = 0 . , µ = 0 .
5. More mean-field results andcomparison with simulations will be given in the next Section.
To explore the phase structure of the model, we simulate numerically the partitionfunction (7). Important ingredient of such simulations is the way we treat thetriality constraint in (10). In the pure gauge theory, h = 0, this constraint issolved exactly in terms of genuine dual variables and the resulting theory can besimulated 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 formulationevery 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 insteadof more complicated worm-like algorithms usually adopted to probe dual modelformulations. The values of the function Q ( n, p ) are computed beforehand andstored in an array, which is then used in simulations. An extra-benefit of theabsence of the triality condition is the possibility to calculate a correlation functionas the expectation value of a product of one-site observables instead of a productover the path connecting the sources (see Section 4).Let us mention that our dual formulation allows to simulate all SU ( N ) modelson 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 presentapproach can be easily extended to any values of N .For simulations above h = 0 .
01, we performed 2 · thermalization updates,and then made measurements every 100 whole lattice updates, collecting a statis-tic of 10 measurements. To estimate statistical error, a jackknife analysis wasperformed at different blocking over bins with size varying from 10 to 10 .When we move to smaller values of h , we see that transition probabilities ofthe one-spin Metropolis update become smaller. To counter this, while keepingthe simulation time reasonable, we used the combination of a multihit update with8 .0 0.1 0.2 0.3 0.4 0.50.51.01.52.02.53.0 Β M 0.0 0.1 0.2 0.3 0.4 0.501234567 Β E0.0 0.1 0.2 0.3 0.4 0.50.70.80.91.01.1 Β B 0.0 0.1 0.2 0.3 0.4 0.50.70.80.91.01.11.2 Β Q Figure 1: Comparison of the observables obtained from mean-field analysis (solidlines), strong coupling expansion up to third order in β (dashed lines) and simu-lation results (dots) at h = 0 . µ = 0 . lattice. Top left: magnetization(blue and red lines denote magnetization and its conjugate, respectively); top right:energy density; bottom left: baryon density; bottom right: quark condensate.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 probabilitythat the update is rejected and generate the number n rej of rejects before the firstacceptance (from the corresponding geometric distribution); in this way we canskip n rej updates and therefore save on the number of updates left to perform onthe current link. Then, if we still need to do some updates, we choose the updatefor the link using the probabilities calculated in presumption that the updateis accepted. Using this technique allows us to vary n hit for different simulationparameters in the range 5 - 2000, while keeping the simulation time constant fora fixed acceptance rate.Below h = 0 .
01 we performed 2 · thermalization updates, with measurementstaken every 50 whole lattice updates, collecting a statistics of 10 measurements.The bins size in the jackknife analysis varied from 10 to 10 .9 .2 Critical behavior A clear indication of the presence of different phases can be seen from the inspectionof the distribution and the scatter plot of the Monte Carlo equilibrium values forthe absolute value of the magnetization near a transition value of β : for somechoice of the parameters h and µ , we observe two separate peaks and two separatespots, respectively, which is suggestive of a first order transition; for other choiceswe see instead just one single peak and a single spot, respectively, as expected fora second order transition or a crossover.Another indication comes from the behavior of the magnetization and its sus-ceptibility versus β near the transition for different values of ( h, µ ). Figs. 2-4 showthat, when h is kept fixed and µ is varied from µ = 0 . to µ = 2 . , the transition soft-ens, the “jump” of the magnetization at the pseudo-critical β becoming less steepand the peak of the susceptibility less pronounced, which suggests that differentregimes are being explored. According to the finite-size scaling analysis discussedbelow, the three regimes represented in Figs. 2-4 correspond to first order, secondorder and crossover, respectively. Β M 0.268 0.269 0.270200400600800 ΒΧ Figure 2: Magnetization (left) and magnetization susceptibility (right) versus β at h = 0 .
01 and µ = 0 on a 16 lattice, in the vicinity of a first order phase transition.The solid lines represent the mean-field estimates.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 valueof the magnetization susceptibility χ . Since we had to explore a three-parameterspace, we could not afford to perform high-statistic simulations if not on latticeswith linear sizes L = 10 and 16. We determined χ for several β values in thetransition region and fitted them to a Lorentzian, thus getting the position of thepeak, which gives the pseudocritical coupling β pc , and its height. Comparing thedependence of the peak height on the lattice size L with the scaling law χ L ( β pc ) = AL γ/ν , (28)we estimated the critical-index ratio γ/ν and collected all our determinations, asmany as 202, in Table 1. 10 .22 0.26 0.3 0.340.00.51.01.52.02.5 Β M 0.2663 0.2664 0.2665 0.2666 0.2667200250300350 ΒΧ Figure 3: The same as Fig. 2 for h = 0 .
01 and µ = 0 . Β M 0.2535 0.2545 0.2555 0.2565141618202224 ΒΧ Figure 4: The same as Fig. 3 for h = 0 .
01 and µ = 2 .
0, where a crossover transitionoccurs. 11 µ β p c γ / ν . . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) - . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) . . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . ( ) . ( ) h µ β p c γ / ν . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) h µ β p c γ / ν . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) - . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) h µ β p c γ / ν . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) . . ( ) . ( ) T a b l e : Su mm a r y o f t h e d e t e r m i n a t i o n s o f β p c a nd γ / ν f o r d i ff e r e n t v a l u e s o f h a nd µ .
12e can see that, within uncertainties, the values of γ/ν are spread in a range be-tween 3., which implies a first order transition, and zero, which holds for crossover,passing through the second order 3-dimensional Ising value, γ/ν = 1 . γ/ν are evidently an artifact of the relatively small latticesizes we could simulate. If we could approach the thermodynamic limit, we wouldsee 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 expectedsince we know that at µ = 0 in the pure gauge limit of QCD or for heavy enoughquark masses there is a whole region of first order deconfinement transitions inthe m u , d - m s plane (the famous Columbia plot), delimited by a line of second ordercritical points in the 3-dimensional Ising class [28]: thereafter, for lower quarkmasses, the crossover region is met. In the simulations of our effective Polyakovloop model at non-zero density toward the thermodynamic limit we should see thecontinuation of the line of second order critical points to non-zero values of thechemical potential. For the lattice volumes considered in our study, we are notable to make a clear-cut assignment of each choice of the parameters h and µ toone of the three transition regions. Using the determination of γ/ν , we tried any-how 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 thethermodynamic limit, but can be helpful in the present context. In Table 2 weintroduced a color code, to help identifying at a glance in which of these regions agiven parameter pair ( h, µ ) falls. γ/ν color phase γ/ν ≥ ≤ γ/ν < < γ/ν < ≤ γ/ν ≤ ≤ γ/ν < ≤ γ/ν < ≤ γ/ν < γ/ν .In Fig. 5(left) each parameter pair ( h, µ ) considered in our simulations is rep-resented by a colored dot in the ( h, µ ) plane, according to the color code defined inTable 2, allowing us to sketch a tentative phase diagram in the right panel of thesame figure. A different visualization of the general phase diagram is presented inFig. 6, where in a 3 d plot γ/ν values are reported in correspondence of each ( h, µ )pair.Fig. 7 shows the values of β pc obtained for each considered pair ( h, µ ). It isinteresting to note that the red points (close to the second order transition) lieapproximately in a one plane, while all points associated with other phases lie ona 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 findthe critical value β g at h = 0, i.e. for the pure gauge theory, and at µ = 0 . ,13 .001 0.01 0.0156 0.03 0.04 μ CROSSOVERFIRST ORDERSECOND ORDER LINE μ h h Figure 5: (Left) Assignment of each parameter pair ( h, µ ) to a transition regionaccording to the color code of Table 2. (Right) Estimated phase diagram.we performed a simple fit of the form β c ( h ) = β g + ah , as suggested in [13]. Weobtained following values β g = 0 . a = − . β g agrees verywell with the value quoted in the literature, β g = 0 .
274 and reasonably well withthe mean-field result β g = 0 . h → µ → ∞ . From the data wehave, one cannot make unambiguous conclusions about its behavior. Nevertheless,data are well fitted by the function µ c = − a ln h + c − bh , with a = 0 . c = − . b = 2367. This shows that the line of second order phase transition might persistin the heavy-dense limit of QCD. In this subsection we study the behavior of the quark condensate and the baryondensity in different phases and compare numerical results with mean-field predic-tions. In the static approximation for the quark determinant one cannot observethe phenomenon of spontaneous breaking of the chiral symmetry. Indeed, even inthe strong coupling region, using Eq. (18), one can easily obtain for the quark con-densate Q = 0 for all µ in the massless limit h = 1. The same result in this limitcan be obtained within mean-field approach and from numerical simulations whichwe performed for various values of β and µ . Nevertheless, we think it might be in-structive to see the behavior of the condensate in the three regimes correspondingto first and second order transitions and to crossover. A more traditional observ-able to study in the effective Polyakov loop models is the baryon density. In thedual formulation the baryon density B and the quark condensate Q are given byEqs. (16) and (17), respectively. We have computed the right-hand sides of theseequations both numerically and using the mean-field approach, for the same valuesof h and µ .The behavior of baryon density and quark condensate as functions of β dependsstrongly on the phase of the system. Left panels of Figs. 8 and 9 show the typical14 / ν = ( ) μ γ / ν 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 . γ/ν value for a second order phase transition in the 3-dimensional Isingclass [27] .behavior of Q and B for h = 0 .
01 and various values of the chemical potentialcorresponding to different phases of the model. These phases are characterized asbefore by the values of γ/ν and can be approximately read off from Fig. 6 or, moreprecisely from the Table 1. One observes a rapid change of both Q and B at thefirst order phase transition. This rapid change becomes smoother and smootherwhen parameters are gradually changed toward the second order line and then tothe crossover regime. The right panels of the same Figures compare Monte Carlodata with the mean-field predictions in the three regions. One can conclude thatmean-field reproduces numerical simulations with good accuracy.The quark condensate as a function of β is shown in Fig. 10(left) for vanishingvalue of the chemical potential and several values of h . The right panel of thesame figure demonstrates the approach to the saturation of the baryon density atzero quark mass, h = 1.As follows from Fig. 1, the baryon density B is a decreasing function of β atsufficiently 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 methodsof calculations used in this paper. Therefore, there should exists some value of thequark mass where the density changes its qualitative behavior. We have not triedto estimate this value. To see the impact of a non-zero chemical potential on the correlation functionbehavior, we calculated the two-point correlation functions for several values ofparameters. We considered six kinds of the correlation functions:15 β pc h Figure 7: Values of β pc for each considered choice of the pair ( h, µ ). The colorcode is as in Table 2. The red points have been made bigger than the others tohighlight them (see comment in the text).Γ nn ( r ) = (cid:104) Tr U (0) Tr U ( r ) (cid:105) , Γ rr ( r ) = (cid:104) Re Tr U (0) Re Tr U ( r ) (cid:105) , Γ na ( r ) = (cid:10) Tr U (0) Tr U † ( r ) (cid:11) , Γ ri ( r ) = (cid:104) Re Tr U (0) Im Tr U ( r ) (cid:105) , Γ aa ( r ) = (cid:10) Tr U † (0) Tr U † ( r ) (cid:11) , Γ ii ( r ) = (cid:104) Im Tr U (0) Im Tr U ( r ) (cid:105) . (29)In the dual formulation the correlation functions can be written asΓ nn ( r ) = (cid:28) R ( n (0) + 1 , p (0)) R ( n (0) , p (0)) R ( n ( r ) + 1 , p ( r )) R ( n ( r ) , p ( r )) (cid:29) Γ na ( r ) = (cid:28) R ( n (0) + 1 , p (0)) R ( n (0) , p (0)) R ( n ( r ) , p ( r ) + 1) R ( n ( r ) , p ( r )) (cid:29) Γ aa ( r ) = (cid:28) R ( n (0) , p (0) + 1) R ( n (0) , p (0)) R ( n ( r ) , p ( r ) + 1) R ( n ( r ) , p ( r )) (cid:29) . (30)These formulas work for r >
0. For r = 0 both shifts to the n , p variables happenat one point, so only one ratio remains. The correlations Γ rr , Γ ri and Γ ii can beobtained as linear combinations of Γ nn , Γ na and Γ aa .The expressions (30) become unusable when h = 0, and can have a bad con-vergence properties for very small h , or very large µ values. We have checked bycomparing the numerical results with the strong coupling expressions for small β values, that the results can be relied on for h > .
005 and µ < h , the average traces can become non-zero, mak-ing the correlation function constant even in the disordered phase. Due to that,16 µ β pc γ/ν β pc and γ/ν for h and µ belonging to the region “very close tosecond order”. These values correspond to the bigger red points in Fig. 7.we introduce the subtracted correlation functions, subtracting the correspondingaverage trace inside the correlations:Γ nn, sub ( r ) = Γ nn ( r ) − (cid:104) Tr U (cid:105) , Γ na, sub ( r ) = Γ na ( r ) − (cid:104) Tr U (cid:105) (cid:10) Tr U † (cid:11) , Γ aa, sub ( r ) = Γ aa ( r ) − (cid:10) Tr U † (cid:11) . (31)For these subtracted correlations we expect an exponential decay,Γ( r ) = A exp( − mr ) r , (32)at least in the disordered phase.Samples of correlation function behavior in different regions of the phase dia-gram are shown in Figs. 11, 12 and 13. One can see that, indeed, both in disorderedand ordered phases, the correlations decay exponentially. While the mass gap, cor-responding 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 the17 ●●●●●●●●● ●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● β μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) h = Q ● ● ●●●●●●●●●●●●●●●●●●●●●●●●●●● ●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● β MF μ = MF μ = MF μ = μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) h = Q Figure 8: Behavior of the quark condensate versus β for fixed h = 0 .
01 and fordifferent values of µ on a 16 lattice. Left panel: results of simulations in thevicinity of a phase transition. Right panel: comparison of mean-field analysis andnumerical results in three different phase regimes. The color code of Table 2 isused here and below to differentiate the phases of the system.mass gap for the Γ ii correlation remains more or less constant, and in particulardoes 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 differentmass gaps m M and m E , m M < m E (see [17]). While at non-zero µ these twosectors should mix, so all the correlators we study should decay with m M , it ispossible that in our case the mixing is small causing the real large distance massgap for the color-electric sector to be visible only on distances larger than the onesfor 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 growswith µ . This means increasing of the screening effects: at finite density non-zero µ pushes system deeper in the deconfined phase. This is in qualitative accordancewith the results of [17] obtained with imaginary µ .We will address these questions in more detail in a future work, which is underpreparation. For now, we can note that at least in the disordered phase alsothe second moment correlation length for the imaginary-imaginary correlations issubstantially different from the one for the real-real correlations. Revealing the phase diagram of QCD at finite temperature and non-zero baryonchemical potential, as well as the nature of strong interacting matter at high tem-peratures, remains one of the important challenges of high energy physics. In spiteof enormous efforts, many aspects of these problems are still far from unambigu-ous resolution. The several methods developed to solve the sign problem in finitedensity QCD have their own advantages and drawbacks. Dual formulations, asthe one used here, solve the sign problem completely, but are restricted so far tostrong coupled regions in the case of non-Abelian gauge theories. Even in this18 ●●●●●●●●● ●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● β μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) h = B ●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ● β MF μ = MF μ = MF μ = μ = γ / ν =( ± ) μ = γ / ν =( ± ) μ = γ / ν =( ± ) h = B Figure 9: Behavior of the baryon density versus β for fixed h = 0 .
01 and fordifferent values of µ on a 16 lattice. Left panel: results of simulations in thevicinity of a phase transition. Right panel: comparison of mean-field analysis andnumerical results in three different phase regimes.region one can get valuable information about the phase structure of non-Abelianmodels and behavior of various observables. The study presented here is in thesame spirit of Refs. [12, 13, 15], though we have used a different dual formulationof the Polyakov loop model, built in Ref. [14]. To corroborate our findings we havealso compared in many cases simulation results with the strong coupling expansionof the dual model and with the mean-field analysis. Let us briefly recapitulate ourmain results. • The phase diagram of the model was studied in great details. We haveclassified three regions in the parameter space of the model ( β, h, µ ) accordingto the type of the critical behavior: first or second order phase transition,or crossover. The values of the ratio of critical indices γ/ν is different indifferent regimes. • As main observables we computed expectation values of the Polyakov loopand its conjugate, the baryon density and the quark condensate. • Our dual formulation allows us to compute correlation functions of thePolyakov loops. We have presented some preliminary results for such corre-lations at non-zero chemical potential. • It is interesting to note that the mean-field results agree very well withnumerical simulations both at zero and non-zero µ . Also, the mean-fieldresults are in good qualitative agreement with a similar analysis in Ref. [24].The overall qualitative picture of the phase diagram and the behavior of allobservables fully agree with the picture described in Refs. [13, 15].All observables considered in this work have shown sensitivity to the chemicalpotential. The general trend is that when µ is increased, they exhibit a lesssteep variation across transition when the coupling β (which corresponds to thetemperature in the underlying QCD theory) is increased. Qualitatively, one can19 ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●● ● ● ● ● ● ● ● ● ● ●●●●●●●●●●●●●●●● β h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) h = γ / ν =( ± ) μ = Q Μ B Figure 10: Left panel: behavior of the quark condensate versus β for fixed µ = 0and for different values of h on a 16 lattice. Right panel: baryon density for h = 1and β = 0 . lattice.say that increasing µ plays effectively the same role as a reduction of the quarkmass.The most important direction for the future work is a detailed study of thedifferent correlations of the Polyakov loops and the extraction of screening (electricand magnetic) masses at finite chemical potential. Also, an investigation of theoscillating phase and the related complex masses can be accomplished within ourdual formulation. All these problems will be addressed in a companion paperwhich is currently under preparation. Acknowledgements
Numerical simulations have been performed on the ReCaS Data Center ofINFN-Cosenza. O. Borisenko acknowledges support from the National Academy ofSciences of Ukraine in frames of priority project ”Fundamental properties of matterin the relativistic collisions of nuclei and in the early Universe” (No. 0120U100935).E.M. was supported in part by a York University Graduate Fellowship Doctoral -International. A.P. acknowledges support from the INFN/NPQCD project.
References [1] O. Philipsen, PoS
LATTICE2019 , 273 (2019) [arXiv:1912.04827 [hep-lat]].[2] C. Schmidt, PoS
LAT2006 , 021 (2006) [arXiv:hep-lat/0610116 [hep-lat]].[3] P. de Forcrand, PoS
LAT2009 , 010 (2009) [arXiv:1005.0539 [hep-lat]].[4] E. Seiler, EPJ Web Conf. , 01019 (2018) [arXiv:1708.08254 [hep-lat]].[5] C. Gattringer, K. Langfeld, Int. J. Mod. Phys. A , no.22, 1643007 (2016)[arXiv:1603.09517 [hep-lat]]. 206] O. Borisenko, V. Chelnokov, S. Voloshyn, EPJ Web Conf. , 11021 (2018)[arXiv:1712.03064 [hep-lat]].[7] C. Marchis, C. Gattringer, Phys. Rev. D , no.3, 034508 (2018)[arXiv:1712.07546 [hep-lat]].[8] G. Gagliardi, W. Unger, Phys. Rev. D , no.3, 034509 (2020)[arXiv:1911.08389 [hep-lat]].[9] F. Karsch, K.H. M¨utter, Nucl.Phys. B (1989) 541.[10] M. Klegrewe, W. Unger, Phys. Rev. D , no.3, 034505 (2020)[arXiv:2005.10813 [hep-lat]].[11] C. Gattringer, Nucl. Phys. B , 242-252 (2011) [arXiv:1104.2503 [hep-lat]].[12] Y. D. Mercado, C. Gattringer, Nucl. Phys. B , 737-750 (2012)[arXiv:1204.6074 [hep-lat]].[13] M. Fromm, J. Langelage, S. Lottini, O. Philipsen, JHEP , 042 (2012)[arXiv:1111.4953 [hep-lat]].[14] O. Borisenko, V. Chelnokov, S. Voloshyn, Phys. Rev. D , no.1, 014502(2020) [arXiv:2005.11073 [hep-lat]].[15] Y. Delgado, C. Gattringer, Acta Phys. Polon. Supp. , 1033-1038 (2012)[arXiv:1208.1169 [hep-lat]].[16] A. Bazavov and J. H. Weber, “Color Screening in Quantum Chromodynam-ics”, [arXiv:2010.01873 [hep-lat]], to appear in Progress in Particle and Nu-clear Physics.[17] M. Andreoli, C. Bonati, M. D’Elia, M. Mesiti, F. Negro, A. Rucci, F. Sanfil-ippo, Phys. Rev. D , no.5, 054515 (2018) [arXiv:1712.09996 [hep-lat]].[18] P. N. Meisinger, M. C. Ogilvie, T. D. Wiser, Int. J. Theor. Phys. (2011)1042 [arXiv:1009.0745 [hep-th]].[19] H. Nishimura, M. C. Ogilvie, K. Pangeni, Phys. Rev. D , no.9, 094501(2016) [arXiv:1512.09131 [hep-lat]].[20] O. Akerlund, P. de Forcrand, T. Rindlisbacher, JHEP , 055 (2016)[arXiv:1602.02925 [hep-lat]].[21] O. Borisenko, V. Chelnokov, E. Mendicelli, A. Papa, Nucl. Phys. B ,214-238 (2019) [arXiv:1812.05384 [hep-lat]].[22] M. Bill`o, M. Caselle, A. D’Adda, S. Panzeri, Int. J. Mod. Phys. A , 1783-1846 (1997) [arXiv:hep-th/9610144 [hep-th]].[23] J. Langelage, M. Neuman and O. Philipsen, JHEP , 131 (2014)[arXiv:1403.4162 [hep-lat]].[24] J. Greensite, K. Splittorff, Phys. Rev. D , 074501 (2012) [arXiv:1206.1159[hep-lat]]. 2125] J. Greensite, Phys. Rev. D , no.11, 114507 (2014) [arXiv:1406.4558 [hep-lat]].[26] O. Borisenko, S. Voloshin, V. Chelnokov, Rept. Math. Phys. , 129-145(2020) [arXiv:1812.06069 [hep-lat]].[27] A. Pelissetto, E. Vicari, Phys. Rept. , 549-727 (2002) [arXiv:cond-mat/0012164 [cond-mat]].[28] F. Karsch and S. Stickan, Phys. Lett. B , 319-325 (2000) [arXiv:hep-lat/0007019 [hep-lat]]. 22 r = 0.267 , h = 0.008 , = 0.9635 nnnaaarrriii r = 0.2681 , h = 0.008 , = 0.9635 r = 0.269 , h = 0.008 , = 0.9635 Figure 11: Behavior of the correlation functions on a 20 lattice for different valuesof parameter β ( h = 0 . µ = 0 . r = 0.265 , h = 0.01 , = 0.9635 nnnaaarrriii r = 0.26649 , h = 0.01 , = 0.9635 r = 0.268 , h = 0.01 , = 0.9635 Figure 12: The same as Figure 11 for h = 0 . µ = 0 . r = 0.25 , h = 0.01 , = 2 nnnaaarrriii r = 0.25543 , h = 0.01 , = 2 r = 0.25 , h = 0.01 , = 2 Figure 13: The same as Figure 11 for h = 0 . µ = 2 ..