# Gluon Propagators in 2+1 Lattice QCD with Nonzero Isospin Chemical Potential

V. G. Bornyakov, A. A. Nikolaev, R. N. Rogalyov, A. S. Terentev

GGluon propagators in lattice QCD with nonzero isospin chemical potential

V. G. Bornyakov

NRC “Kurchatov Institute”– IHEP Protvino, 142281, Russian Federation,NRC “Kurchatov Institute”— ITEP, Moscow 117218, Russian Federation,School of Biomedicine, Far Eastern Federal University,Sukhanova 8, Vladivostok, 690950, Russian Federation

A. A. Nikolaev

Department of Physics, College of Science, Swansea University, Swansea SA2 8PP, United Kingdom

R. N. Rogalyov

NRC ”Kurchatov Institute” - IHEP, Protvino 142281, Russian Federation

A. S. Terentev

National University of Science and Technology ”MISIS”, Leninsky Avenue 4, 119991 Moscow, Russia

The static longitudinal and transverse gluon propagators in the Landau gauge are studied in2 + 1 lattice QCD with nonzero isospin chemical potential µ I . Parameterization of the momentumdependence of the propagators is provided for all values of the chemical potential under study.We ﬁnd that the longitudinal propagator is infrared suppressed at nonzero µ I with suppressionincreasing with increasing µ I . It is found, respectively, that the electric screening mass is increasingwith increasing µ I . Additionally, we analyze the diﬀerence between two propagators as a functionof the momentum and thus compare interactions in chromoelectric and chromomagnetic sectors. PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.AwKeywords: gauge ﬁeld theory, gluon propagator

I. INTRODUCTION

To understand the physics of neutron star matter and results of heavy-ion collision experiments it is importantto study the quark matter with isospin (isotopic) asymmetry. QCD phase diagram at nonzero values of the isospinchemical potential µ I has been studied recently quite intensively using lattice QCD [1–11], eﬀective models as NJLmodel [12–20], quark meson model [21–23], chiral perturbation theory [24–26], perturbative QCD [27, 28], randommatrix model [29], hadron resonance gas model [30] and other approaches, see e.g. recent review [31].It is also believed that the study of QCD at nonzero µ I might help, in one way or another, to solve more diﬃcultproblem of QCD at nonzero baryon chemical potential. First, QCD with nonzero µ I is an ideal test system for somelattice approaches to QCD with nonzero baryon chemical potential such as Taylor expansion, analytical continuationand others. Second, lattice results obtained at nonzero µ I by means of ﬁrst principle calculations can be comparedwith respective results of the eﬀective ﬁeld theories of QCD and phenomenological models. Such comparison allowsto test the eﬀectiveness of QCD models and thus to estimate their applicability to the study of QCD at nonzero µ B .In this paper we make the ﬁrst study of the inﬂuence of µ I on the Landau gauge gluon propagators using latticeQCD approach. We use the same lattice action as in Ref. [10] and in fact the same set of the lattice conﬁgurations. Ourgoal is to study how the gluon propagators change with an increase of µ I . In particular, we ﬁnd strong suppression ofthe (color)-electric (longitudinal) gluon propagator with increasing µ I which is reﬂected in increasing of the respectivescreening mass. We parameterize the propagators in the wide range of momentum values using the ﬁt function whichis the tree level prediction of the Reﬁned Gribov-Zwanziger approach. We also compute an observable introducedrecently in [32], a diﬀerence between the (color-)electric and (color-)magnetic propagators, and obtain its dependenceon the momentum and isospin chemical potential.The gluon propagators are among important quantities to study, e.g. they play crucial role in the Dyson-Schwingerequations approach to QCD. Landau gauge gluon propagators in the non-Abelian gauge theories at zero and nonzerotemperature were extensively studied in the infrared range of momenta by various methods. We shall note latticegauge theory, Dyson-Schwinger equations, Gribov-Zwanziger approach. At the same time the studies at nonzero quarkchemical potential are restricted to a few papers only. Here we close this gap for the case of nonzero isospin chemicalpotential. Let us note that for the lattice QCD at nonzero baryon chemical potential the results can be obtained onlyfor small values of µ B because of the sign problem [33].It is known that in QCD with nonzero µ I there is a phase transition from the normal to the superﬂuid phase at µ I = µ π /

2. Recently results were obtained [10] indicating that at large µ I the theory undergoes smooth transition to a r X i v : . [ h e p - l a t ] F e b the deconﬁnement phase. It is interesting to check if the gluon propagators change at respective values of µ I .The paper is organized as follows. In Section II we specify details of the lattice setup to be used: lattice action,deﬁnition of the propagators and details of the simulation. In Section III we consider the propagators at low momentato compute respective screening masses and determine their dependence on µ I . In Section IV the ﬁt over a widerange of momenta is performed using Gribov-Stingl ﬁt function motivated by the Reﬁned Gribov-Zwanziger eﬀectiveaction. We demonstrate that for the longitudinal propagator this ﬁt works well over the whole momenta domainunder study for all values of µ I considered, whereas for the transverse propagator it works only zero and small µ I .In Section V the dependence of the diﬀerence between the electric and magnetic propagators on the momentumand isospin chemical potential is discussed. We show that this diﬀerence decreases exponentially fast with momentasimilarly to the case of SU (2) QCD [32]. The last section is devoted to the discussion of the results and to conclusionsto be drawn. II. SIMULATION DETAILS

We carry out our study of the 2 + 1 lattice QCD with light quark mass m u = m d ≡ m ud and strange quark mass m s using 28 lattices for a set of the chemical potentials in the range aµ I ∈ (0 , . Z = (cid:90) DU e − S G ( U ) det [ M † ud M ud + λ ] / det [ M s ] / (1)where λ is coupling of the pionic source term, light quarks lattice operator M ud and strange quark lattice operator M s are M ud = /D ( µ I ) + m ud , M s = /D (0) + m s , (2) /D ( µ I ) is the staggered lattice Dirac operator with quark chemical potential µ I = µ u = − µ d . The unphysicalpionic source term which breaks U (1) symmetry explicitly is necessary to observe the spontaneous breaking of thissymmetry in a ﬁnite volume in the limit λ → µ I [1]. The lattice conﬁgurations were generated [10]at β = 4 . , am ud = 0 . , am s = 0 . , λ = 0 . am ud . The lattice spacing of the ensemble in physical units wasdetermined via the Sommer scale value r = 0 .

468 fm [35] to be a = 0 . m π = 380 MeV [10].To reach large quark densities without lattice artifacts one needs suﬃciently small lattice spacing to satisfy condition aµ I (cid:28)

1. At the same time to study the gluon propagators in the infrared region it is necessary to employ largephysical volume. As a result of a compromise between these two requirements our lattice size is rather moderate : L = 28 a = 1 .

92 fm. This implies a potential problem of large ﬁnite volume eﬀects at small momenta.In our study of the gluon propagators we employ the standard deﬁnition of the lattice gauge vector potential A x,µ [36]: A x,µ = 12 iag (cid:16) U xµ − U † xµ (cid:17)(cid:12)(cid:12)(cid:12)(cid:12) traceless . (3)The lattice Landau gauge ﬁxing condition is( ∇ B A ) x ≡ a (cid:88) µ =1 ( A x,µ − A x − a ˆ µ,µ ) = 0 , (4)which is equivalent to ﬁnding an extremum of the gauge functional F U ( ω ) = 14 V (cid:88) x,µ

13 Re Tr U ωxµ , (5)with respect to gauge transformations ω x . To ﬁx the Landau gauge we use the simulated annealing (SA) algorithmwith ﬁnalizing overrelaxation [37, 38].The gluon propagator D abµν ( p ) is deﬁned as follows: D abµν ( p ) = 1 V a (cid:104) (cid:101) A aµ ( q ) (cid:101) A bν ( − q ) (cid:105) , (6)where (cid:101) A bµ ( q ) = a (cid:88) x A bx,µ exp (cid:16) iq ( x + ˆ µa (cid:17) , (7) q i ∈ ( − N s / , N s / q ∈ ( − N t / , N t /

2] and the physical momenta p µ are deﬁned by the relations ap i = 2 sin ( πq i /N s ), ap = 2 sin ( πq /N t ).At nonzero µ I the O (4) symmetry is broken and there are two tensor structures for the gluon propagator [39] : D abµν ( p ) = δ ab (cid:0) P Tµν ( p ) D T ( p ) + P Lµν ( p ) D L ( p ) (cid:1) . (8)We consider the soft modes p = 0 and use the notation D L,T ( p ) = D L,T (0 , | (cid:126)p | ). III. SCREENING MASSES

Various deﬁnitions of the screening masses were discussed in great detail in our previous studies [32, 40]. Here weuse the deﬁnition in terms of the Taylor expansion of D − L,T ( p ) in p at low momenta: D − L,T ( p ) = Z − ( ˜ m E,M + p + c · ( p ) + ... ) . (9) This method was used in [41] in the studies of lattice QCD at ﬁnite temperature and was applied to QC D in [40].One can note that the Yukawa-type ﬁt-function D − L,T ( p ) = Z − ( ˜ m E,M + p ) (10)successfully employed in the previous studies of lattice gluodynamics at zero and ﬁnite temperatures in Refs.[42–44]does not work in the case under consideration because the minimal nonzero momentum on our lattices is rather large: p min ≈

640 MeV. Thus we have no data points in deep infrared domain, where higher-order terms in formula (9) canbe neglected. As explained below, we have to take into account up to the O (cid:16) ( p ) (cid:17) terms of the Taylor expansion ofthe inverse propagators.It should be noted that the screening masses ˜ m E,M are related to the correlation lengths:˜ m E,M = ξ − E,M , (11)where ξ E,M are usually deﬁned in terms of propagators as follows [45]: ξ = 12 (cid:82) V dx d(cid:126)x ˜ D ( x , (cid:126)x ) | (cid:126)x | (cid:82) V dx d(cid:126)xD ( x , (cid:126)x ) = − D (0 ,(cid:126) (cid:88) i =1 (cid:18) ddp i (cid:19) (cid:12)(cid:12)(cid:12) (cid:126)p =0 D (0 , (cid:126)p ) . (12)To extract the screening masses from the data, we employed the ﬁt function D − L,T ( p ) (cid:39) ζ E,M (cid:16) m E,M + p + r E,M ( p ) + q E,M ( p ) (cid:17) . (13)keeping terms up to p or ( p ) or ( p ) , depending on the ﬁt range. We will call below respective ﬁt functions as ﬁtfunction with two terms, with three terms or with four terms. We perform the ﬁt over the range p low ≤ p < p cut,where p low is equal to either zero or p min and p cut varies from 1.3 GeV to 2 GeV.We ﬁnd that, in the general case, both m E and m M depend substantially on the choice of p low, p cut and the ﬁtfunction. To improve stability of the results, we consider the propagators as functions of 2 variables, D − L,T ( p , µ I ).In so doing, we assume smooth dependence of a propagator on both variables and ﬁt the overall set of data by apolynomial in p and µ I thus reducing drastically the total number of estimated parameters. We choose the polynomialby selecting a ﬁnite number of terms in the series D − L,T ( p , µ I ) (cid:39) ∞ (cid:88) n =0 ∞ (cid:88) k =0 c n, k ( p ) n µ kI , (14)we omit the superscript L or T of the coeﬃcients c i,j . Odd powers of µ I are prohibited by isospin symmetry of stronginteractions and, in particular, gluon propagators are invariant under the change µ I → − µ I . Such polynomialsproviding a reliable quality of the ﬁt are found by trial-and-error method in the longitudinal case, the resultingexpressions for each ﬁt range are given below. The results of this (that is, combined) and the direct (that is, byformulas like 9) ﬁts for the inverse propagator at a given value of µ I are compared in Fig.1. It should be rememberedthat an increase of the number of estimated parameters gives rise to increasing of uncertainties in their values, whereasa decrease of this number results in either a dramatic decrease of the ﬁt p -value or in a substantial shrinking of theﬁt range. / D L ( p ) p datacombined fitthe direct fitthe combined fit µ I = 0.0 / D L ( p ) p datathe combined fitthe direct fitthe combined fit µ I = 1.0 FIG. 1: The inverse longitudinal propagator interpolated by formula (16) (the combined ﬁt) or by (9) and (13) (the direct ﬁt)at two diﬀerent values of µ I . The corridors of errors in the former case are shown by the respective error bars. Data are shownby purple squares. In the case of transverse propagator both ﬁtting procedures are not stable with respect to eliminating the zeromomentum; data on larger lattices and increased statistics are needed to arrive at reliable results.

A. Chromoelectric screening mass

We ﬁt our data by the function (13) with three or four terms using p cut = 1 . , . . p low = 0 or p min. The results of the ﬁts are stable with respect to exclusion of the zero momentum only for p cut = 2 GeV. Inthis case, the ﬁt function (13) with three terms works well at µ I ≥ . ≤ µ I ≤ . µ I are taken into account).As we explained above one can substantially reduce the number of ﬁt parameters by considering the inverse propa-gator as the function of p and µ I . This approach also makes it possible to derive m E as function of µ I and estimaterespective conﬁdence interval for each value of µ I ∈ (0 , . (cid:16) χ /

34 = 1 . p -value ≈ . (cid:17) can be achieved for the momentum range 0 ≤ p < . D − L ( p , µ I ) (cid:39) c + c µ I + c µ I + c µ I (15)+ c p + c p µ I + c p µ I + c p µ I + c p + c p µ I + c p µ I + c p + c p µ I + c p µ I + c p . m E , G e V µ I , GeV m E2 , fit Am E2 , fit B m E , G e V µ I , GeV m E2 , fit Am E2 , fit B FIG. 2: The chromoelectric squared screening masses as a function of the isospin chemical potential. Fit A is related to the ﬁtfunctions (9) and (13); ﬁt B—to (15); the momentum ranges are 0 ≤ p < . ≤ p < to lower correlations between the ﬁt parameters. The same ﬁt function provides the best ﬁt quality ( χ /

27 = 1 . p − value ≈ .

22) when the data at p = 0 are excluded: p min < p < . ≤ µ I ≤ . p cutresults in some decrease of the p -value ( χ /

55 = 1 . p − value ≈ . m E ( µ I ) = c + c µ I + c µ I + c µ I c + c µ I + c µ I + c µ I (16)The error in this quantity which takes into account the correlations between the ﬁt parameters was calculated withthe use of the REDUCE computer algebra system. The respective corridor of errors is shown by the error bars inFig.2.It should be emphasized that ﬁt B at both values of p cut and ﬁt A at p cut = 2 . p = 0 data. This indicates that the dynamics below p min has only a little eﬀect on the screening in thechromoelectric sector.We show the results of ﬁt A (13) with huge errors (for p cut = 1 . µ I ∼ . p cut = 2 . ∼ p is not taken into account when p cut = 2 . ∼ p and ∼ p are not taken into account when p cut = 1 . p > m E is a smooth function of µ I . The systematic error of ﬁt B canbe estimated by a comparison of the error corridors for the cutoﬀ momenta p cut = 1 . p cut = 2 . ≤ µ I ≤ . µ I at its greater values. The dependence of m E on µ I found here is qualitatively similar to itsdependence on the temperature at T > T c both in pure gluodynamics and QCD as was demonstrated by latticesimulations in [42, 44]. It is also similar to the dependence of m E on the baryon chemical potential in QC D [32, 40].

IV. DEPENDENCE ON THE MOMENTUM AND ISOSPIN CHEMICAL POTENTIAL

First we consider the momentum dependence of the gluon propagators for various values of µ I in more detail. Thepropagators are renormalized according to the MOM scheme to satisfy the condition D L,T ( p = κ ) = 1 /κ (17)at κ = 3 GeV. D L , G e V - p, GeV µ I =0.0 GeV µ I =0.2 GeV0.4 GeV0.6 GeV1.0 GeV1.2 GeV 1 10 0 0.5 1 1.5 2 D T , G e V - p, GeV µ I =0.0 GeV0.4 GeV1.0 GeV1.2 GeV FIG. 3: The propagators D L (left) and D T (right) as functions of p at diﬀerent values of µ I . The curves correspond to theformula (18). In Fig 3 (left) we present the momentum dependence for the longitudinal propagator D L ( p ) for some values of µ I .In the infrared domain it clearly decreases with an increase of µ I . This dependence on µ I is similar to the dependenceof D L ( p ) on T at T > T c in both gluodynamics and QCD [46, 47, 49] as well as on µ q in QC D [32].In Fig. 3(right) the momentum dependence for the transverse propagator D T ( p ) for the same values of µ I is shown. Itis clearly seen that D T ( p ) is substantially less sensitive to an increase of µ I . In the deep infrared region µ I dependenceof the transverse propagator diﬀers substantially from that of the longitudinal propagator. It increases with µ I atsmall µ I until it reaches a weakly pronounced peak (at µ I (cid:39) . ≤ p ≤ . µ I (cid:39) . p = 0 . p > µ I : 0 . ≤ µ I ≤ . µ I . A typical dependence of this type is shown in Fig.4 (right panel) by purple triangles; this behavior is also incontrast to that of the longitudinal propagator, which decreases with µ I at all momenta as is seen in the left panel ofFig.4It is known that at a ﬁnite temperature the propagator D T ( p ) has a clear maximum at the value of momentumincreasing with temperature. Our data give no evidence for such maximum at a small momentum, however, we donot rule out its existence.Now we proceed to an interpolation formula for our data. It was demonstrated many times [49–53] that the infraredbehavior of the gluon propagators at zero and ﬁnite temperature can be well described by the ﬁt function proposedin [54] as well as the tree level prediction of the Reﬁned Gribov-Zwanziger approach [55], D ( p ) = Z M + p p + M p + M . (18)This ﬁt function is in fact identical to that used in our previous study [32], the relations between the two sets of ﬁtparameters are obvious.Since the minimal nonzero momentum p min for the lattice under study is rather large our ﬁt over the infrareddomain may suﬀer from ﬁnite-volume eﬀects. Nevertheless, we believe that we ﬁnd qualitatively correct dependenceof the propagators on µ I , at least for D L at large µ I . D L , G e V - µ I , GeV p= 0.00 GeV 0.64 GeV 0.91 GeV 1.43 GeV D T , G e V - µ I , GeV p=0.00 GeV0.64 GeV0.91 GeV1.43 GeV FIG. 4: The propagators D L (left) and D T (right) as functions of µ I at diﬀerent values of p . Note logarithm scale on the yaxis. M , G e V µ I , GeV Dudal et al., p cut =1 GeVDudal et al., p cut =3 GeVM M M , G e V µ I , GeV Dudal et al., p cut =1 GeVDudal et al., p cut =3 GeVM M FIG. 5: Parameters M (left) and M (right) of the ﬁt (18) as functions of µ I . Yellow strips show the results obtained in [53]with p cut = 1 GeV, brown strips—with p cut = 3 GeV. To lower uncertainties in the ﬁt parameters M i , we perform ﬁt for the ratio D ( p ) D ( p ) (cid:39) M + p p + M p + M p + M p + M M + p . (19)where p should lie in the ﬁt domain; our choice is p = 2 GeV. In so doing, the normalization factors Z L,T are notdetermined and thus should be computed by the formula Z = D bare ( p ) κ D bare ( κ ) M + p p + M p + M . (20)With this Z factor, formula (18) gives the propagator renormalized at p = κ .The ﬁt domain in the longitudinal case ranges up to p cut = 5 GeV; in the transverse case the cutoﬀ momentum foreach value of µ I is indicated in Table II.A few words about the ﬁt stability should be said. In the longitudinal case, the ﬁt parameters are unaﬀected byelimination of zero momentum, in the transverse case the ﬁt parameters change signiﬁcantly at 4 of 7 values of µ I when zero momentum is excluded. M , G e V µ I , GeV Dudal et al., p cut =1 GeVDudal et al., p cut =3 GeVM M FIG. 6: Parameter M of the GS ﬁt (18) as a function of µ I . Yellow strip shows the result obtained in [53] with p cut = 1 GeV,brown strip—with p cut = 3 GeV. The ﬁt parameters for D L ( p ) depend only weakly on the cutoﬀ momentum, as is shown, as an example, in Table IIIfor µ I = 0; the range 2 . < p cut < . p cut the ﬁt parameters are poorlydetermined). At the other values of µ I the p cut -dependence of M i is also insigniﬁcant, however, the minimum valueof p cut at which the ﬁt parameters can be determined increases with µ I and reaches 4.5 GeV at µ I = 1 . D T ( p ) depend on p cut as follows. At µ I = 0 . p -value for 1 . < p cut < . M i depend on p cut only weakly. At µ I = 0 . p -valuefor 1 . < p cut < . M i depend signiﬁcantly on p cut at p cut > . ≤ µ I ≤ . p -value ≥ .

05 only at p cut ≤ . . ≤ p cut ≤ . M i vary within their errors and at p cut ≤ . M , M , and M are given in the respective powers of GeV for the sake of comparison with[53]. Results of this ﬁt are also shown in Fig.3 together with our lattice data.The dependence of these ﬁt parameters on the isospin chemical potential is presented in Figs.5 and 6, where ourresults are compared with those obtained in [53] in SU (3) gluodynamics at T = 0 on a large lattice with high statisticalprecision.The parameters M ( L )1 and M ( L )2 have a similar dependence on µ I : they jump between µ I = 0 . . µ I > . M ( L )3 shows approximately linear growth through the whole range µ I under consideration.In the case of the transverse propagator all M ( T ) i approach a shallow minimum at µ I ≈ . µ I > . M ( L ) i and M ( T ) i is seen at µ I = 0 . µ I ∼ . J L,T ( p ) deﬁned as J L,T ( p ) = p D L,T ( p ) (21)It is seen in Figure 7 (left) that with increasing µ I the maximum of the longitudinal dressing function goes downand shifts to the right, thus approaching dressing function of a massive scalar particle. We note once more that thisdependence on µ I is very similar to dependence on the temperature, see e.g. Ref. [47, 48].The transverse dressing function increases at µ I < . p ∼ . ÷ . µ I the height of the peak decreases with increasing µ I and the position shifts to p ∼ . ÷ . J = J T − J L is shown in Fig.8, rightpanel. It equals zero at µ I = 0 and then it rapidly increases over the range 0 < µ I < . µ I increases J L p, GeV µ I =0.0 GeV0.2 GeV0.4 GeV0.8 GeV1.2 GeV J T p, GeV µ I =0.0 GeV0.2 GeV0.4 GeV0.8 GeV1.2 GeV FIG. 7: The longitudinal (left)and transverse (right) gluon dressing functions at various isospin chemical potentials. The curvesshow results of the ﬁt by (18). from 0 . . < ∆ J < .

1, the width increases and the position shiftsto higher values of p .It should be emphasized that, at 0 < µ I < . J shows arapid change. V. FEATURES OF D L − D T BEHAVIOR

In the previous two sections we demonstrated that the diﬀerence between the longitudinal and transverse propa-gators increases with an increase of the isospin chemical potential. Therewith, they come close to each other at highmomenta for a ﬁxed value of µ I . In this section we study the rate of the decrease of ∆( p ) = D T ( p ) − D L ( p ) withincreasing momentum and how the picture changes with increasing µ I . A similar comparison of these two propa-gators was made in [32] in QC D as well as in [44] in SU (3) gluodynamics at ﬁnite temperatures where the ratio D L ( p ) /D T ( p ) was considered.Here we show that, in the theory under study, the diﬀerence between transverse and longitudinal propagators,∆( p ) = D T ( p ) − D L ( p ) at p = 0 shows clear exponential dependence on p .Our numerical results for ∆( p ) are presented in Fig.8(left), together with the ﬁt function∆( p ) = c exp( − ν · p ) . (22)The exponential decreasing is well established starting from some momentum p depending on µ I . We found that0 < p < p min for 0 . ≤ µ I ≤ . p ≈ . µ I = 1 . p < p < p cut is determined by the requirement that ∆( p ) diﬀers from zero by more than 3 statistical errors (seeTable IV).As a check we compare the ﬁt function (22) with the power ﬁt function∆( p ) = Cp v (23)motivated by a power-like behavior of both gluon propagators when p → ∞ .At 0 . ≤ µ I ≤ . p ) does not vanish for only a few momenta. At µ I > . p ) (cid:39) C ( p + d ) v ; (24)0 ∆ , G e V - p, GeV µ I = 0.2 GeV0.4 GeV0.6 GeV1.0 GeV ∆ J , b a r e p GeV µ I = 0.2 GeV0.4 GeV0.8 GeV1.2 GeV FIG. 8: The diﬀerence between the propagators D T − D L (left) and dressing functions J T − J L (right) as a function of themomentum. The curves show results of the ﬁt to (22). it also works only for 0 . ≤ µ I ≤ .

4, the respective p -values as well as the ﬁtting parameters in formula (22) arepresented in Table IV. We also tried several other power-like ﬁt functions which also have failed.Since our data can be well approximated by the exponential function and poorly — by a power-like ﬁt functioneven with 3 parameters, we arrive at the conclusion that the Gribov-Stingl ﬁt cannot work for D L ( p ) and D T ( p )simultaneously at large µ I . Assuming that they do work simultaneously, we obtain the power-like dependence ofthe diﬀerence D T ( p ) − D L ( p ), whereas, as we found, it cannot be ﬁtted by a power-like function. Of course, theabove reasoning is valid provided that the ﬁt domain is suﬃciently large and statistical precision is suﬃcientlyhigh to discriminate between power-like and exponential behavior. Given statistical precision as in the case underconsideration, power-like and exponential behavior are diﬀerentiated when the ﬁt domain is larger than p min ≤ p < . µ I ≥ . p ; µ I ) (cid:39) c exp( − ν · p ) cannot be ﬁtted by a power-like ﬁt function at µ I ≥ . D T ( p ) does not work at these values of µ I when p cut > . D L ( p ) even at p cut = 5 GeV for all µ I underconsideration.We show the curves resulting from our ﬁt (22) to ∆( p ) for some values of µ I in Fig.8 (left panel) as well as thediﬀerence between the dressing functions ∆ J ( p ) = p ∆( p ) = J T ( p ) − J L ( p ) (25)(right panel).The dependence of the parameters c and ν on the isospin chemical potential is shown in Fig. 9. The position andthe height of the peak of the diﬀerence between the dressing functions (25) is given by the formulas p max = 2 ν , ∆ J max = 4 e c ν , (26)the respective values are plotted in Fig.10.It is interesting to notice that, at µ I ≤ . J max rapidly increases, whereas p max does not change. A com-pletely diﬀerent situation occurs at µ I ≥ . J max changes only slightly, whereas p max increases substantiallyand shows a clear linear dependence on µ I at µ I ≥ . p max is similar to the behaviorof the electric screening mass (which is associated with the interaction radius), whereas ∆ J max (which is associatedwith the strength of interaction) changes in an opposite manner.1 c , G e V - µ I , GeV c ∆ (0) ν , G e V - µ I , GeV FIG. 9: The parameters of the ﬁt function (22) as functions of µ I . In the left panel the parameter c representing the value of∆(0) extrapolated by formula (22) from the ﬁt domain is compared with its actual value. p m ax , G e V µ I , GeV ∆ J m ax µ I , GeV FIG. 10: Dependence of the position (left panel) and the height (right panel) of the peak of ∆ J ( p ) on the isospin chemicalpotential. VI. CONCLUSIONS

We presented results of our study of the static longitudinal and transverse propagators in the Landau gauge of thelattice QCD with staggered quark action at nonzero isospin chemical potential.Our main observations are as follows. We found that the longitudinal propagator D L ( p ) is more and more suppressedin the infrared with increasing µ I . This is reﬂected, in particular, in the increase of the chromoelectric screening mass.Such dependence of D L ( p ) on µ I is analogous to its dependence on temperature at T > T c [46, 47, 49]. It is alsosimilar to dependence of D L ( p ) on the quark chemical potential µ q in QC D [32]. We found much weaker dependenceon µ I for the transverse propagator D T ( p ).We determined the parameters of the Gribov-Stingl ﬁt function (18) associated with the Gribov-Zwanziger conﬁne-ment scenario. The parameters for the transverse propagator depend only weakly on µ I and agree well with those2in pure gluodynamics at T = 0, the parameters M and M for the longitudinal propagator increase with increasing µ I . The diﬀerence between the parameters of the longitudinal and transverse propagators approaches a local peak at µ I = 0 . µ I ∼ . µ I the longitudinal dressing function is greater than the trans-verse one indicating that the chromoelectric interactions are screened at shorter distances than the chromomagneticinteractions.Finally, we studied the diﬀerence ∆( p ) = D L ( p ) − D T ( p ) and found that it decreases exponentially with themomentum at large p similarly to the case of the QC D theory [32]. We also investigated the µ I dependence of theparameters of the exponential ﬁt function (22) as well as the parameters p max and δJ max characterizing the behaviorof the diﬀerence between the longitudinal and transverse dressing functions. This dependence indicates that withincreasing µ I the domain where chromomagnetic forces dominate over chromoelecric ones extends to ever shorterdistances. Acknowledgments

The authors are grateful to V.V. Braguta and A.Yu. Kotov for very fruitful discussions. The work was completeddue to support of the Russian Foundation for Basic Research via grant 18-32-20172 mol-a-ved and grant 18-02-40130mega. The research was carried out using the Central Linux Cluster of the NRC ”Kurchatov Institute” - IHEP, theLinux Cluster of the NRC ”Kurchatov Institute” - ITEP (Moscow) and the equipment of the shared research facilitiesof HPC computing resources at Lomonosov Moscow State University. In addition, we used computer resources ofthe federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at NRC“Kurchatov Institute”, http://ckp.nrcki.ru/. [1] J. B. Kogut and D. K. Sinclair, Phys. Rev. D , 034505 (2002) [arXiv:hep-lat/0202028 [hep-lat]].[2] J. B. Kogut and D. K. Sinclair, Phys. Rev. D , 094501 (2004) [arXiv:hep-lat/0407027 [hep-lat]].[3] P. de Forcrand, M. A. Stephanov and U. Wenger, PoS LATTICE2007 , 237 (2007) [arXiv:0711.0023 [hep-lat]].[4] W. Detmold, K. Orginos and Z. Shi, Phys. Rev. D , 054507 (2012) [arXiv:1205.4224 [hep-lat]].[5] P. Cea, L. Cosmai, M. D’Elia, A. Papa and F. Sanﬁlippo, Phys. Rev. D , 094512 (2012) [arXiv:1202.5700 [hep-lat]].[6] G. Endr¨odi, Phys. Rev. D , no.9, 094501 (2014) [arXiv:1407.1216 [hep-lat]].[7] B. B. Brandt and G. Endrodi, PoS LATTICE (2016) 039;[8] B. B. Brandt, G. Endrodi and S. Schmalzbauer, Phys. Rev. D , no.5, 054514 (2018) [arXiv:1712.08190 [hep-lat]].[9] B. B. Brandt, G. Endrodi and S. Schmalzbauer, PoS Conﬁnement2018 , 260 (2018) [arXiv:1811.06004 [hep-lat]].[10] V. V. Braguta, A. Y. Kotov and A. A. Nikolaev, JETP Lett. , no.1, 1-4 (2019)[11] B. B. Brandt, F. Cuteri, G. Endr˝odi and S. Schmalzbauer, Particles , no.1, 80-86 (2020) [arXiv:1912.07451 [hep-lat]].[12] D. Toublan and J. B. Kogut, Phys. Lett. B , 212-216 (2003) [arXiv:hep-ph/0301183 [hep-ph]].[13] A. Barducci, R. Casalbuoni, G. Pettini and L. Ravagli, Phys. Rev. D , 096004 (2004) [arXiv:hep-ph/0402104 [hep-ph]].[14] L. y. He, M. Jin and P. f. Zhuang, Phys. Rev. D , 116001 (2005) [arXiv:hep-ph/0503272 [hep-ph]].[15] L. He and P. Zhuang, Phys. Lett. B , 93-101 (2005) [arXiv:hep-ph/0501024 [hep-ph]].[16] D. Ebert and K. G. Klimenko, J. Phys. G , 599-608 (2006) [arXiv:hep-ph/0507007 [hep-ph]].[17] G. f. Sun, L. He and P. Zhuang, Phys. Rev. D , 096004 (2007) [arXiv:hep-ph/0703159 [hep-ph]].[18] T. Xia, L. He and P. Zhuang, Phys. Rev. D , no.5, 056013 (2013) [arXiv:1307.4622 [hep-ph]].[19] T. G. Khunjua, K. G. Klimenko and R. N. Zhokhov, Phys. Rev. D , no.5, 054030 (2018) [arXiv:1804.01014 [hep-ph]].[20] T. Khunjua, K. Klimenko and R. Zhokhov, Symmetry , no.6, 778 (2019) [arXiv:1912.08635 [hep-ph]].[21] K. Kamikado, N. Strodthoﬀ, L. von Smekal and J. Wambach, Phys. Lett. B , 1044-1053 (2013) [arXiv:1207.0400[hep-ph]].[22] Z. Wang and P. Zhuang, Phys. Rev. D , no.1, 014006 (2017) [arXiv:1703.01035 [hep-ph]].[23] P. Adhikari, J. O. Andersen and P. Kneschke, Phys. Rev. D , no.7, 074016 (2018) [arXiv:1805.08599 [hep-ph]].[24] D. T. Son and M. A. Stephanov, Phys. Rev. Lett. , 592-595 (2001) [arXiv:hep-ph/0005225 [hep-ph]].[25] K. Splittorﬀ, D. T. Son and M. A. Stephanov, Phys. Rev. D , 016003 (2001) [arXiv:hep-ph/0012274 [hep-ph]].[26] P. Adhikari and J. O. Andersen, Phys. Lett. B , 135352 (2020) [arXiv:1909.01131 [hep-ph]].[27] J. O. Andersen, N. Haque, M. G. Mustafa and M. Strickland, Phys. Rev. D , no.5, 054045 (2016) [arXiv:1511.04660[hep-ph]].[28] T. Graf, J. Schaﬀner-Bielich and E. S. Fraga, Phys. Rev. D , no.8, 085030 (2016) [arXiv:1511.09457 [hep-ph]].[29] B. Klein, D. Toublan and J. J. M. Verbaarschot, Phys. Rev. D , 014009 (2003) [arXiv:hep-ph/0301143 [hep-ph]].[30] D. Toublan and J. B. Kogut, Phys. Lett. B , 129-136 (2005) [arXiv:hep-ph/0409310 [hep-ph]]. [31] M. Mannarelli, Particles , no.3, 411-443 (2019) [arXiv:1908.02042 [hep-ph]].[32] V. G. Bornyakov, V. V. Braguta, A. A. Nikolaev and R. N. Rogalyov, Phys. Rev. D (2020), 114511 [arXiv:2003.00232[hep-lat]].[33] S. Muroya, A. Nakamura, C. Nonaka, and T. Takaishi, Prog. Theor. Phys. , 615 (2003), hep-lat/0306031.[34] P. Weisz, Nucl. Phys. B212 , 1 (1983).[35] A. Bazavov et al., Phys. Rev.

D85 , 054503 (2012), 1111.1710.[36] J. E. Mandula and M. Ogilvie, Phys. Lett.

B185 , 127 (1987).[37] I. L. Bogolubsky, V. G. Bornyakov, G. Burgio, E. M. Ilgenfritz, M. Muller-Preussker and V. K. Mitrjushkin, Phys. Rev. D , 014504 (2008) [arXiv:0707.3611 [hep-lat]].[38] V. G. Bornyakov, V. K. Mitrjushkin, and M. Muller-Preussker, Phys. Rev. D81 , 054503 (2010), 0912.4475.[39] J. I. Kapusta and C. Gale,

Finite-temperature ﬁeld theory: Principles and applications , Cambridge Monographs on Math-ematical Physics (Cambridge University Press, 2011), ISBN 9780521173223, 9780521820820, 9780511222801.[40] V. Bornyakov, A. Kotov, A. Nikolaev and R. Rogalyov, Particles (2020) no.2, 308-319 [arXiv:1912.08529 [hep-lat]].[41] V. G. Bornyakov and V. K. Mitrjushkin, Int. J. Mod. Phys. A27 , 1250050 (2012), 1103.0442.[42] V. G. Bornyakov and V. K. Mitrjushkin, Phys. Rev.

D84 , 094503 (2011), 1011.4790.[43] O. Oliveira and P. Bicudo, J. Phys.

G38 , 045003 (2011), 1002.4151.[44] P. J. Silva, O. Oliveira, P. Bicudo, and N. Cardoso, Phys. Rev.

D89 , 074503 (2014), 1310.5629.[45] S. Ma,

Modern Theory of critical phenomena (W. A. Benjamin, Advanced Book Program, Minnesota University, 1976).[46] A. Cucchieri, A. Maas and T. Mendes, Phys. Rev. D (2007), 076003 doi:10.1103/PhysRevD.75.076003 [arXiv:hep-lat/0702022 [hep-lat]].[47] C. S. Fischer, A. Maas, and J. A. Muller, Eur.Phys.J. C68 , 165 (2010), 1003.1960.[48] A. Maas, Phys. Rept. , 203 (2013), 1106.3942.[49] R. Aouane, V. G. Bornyakov, E. M. Ilgenfritz, V. K. Mitrjushkin, M. Muller-Preussker, and A. Sternbeck, Phys. Rev.

D85 , 034501 (2012), 1108.1735.[50] D. Dudal, O. Oliveira, and N. Vandersickel, Phys. Rev.

D81 , 074505 (2010), 1002.2374.[51] A. Cucchieri, D. Dudal, T. Mendes, and N. Vandersickel, Phys. Rev.

D85 , 094513 (2012), 1111.2327.[52] O. Oliveira and P. J. Silva, Phys. Rev.

D86 , 114513 (2012), 1207.3029.[53] D. Dudal, O. Oliveira, and P. J. Silva, Annals Phys. , 351 (2018), 1803.02281.[54] M. Stingl, Z. Phys.

A353 , 423–445 (1996), hep-th/9502157.[55] D. Dudal, J. A. Gracey, S. P. Sorella, N. Vandersickel, and H. Verschelde, Phys. Rev.

D78 , 065047 (2008), 0806.4348.

Appendix A: Fit results

Parameters of the GS Fit to D renL µ I Z M M M p − value0.0 0.864(11) 2.55(18) 0.88(11) 0.165(14) 0.260.2 0.865(10) 2.36(23) 0.73(15) 0.334(26) 0.560.4 0.826(10) 5.86(95) 3.01(58) 0.86(11) 0.280.6 0.838(10) 5.0(1.0) 2.41(63) 1.84(24) 0.180.8 0.840(10) 5.9(1.1) 2.87(60) 4.03(57) 0.891.0 0.840(11) 7.7(1.8) 3.95(94) 7.3(1.3) 0.151.2 0.819(17) 13.5(3.3) 7.2(1.5) 16.3(3.0) 0.63TABLE I: Values of the parameters of the Gribov-Stingl ﬁt (18) for the longitudinal propagator D L ( p ) at various µ I . Thecutoﬀ momentum p cut = 5 GeV for all values of µ I Parameters of the GS Fit to D renT µ I Z M M M p − value p cut , GeV0.0 0.837(8) 2.96(16) 1.02(10) 0.189(11) 0.47 5.00.2 0.861(8) 2.01(14) 0.47(9) 0.125(9) 0.04 5.00.4 0.914(9) 1.36(14) 0.23(6) 0.087(8) 0.91 3.00.6 0.687(7) 4.7(1.3) 1.05(28) 0.259(47) 0.68 2.20.8 0.473(5) 8.9(4.9) 1.03(40) 0.326(79) 0.08 2.21.0 0.619(7) 6.4(2.4) 1.07(34) 0.381(78) 0.18 2.21.2 0.724(11) 6.2(5.0) 1.40(92) 0.72(31) 0.05 1.9TABLE II: Values of the parameters of the Gribov-Stingl ﬁt (18) for the transverse propagator D T ( p ). Stability of the GS ﬁt to D renL with respect to the momentum cutoﬀ p cut M M M p − value2.2 3.59(83) 1.20(27) 0.214(36) 0.922.6 2.67(50) 0.90(20) 0.172(26) 0.743.0 2.04(23) 0.67(11) 0.139(14) 0.863.5 2.10(18) 0.69(10) 0.143(13) 0.704.0 2.26(18) 0.75(11) 0.151(13) 0.394.5 2.43(19) 0.82(12) 0.159(14) 0.205.0 2.55(18) 0.88(12) 0.165(14) 0.25TABLE III: Dependence of the parameters M i , determined by a ﬁt of the formula (18) to D L ( p ) over the range 0 ≤ p < p cut on the cutoﬀ momentum p cut at µ I = 0. Results of the ﬁt of ∆( p ) by the exponential function µ I p cut , GeV c, GeV − ν, GeV − ∆ J max p − value p − value(formula (22)) (formula (24))0.2 1.3 17.9(10.4) 4.76(66) 0.52(21) 0.39 0.570.4 1.6 38.7(5.0) 4.839(80) 0.91(13) 0.96 0.430.6 2.2 17.4(2.1) 3.57(10) 0.738(76) 0.086 0.0020.8 2.5 19.3(1.6) 3.135(64) 1.066(65) 0.056 0.00041.0 3.1 15.07(72) 2.821(54) 1.014(64) 0.039 10 − c and ν as well as ∆ J max determined by ﬁtting formula (22) to our data for ∆( pp