Lee–Yang Zeroes in the Baryon Fugacity Plane: The Role of High Densities

: We compute the canonical partition functions and the Lee–Yang zeros in N f = 2 lattice QCD at temperature T = 1.20 T c lying above the Roberge–Weiss phase transition temperature T RW . The phase transition is characterized by the discontinuities in the baryon number density at speciﬁc values of imaginary baryon chemical potential. We further develop our method to compute the canonical partition functions using the asymptotic expression for respective integral. Then, we compute the Lee–Yang zeros and study their behavior in the limit of high baryon density.


Introduction
The properties of the QCD phase diagram in the temperature-baryon chemical potential plane are studied both experimentally and theoretically very intensively; see, e.g., [1][2][3][4].
The current experiments at RHIC (BNL) and LHC (CERN), and the future experiments at FAIR (GSI) and NICA (JINR), are devoted to such studies.Lattice QCD is a first principles theoretical approach to study the nonperturbative properties of QCD.It proved to be a very powerful method to study QCD at zero baryon chemical potential µ B .At nonzero µ B , the standard Monte Carlo methods are not applicable due to a severe sign problem.Various approaches to bypass this problem are under development with partial successes.These include reweighting, Taylor expansion, analytic continuation, canonical approach, strong coupling/dual methods, the density of states, and complex Langevin; see for review, e.g., [5,6].Thanks to recent developments, one can access the regions at finite T and µ B up to µ B /T ∼ 3; see, e.g., [7][8][9].But the task of locating the critical endpoint still remains unsolved.
In our work, we follow the canonical approach suggested in Ref. [10].For imaginary chemical potential, i.e., for µ B = iµ I , the sign problem is absent since the fermion determinant is positive.For this reason, it is possible to use standard Monte Carlo algorithms to compute (imaginary) baryon density ρ B and other observables; see, e.g., [11][12][13][14].Then, the canonical partition functions Z C (n, T, V) can be computed and, as a result, the grand partition function, Z GC , can be evaluated at any complex value of the chemical potential, including physically relevant real values.The properties of QCD at imaginary chemical potential were first formulated in Ref. [15].It was shown that Z GC is a periodic function of µ I and at high enough temperature (above T = T RW ) the free energy has cusps at µ I /T = π + 2πk, k = 0, ±1, ....This leads to discontinuities in the imaginary baryon number density signaling the first-order phase transition.
Discontinuities in the baryon-number density emerge at the points where the pressure p(θ) loses analyticity and, therefore, Z GC (ξ) vanishes.These points are accumulated on the lines, which tend to the cuts of the baryon-number density in the complex chemicalpotential plane; evidence for this scenario comes from the same arguments as in [16,17].
The values of ξ where Z GC (ξ) = 0 are called Lee-Yang zeros (LYZ); they were introduced in Refs.[16,17] by T.D. Lee and C.-N. Yang.They discussed the idea of how the non-analyticity of phase transitions can appear in the thermodynamic limit.The idea was to consider the truncated series of Z GC (θ, T, V) in terms of fugacity ξ B = e µ B /T = e θ .The truncated series becomes a polynomial whose zeros (Lee-Yang zeros) completely determine the phase structure of the theory.For this reason, the study of LYZ in QCD is an interesting and important task, and their properties were studied in a number of works by means of lattice QCD approach [18,19] or in effective QCD models [20][21][22].
After the pioneering work by Roberge and Weiss [15], where they pointed out the relation between imaginary chemical potential and the phases of QCD, the phase structure of the complex chemical-potential plane was thoroughly studied numerically [11], and the location of discontinuities of the baryon-number density as well as the consistency of the respective canonical partition functions with the experimental data were found [23].
The underlying reason for a particular arrangement of Lee-Yang zeroes in QCD and its dependence on the temperature is the subject of further studies.
In a number of papers [19,24,25], we developed methods to compute the canonical partition functions Z C (n, T, V) using the lattice data for imaginary baryon density ρ B at both low and high temperatures.In particular, in Ref. [25] using the stationary phase method we suggested the solution of the problem of computation of Z C (n, T, V) for temperatures above the Roberge-Weiss transition [15] when the imaginary baryon density has discontinuity and successfully applied this new method to the case of T = 1.35T c .At this temperature, the baryon density is described by the polynomial of degree three.
We demonstrated in Ref. [25] that Z C (n, T, V) computed using this method allow one to reproduce the input expression for the baryon number density and, moreover, allow one to obtain the correct Lee-Yang zeros pattern corresponding to the first-order Roberge-Weiss transition.Here, we improve our method using higher-order approximation for the stationary phase method to compute Z C (n, T, V) and apply it to the case of T = 1.20 T c > T RW when one needs to use a higher degree polynomial to describe lattice data for imaginary baryon density.
In this work, we study the arrangement of the Lee-Yang zeroes in the complex fugacity plane and the respective discontinuities of the net baryon-number density in the complex chemical-potential plane (along the lines θ I = 2πn, n ∈ Z) in more detail.
The paper is organized as follows.In Section 2, we describe the improvements in computation of Z C (n, T, V) and demonstrate their effect.In Section 3, the results for the Lee-Yang zeros are presented.We discuss the difference in their behavior compared to the higher temperature case.In Section 4, we present our conclusions.

Computation of the Canonical Partition Functions
It is useful to introduce notation Z n for the normalized canonical partition function, Z C (n, T, V)/Z C (0, T, V).The pressure p and the baryon number density ρ are defined as follows: where θ = µ B /T, N B is the net baryon number in the lattice volume.In the complex θ plane, it is convenient to use notations θ = θ R + ıθ I and ρ = ρ R + ıρ I .We also need the dimensionless variables ρ = ρ/T 3 = ρR + ı ρI .It is known that T > T RW ρ I (θ I ) can be well approximated by an odd-power polynomial [24,26].In this paper, we use numerical results for imaginary baryon number density at T = 1.20 T c > T RW .These numerical results are slightly different from those we presented before in Ref. [24] due to increased statistics.In Ref. [24], we fitted the lattice data for imaginary baryon number density by the polynomial of degree five For this fit, we find a 1 = 0.4901 (6), a 3 = 0.0132(3), a 5 = −0.00016(2).For real chemical potential, Equation (3) Negative a 5 implies unphysical behavior at large chemical potential if higher terms in (4) are absent.For this reason, we decided to also use another fit of the lattice data: We obtained positive a 7 = 1.43 × 10 −5 (13) for the fit (5) and slightly different a 1 and a 3 : a 1 = 0.49006(42), a 3 = 0.01360 (13).
For this fit, the real baryon number density in Equation ( 6) is positive for all values of θ R .Below, we present our results for the fit (5) and make a brief comparison with the results obtained for the fit (3).We show our data for the imaginary baryon number density, which fit (3) and ( 5) in Figure 1.The difference between the fits is not seen by eye.Let us note that the behavior of ρ I seen in Figure 1 implies that it has discontinuity at θ I = π since it is an odd periodic function of θ I with period 2π, as we discussed in the Introduction.Given coefficients a i , we can find the pressure up to the integration constant and, therefore, Z GC (θ, T, V) up to a factor.Then, we can compute the canonical partition functions by performing the Fourier transform numerically.
Before we proceed, we need to introduce some notations.We introduce the dimensionless quantity ν = VT 3 characterizing the volume (the expression for ν in terms of lattice parameters has the form ν = N 3 s N −3 t ).We also introduce the variable = n VT 3 , where n is the net baryon number of a particular physical state.This variable can help one to remember that the canonical partition functions Z C (n, T, V) correspond to the respective baryon number density .In what follows, is called a particular density, in contrast to the grand-canonical ensemble average density ρ.
In Ref. [24], we demonstrated that using Equations ( 1)-( 7) it is possible to compute Z n at temperature above T RW only for n below some n max .The value of n max is estimated from the following condition: the sign of the canonical partition functions Z C (n, T, V) calculated by Formula (7) alternates at n > n max .Negative values of Z C (n, T, V) (and, therefore, Z n ) are unphysical and indicate that our fit-function is not suitable.That is, it cannot describe high-density (n/ν > max = n max /ν) contributions to Z GC (θ, T, V).It should be noticed that n max depends significantly on the volume, whereas the respective density max depends on the volume only weakly and is in the region of 1.6 at T = 1.35T c .
In Ref. [25], we explained the reason for this problem and found its solution.The unphysical negative values of Z n stem from the unphysical condition that the finite-volume system under study undergoes during a phase transition, implicitly imposed by combining the 2π−periodicity of ρ I (θ I ) with the Equation (3).Actually, the function ( 3) is discontinuous at θ I = π + 2π k, k = ±1, ±2, ... .It should be emphasized that with Z n , n < n max evaluated in [24], we reproduced the dependence of the net baryon-number density ρI on θ I for all θ I except for small neighborhoods of the points θ I = 2πn, n ∈ Z.
For this reason, the right-hand side of Equation ( 3) has to be modified near the Roberge-Weiss transition at θ I = π.This modification should depend on the volume and involve discontinuities of the type (3) in the infinite-volume limit only.To evaluate the net baryon-number density ρ in a very close vicinity of the transition point θ I = π presents a challenge.In Ref. [25], we suggested another approach.We employed an approximate analytical expression for Z n originally introduced in [15].In this study, we further develop our ideas formulated in Ref. [25] and demonstrate that the improved approximation has the needed behavior near θ I = π and adequately reproduces density (3) in the infinite-volume limit.
It should be emphasized that the problem of negative canonical partition functions Z C (n, T, V) < 0 is associated with the asymptotic regime when ν is fixed, n → ∞ ( → ∞); however, a detailed study of high densities in the finite-volume case should be the subject of a separate study.Here, we consider two other asymptotic regimes: In the latter regime, the approach to obtain approximate expressions for the canonical partition functions in the case when a i = 0 for i > 3 was formulated in Ref. [15].Here, we consider the case when a 1 = 0, a 3 = 0, and a 7 = 0. Thus, we evaluate the quantities , where applying the stationary phase method [27].Let us note that the integrals in Equation ( 8) are computed along the imaginary axis, while the saddle point θ s is real, as is shown in Figure 2.For this reason, we deform contour AD to ABCD.In our case, we use the stationary phase method for ν → ∞.The saddle point θ = θ s = θ s ( ) is determined from the equation We note that dependence of θ s on is implied below though not explicitly indicated.
That is, we rearrange the integral in the numerator of (8) as follows: We use the substitution θ = θ s + ıθ I to rearrange the integral over BC into the second term in the lower row of formula (10), where the function f (θ I , ) has the form with Furthermore, we transform the integrals over AB and CD using, respectively, the substitutions θ = θ R − ıπ and θ = θ R + ıπ and obtain for functions h ± (θ R , ) The sum of the first and the third terms in Equation ( 10) should be omitted for the following reason.The baryon number density and, therefore, F n (θ) are holomorphic functions over the interior of the rectangle ABCD and, to employ the stationary phase method, we should keep them continuous over the rectangle ABCD with its boundary (the periodicity of the baryon-number density under consideration requires that discontinuities emerge along the lines AB and CD).Thus, we have to choose the values F n (θ) on the segments AB and CD as limiting values when θ tends to the boundaries AB and CD from the interior of the rectangle ABCD so that F n (θ) is continuous and thus non-periodic.This being so, F n (θ) = h + (θ R , ) on CD and F n (θ) = h − (θ R , ) on AB, where h ± (θ R , ) are given by formula (13).It is evident from Equation ( 13) that Re h This means that the first and the third terms in the upper row of formula (10) do not cancel each other, their sum is the product of the exponent exp(νRe h + (θ R , )) and the rapidly oscillating function 2ı sin(νIm h + (θ R , )).Thus, the sum of integrals The second term in Equation ( 10) can be estimated using the asymptotic expansion of the integral where in some neighborhood of the origin the function ϕ(x; λ) meets the Taylor expansion ϕ(x; λ) = ∑ ∞ n=0 c n (λ)x n , with coefficients c n (λ) remaining finite in the limit λ → ∞.The validity of Equation (15) for the estimate under consideration can be proven by replacing θ I = xν −1/3 .This substitution gives rise to where × exp 1 and the function ϕ from Equation ( 15) is provided by the second and third rows of Equation ( 17).Now, it is a straightforward matter to obtain the sought-for asymptotic expansion.

R(
where Finally, our approximation Z nA for the normalized canonical partition functions is Equation ( 20) corresponds to the next-to-next-to-leading approximation in the asymptotic expansion (16).It provides better approximation for Z nA than the values obtained in the leading approximation (first factor in Equation ( 20)) or in the next-to-leading approximation (first term in Equation ( 18)), as follows from the results for the baryon number density presented below.These approximations are denoted below as third-order, second-order, and first-order approximations.
To demonstrate the quality of the approximation Z nA , we computed the real and imaginary baryon number density using expressions and compared them with input densities ( 6) and ( 5).We present respective relative deviations in Figures 3-5.One can see from the left panels of Figures 3 and 4 that the relative deviation of ρI is substantially reduced when we improve the approximation for Z nA from the first approximation to the third approximation.We should note that we did not include the very vicinity of θ I = π in these figures.When θ I is very close to π, the functions ( 5) and ( 22) are intrinsically different for any finite volume.The agreement between the input ρI and approximated ρI becomes even better when we move from V = 16 3 to V = 40 3 , as follows from comparison of Figures 3 and 4.This effect is emphasized in the left panel of Figure 5.  5) and ( 22) for θ = ıθ I (left panel) and via Equations ( 6) and ( 21) for θ = θ R (right panel) for fixed volume V = 16 3 .Three approximations for Z nA are used in Equations ( 22) and ( 21), as explained after Equation (20).For the relative deviation of ρR shown in the right panels of same figures, we see qualitatively the same behavior with some differences.The second and third approximations provide more or less the same precision, which is much better than the precision of the first approximation.Increasing the volume also improves precision very much.
In the Figure 6, we show the approximated ρI in the vicinity of θ I = π to emphasize its strong dependence on the volume V.  22) near the θ I = π for three volumes V = 16 3 , 20 3 , 40 3 and best (third) approximation for Z nA .

Evaluation of Lee-Yang Zeros
The grand partition function Z GC (θ, T, V) in the finite volume V can be represented as a polynomial of the baryon fugacity where N tot = 2N f N 3 s L t for Wilson fermions.Zeros of Equation ( 23) are Lee-Yang zeros.The computation of zeros of the high degree polynomial, Equation (23), is in general a very difficult task.In this work, we employ a very efficient package, MPSolve v.3.1.8(Multiprecision Polynomial Solver) [28,29], which provides the calculation of polynomial roots with arbitrary precision.In our case, it was enough to use Z nA with an accuracy 12,000 digits to correctly find all LYZ.
Before discussing the results obtained for Equation (5), we would like to describe the observations made for Equation (3).We compare the results obtained for different volumes for the same values of max = N max /(VT 3 ).We start with max ≈ 0.39 (corresponding to N max = 25 for N s = 16) and increase it to max ≈ 5.47.The latter density value corresponds to the maximal value that we can obtain at T/T c = 1.20 with a 5 = 0 since respective Z nA become negative for a high enough density, which is related to the fact that ρR becomes negative for high density in the case of Equation (3).Examining LYZ distribution in detail, we observed a distribution of LYZ similar to that we reported in [25] for T/T c = 1.35.It should be noted that the first LYZ on the negative real axes appears at T/T c = 1.20 at max ≈ 1.47 in comparison with max ≈ 1.60 for T/T c = 1.35.It is also worth to note that for fixed N s , all LYZ located on the real negative axis for given value of N max will keep the same position when N max is increased, with new real LYZ appearing closer to ξ B = 0.This behavior of LYZ is observed for both temperatures T/T c = 1.35 and T/T c = 1.20.
The approximation with Equation (5) gives us the opportunity to study the behavior of the system in the limit of large values of N max since the density ρR is positive for all values of θ R .
In Figure 7, the results for LYZ are presented in the complex fugacity plane for two volumes with N s = 16, 40 and several values of N max for each volume.The general features of the LYZ distribution are very similar to those observed before in Ref. [25] for temperature T/T c = 1.35.In the bottom insertion of this figure, ξ cl show real LYZ nearest to ξ B = −1.
One can see that the positions of these LYZ do not depend on ρ max but depend on N s .The main difference from the results obtained for T/T c = 1.35 is the appearance of additional LYZ branches for max ≈ 8 for N s = 40 and a bit later for smaller volumes; they can be seen in the right insertion of Figure 7.The reason for the branches is not yet completely clear.We can see few such reasons: the insufficient numerical precision in the computation of LYZ, the approximation made in Equation (20) for Z nA , or the approximation for the baryon density in Equation ( 5).In Ref. [25], we derived the relation connecting the discontinuity of the baryon number density ∆ ρI to the density of LYZ on the real negative semiaxis: where g(θ R ) is the density of the LYZ at θ I = π.This relation is correct in the limit of an infinite volume.g(θ R ) is defined as where N LYZ (θ R ) is the number of LYZ at θ = θ R + iπ.We approximate g(θ R ) as where ∆θ R is the distance between adjacent LYZ on the real negative semiaxis.We found that Equation ( 24) is nicely satisfied at T/T c = 1.35 [25] when only terms up to θ 3 I were used in the approximation of the density ρI .Here, we check it for T/T c = 1.20 with ρI approximated by Equation (5).Using Equation ( 5) we derive for ∆ ρI Our results for the lhs and rhs of Equation ( 24) are presented in the Figure 8.One can see that the numerical value of the discontinuity ∆ ρI agrees quite well with the analytical prediction at small and large |θ R |.In the intermediate region, we observe a discrepancy with the analytical prediction in the form of oscillations.The maximum of this discrepancy is observed near |θ R | = 0.5.By comparing the results for two volumes, one can conclude that with increasing volume the amplitude of the oscillations is increasing but the range of |θ R | where they are observed is decreasing.Let us note that we observed similar behavior for the density g(θ R ), though less pronounced, for the case of Equation (3), but in the limit of large volumes such oscillations disappeared.The reason for these oscillations is still to be clarified.(24) between the LYZ density g(θ R ) defined in Equation ( 26) and the discontinuity ∆ ρI in the average particle-number density at θ I = π.
We also studied the behavior of complex LYZ forming small collapsing loop, see Figure 7.In Figure 9, we show the maximal real part of the complex LYZ for N s = 40 as a function of max for both T/T c = 1.2 and T/T c = 1.35.For both temperatures, this value goes to zero, indicating that complex LYZ are collapsing to a point ξ B = 0, but at T/T c = 1.2 this happens slower than for T/T c = 1.35.In Figure 9, we also show the fit by the function For both temperatures, this fit works quite well.We observed qualitatively similar features of other sides of the loop of complex LYZ, supporting the conclusion that they collapse to a point ξ B = 0 in the limit of large max .

Conclusions
We computed the canonical partition functions Z nA and LYZ in N f = 2 lattice QCD with Wilson fermion action in the deconfinement phase at T/T c = 1.20 above T RW and compared our new results with the results for T/T c = 1.35 presented in Ref. [25].In this paper, we further developed our method to compute the canonical partition functions.The stationary phase method used in computation of Z nA was extended to the third order, see Equation (20).
The method was applied to the case of the imaginary baryon density described by Equation (5), with coefficients a i determined by fits of the lattice results as shown in Figure 1.We demonstrated that this method to compute Z nA works correctly by reconstructing the input baryon densities.We observed that the relative deviation δ ρI / ρI of the approximated ρI from the input ρI becomes substantially better when we improve the approximation for Z nA from the first approximation to the second approximation and then to the third approximation.This agreement improves further when we increase the spatial volume from V = 16 3 to V = 40 3 , as can be seen in the left panel of Figure 5.For the respective relative deviation of ρR shown in the right panels of same figures, we found that the second and the third approximations provide qualitatively similar precision, which is much better than the precision of the first approximation.In this case, the spatial volume increasing also improves precision quite substantially.
Using the canonical partition functions Z nA , we computed LYZ.We found that the general features of LYZ are the same as at T/T c = 1.35.Most importantly, there are LYZ on the negative real semi-axes, signaling the first-order Roberge-Weiss transition.In Figure 8, we depict LYZ density g(θ R ) defined in Equation ( 26) and the discontinuity ∆ ρI in the average particle-number density at θ I = π and show that the derived relation (24) is satisfied.This figure confirms the reliability of our numerical results for LYZ.In this work, we completed computations at very high densities max , reaching values 8.0 (for N s = 40) or even 16 (for N s = 16), and observed some deviations from regular behavior in the properties of LYZ, indicating the limits of applicability of our method.
We think that the approach to compute the canonical partition via the stationary phase method, which we presented in Ref. [25] and developed further in this paper, as well as the results obtained for Z nA and LYZ at high temperatures, will be useful for the computation of these important quantities at lower temperatures where a phase transition or cross-over is present at real values of the chemical potential.

Figure 1 .
Figure 1.Lattice numerical data for imaginary baryon number density ρI (θ I ) at T = 1.20 T c and fitting functions Equations (3) and (5).The coefficients a i are presented in the text below the corresponding explicit expressions (3) and (5).

Figure 2 .
Figure 2. Integration contour for evaluating values of Z nA (8) in the complex θ plane.Explicit partitioning into the sum of integrals over different parts is presented in expression (10) .
imaginary and can be omitted since the integrals both in the numerator and denominator of formula (8) are real.Therefore, the sum B A + D C just cancels the imaginary part of the integral C B .

Figure 3 .
Figure 3. Relative deviation of the baryon number density evaluated via Equations (5) and (22) for

Figure 4 .
Figure 4. Relative deviation of the baryon number density as in Figure 3 but for volume V = 40 3 .

Figure 6 .
Figure 6.Behavior of the baryon number density for the imaginary chemical potential computed via Equation (22) near the θ I = π for three volumes V = 16 3 , 203 , 403 and best (third) approximation for Z nA .

Figure
Figure LYZ distribution in the fugacity plane at T/T c = 1.20.ξ cl in the bottom insertion show real LYZ nearest to ξ B = −1.In the right insertion, we show the vicinity of ξ B = 0.

Figure 9 .
Figure 9. Behavior of the real value of the most right complex LYZ located on a small loop near ξ B = 0 as a function of the maximum density max for temperatures T/T c = 1.35 and 1.20.The lines show results of the fit to the function (28).