Lee-Yang zeros in lattice QCD for searching phase transition points

We report Lee-Yang zeros behavior at finite temperature and density. The quark number densities,, are calculated at the pure imaginary chemical potential, where no sign problem occurs. Then, the canonical partition functions, Z_C(n,T,V), up to some maximal values of n are estimated through fitting theoretically motivated functions to, which are used to compute the Lee-Yang zeros. We study the temperature dependence of the distributions of the Lee-Yang zeros around the pseudo-critical temperature region T/T_c = 0.84 - 1.35. In the distributions of the Lee-Yang zeros, we observe the Roberge-Weiss phase transition at T/T_c>= 1.20. We discuss the dependence of the behaviors of Lee-Yang zeros on the maximal value of n, so that we can estimate a reliable infinite volume limit.


Introduction
Revealing the phase structure at finite temperature and density in quantum chromodynamics (QCD) is one of the most important and interesting subjects in quark-hadron physics. At high temperature and low density the quarkgluon plasma was created by heavy ion colliders at RHIC (BNL) [1] and LHC (CERN) [2]. High density region of the QCD phase diagram will be explored by FAIR (GSI), NICA (JINR), and J-PARC (KEK/JAEA) in the near future. These experiments will provide us valuable information for understanding the early universe and the interior of neutron stars.
On the theoretical side, the lattice QCD is an ideal tool to conduct directly non-perturbative calculations in QCD. The lattice QCD, however, suffers from the sign problem at finite density: the fermion determinant det D(µ q ) at finite quark chemical potential µ q is complex in general, and consequently, it is impossible to apply the conventional Monte Carlo method.
The canonical approach [3], which we use in this paper, is a promising candidate for avoiding the sign problem. Since the fermion determinant at pure imaginary chemical potential µ q = iµ qI (µ qI ∈ R) is real, canonical partition functions Z C (n, T, V ) can be evaluated with the conventional Monte Carlo method at finite µ qI . Note that the canonical partition functions Z C (n, T, V ) depend on a net number of quarks and antiquarks, n, tempera-ture T and volume V , but do not depend on the chemical potential. Although the difficulty associated with the complex fermion determinant remains as the highly oscillating integral, authors of Refs. [4][5][6] pointed out that the integral can be carried out with a multi-precision arithmetic. Once the canonical partition functions are extracted, the grand canonical partition function Z GC (µ q , T, V ) with an arbitrary chemical potential can be constructed via the fugacity expansion, where ξ q = e µq/T is the quark fugacity. The canonical approach has been used in several analyses to reveal the QCD phase diagram [5][6][7][8][9][10][11][12][13][14][15][16].
In this paper, we employ Lee-Yang zeros (LYZs) to clarify the QCD phase structure. Analysis of Lee-Yang zeros (LYZs) is a general and powerful tool to investigate the phase structure. The theorems of Yang and Lee [17,18] state that zeros of the grand canonical partition function in complex fugacity plane contain a valuable information on the phase transitions of a system.
The number of LYZs is finite in a finite system. In the infinite size limit of the system, the number becomes infinite and one-dimensional cuts are made by the zeros, and the edges of the cuts, which are the accumulation points of the LYZs, are going to the phase transition point in the fugacity. By investigating the behavior of the edges with system size, we can discuss the properties of phase transition.
In the literatures, there are several studies of phase structure at finite density in lattice QCD by using LYZs. The pioneering work is that by Barbour and Bell [19], which were performed on lattice volumes, 2 4 and 4 4 . Then the strong coupling calculation at zero temperature using the Glasgow algorithm [20], and one dimensional QCD simulations [21] were reported. Fodor and Katz studies QCD phase diagram by using LYZs on larger lattice sizes. They determined the critical point by distinguishing the first order transition and crossover with the multiparameter reweighting method on realistic quark masses, N t = 4, and N s = 6, 8, 10, and 12 lattices, where N t (N s ) is the number of lattice sites in temporal (spatial) direction [22]. There was, however, a discussion of the reliability of their results [23].
In Ref. [24], LYZs are discussed both in analytic and theoretical approaches; the authors found that the scaling of LYZs with the system size in lattice QCD around the Roberge-Weiss (RW) phase transition point [25] is consistent with that derived analytically in the high temperature limit. In this paper, We will discuss on the behavior of LYZs around RW phase transition point. Recently, LYZs were estimated by the random matrix model [26].
As shown in Ref. [27], the distributions of LYZs can be calculated from the experimental data of the net-proton multiplicity distribution at RHIC [28,29]. Therefore, studies of LYZs of lattice QCD will be useful to compare with ones from data of future experiments.
In this paper, we show the first results of the calculation of LYZs by the canonical approach. We discuss the dependence of distributions of LYZs at temperatures T = 0.84-1.35 on N t = 4 and N s = 8, 16, and 20 lattices with quark mass m π /m ρ = 0.80.
It is difficult to calculate the grand canonical partition function Z GC directly. Therefore we calculated Z GC by the "integration method" [13][14][15][16], in which Z GC is obtained by integrating a fit function of baryon number densities n B in the pure imaginary chemical potential regions. See Sec. 2.
The method has the essential difference in calculating Z n from previous works. Since the number density is fit by some analytic functions and and we construct Z GC from that, Z n can be calculable for infinitely large n in principle.
There is, however, a practical limit of the maximum number of Z n , N max , which is smaller than the size of the system. Consequently, we have to see the scaling of LYZs not only with the system size but also with N max .
The main objective of the present study is to investigate the distributions of LYZs for various temperatures from canonical method, and to extract information of the phase transition in real chemical potential region from the method.
The organization of the paper is as follows. In Sec. 2, the canonical approach is briefly described. In Sec. 3, algorithms used in calculations of LYZs are introduced. Here "cut Baum-Kuchen algorithm" is explained, which is used for finding the roots of the grand canonical partition function. In Sec. 4, the numerical results are shown. We show the properties of distributions of LYZs for temperatures T = 0.84-1.35, and their behavior is examined for each temperature. Section 5 is devoted to the summary.

Canonical approach
The canonical partition functions in Eq. (1) can be written as a Fourier transformation of the grand canonical partition function of the pure imaginary chemical potential [3], where θ q = µqI T . We can construct Z GC at the pure imaginary µ q as where C is an integration constant and n qI is an imaginary number density, n q = in qI . The number density n q is defined as where S G is a gauge action, D(µ q ) is a Dirac operator and N f is the number of flavors. The conventional Monte Carlo method is applicable to a calculation of n q at the pure imaginary µ q since the fermion determinant is real. It was observed that the pure imaginary number density n qI can be well approximated by an odd power polynomial of θ q , at high temperature and by a Fourier series, at low temperature [30][31][32][33][34]. In this paper, we employ Eqs. (5) and (6) with maximal values N poly and N sin . Note that when we fit Eq. (6) to the number density, the integration in Eq. (3) can be performed as 3k and C is an integration constant. Complex Fourier series of the parts of Z GC can be written with the modified Bessel function of the first kind, I n (x): Therefore, the grand canonical partition function follows the equation where all Z C (n) for mod(n, 3) = 0 are zero. In the case of N sin = 1, it corresponds to the Skellam model [35]: that is often used in analyzing energy heavy ion collision experiments.

Partition functions and Lee-Yang zeros
The grand canonical partition function can be rewritten as where we took into account the property of the canonical partition function Z C (n, T, V ) = 0 for mod(n, 3) = 0 because of the Roberge-Weiss symmetry on the number density. Here µ B = 3µ q is the baryon chemical potential and ξ B = ξ 3 q is the baryon fugacity. We truncate the fugacity expansion at |n| = N max . On an N 3 s × N t lattice at finite temperature, where N s (N t ) is the number of lattice sites in spatial (temporal) direction, the exact grand canonical partition function is obtained for N max = 2N 3 s N f . Lee-Yang zeros, α n , are given as roots of the equation, in the complex ξ B plane. We must solve a polynomial equation of high degree of 2N max . This is a famous illposed problem. In Ref. [27], a new algorithm was proposed to overcome the difficulty: the cut Baum-Kuchen (cBK) 1 1 Figure 1: Schematic diagram to show paths in the cut Baum-Kuchen algorithm, taken from Ref. [27]. We perform the Cauchy's integration, Eq. (17), along a closed path, which gives the number of LYZs inside of the path. If the outcome is not zero, we divide the path into four pieces and continue the integration on each path.
algorithm with the multi-precision arithmetic. The multiprecision arithmetic is done with the FMLIB package [36] with 100-300 significant digits. The cBK algorithm is as follows. Since the α n are roots of f (ξ B ), It is easy to see that f (ξ B ) satisfies Then, we can extract the number of LYZs inside a closedintegral path C 0 by Cauchy's integral, Z C (n, T, V ) satisfy from the charge-parity invariance. If α n is a LYZ, then α * n and α −1 n are also LYZs because the Z C (n, T, V ) are real and satisfy Eq. (18). Thanks to these properties, it is enough to search only inside the upper half of the unit circle in the complex ξ B plane. The first closed-integral path C 0 is along the boundary of the plane shape defined by relations 0 ≤ arg (ξ B ) ≤ π and 0 ≤ |ξ B | ≤ 1. If there are LYZs inside a contour, we divide respective annulus sector into four parts bisecting both radial and angular ranges of polar coordinates and perform the Cauchy integrals for each sector. In the cBK algorithm we can locate LYZs by conducting this procedure recursively. See Fig. 1.

Lattice setup
We generate the gauge field configurations in full QCD using the hybrid Monte Carlo method [37] with the Iwasaki gauge field action [38] and the N f = 2 clover improved Wilson fermion action [39]. Simulations are carried out on   ) lattices to check the volume dependence. We use the same parameters (couplings β and hopping parameters κ) as those in Ref. [40]. The pseudo-critical temperature T c was determined from the peak of the Polyakov-loop susceptibility [40]. All parameters on a N 3 s = 16 3 lattice are listed in Table 1. The clover coefficient c SW in the clover improved Wilson fermion action is given by c SW = (1−0.8412β −1 ) −3/4 [38]. We produce 4000 (2000) gauge field configurations taking every 10th trajectory for T /T c = 0.84, 0.93, 0.99, and 1.08 (T /T c = 1.20 and 1.35). The first 200 configurations are discarded for thermalization.

Lee-Yang zeros in the complex ξ B plane
We calculate the number density n qI /T 3 at 19 to 40 values of µ qI depending on temperature. We fit the n qI /T 3 with Eq. (6) with N sin = 1, 1, 3, and 7 for T /T c = 0.84, 0.93, 0.99, and 1.08, and with polynomials Eq. (5) with N poly = 5 and 2 for T /T c = 1.20 and 1.35, respectively. Coefficients, f k and a k , obtained by the fitting are shown in Table 1.
As shown in Table 2, we calculate LYZs with the cBK algorithm for few values of N max for each temperature value. After the cBK recursions are carried out with eight or more times, we can obtain magnitude r B = |ξ B | and argument θ B = arg (ξ B ) of LYZs at centers of obtained annulus sectors. Due to the finite size of the annulus sectors, the LYZ coordinates have systematic errors: δr B < 2.0 × 10 −3 ∼ 1/2 9 and δθ B < 6.2 × 10 −3 ∼ π/2 9 .

Temperature dependence
In Fig. 2, we present the temperature dependence of LYZs. We show in Figs. 2, 3, 6, and 7 only part of the LYZs. The rest of them can be restored using complex conjugation or inversion; See the previous section. For all temperatures except T /T c = 1.20, LYZs are calculated with N max = 100. The calculation of LYZs for T /T c = 1.20 is terminated at N max = 95 because for T /T c = 1.20 negative values of Z C (3n, T, V ) appear for 3n ≥ 294, see also Ref. [16]. Note that all LYZs have the systematic errors within symbols although the errors are not displayed in Fig. 2. The same for the following figures.
Distributions of LYZs above T c and below T c are quite different. Below T c the distributions of the LYZs are nearly circles. And as the temperature becomes lower, the radius of the distribution becomes smaller for the fixed N max .
On the other hand, above T c the distributions are reaching to the value ξ B = −1. In the complex µ q /T plane, this where k is an integer. Thus, ξ B = −1 represents the RW phase transition. At T /T c = 1.35, 1.20, and 1.08, there are LYZs that are very close to ξ B = −1. We do not see, however, a LYZ exactly at ξ B = −1. This is probably due to finite V and finite N max effects.

V and N max dependences
Phase transitions of QCD shall be seen in large-V and large-N max limits. Consequently, it is important whether a LYZ appears on the real positive axis in the complex ξ B plane as V and N max go to large values.
In Fig. 3, we show the V and N max dependences of LYZs at T /T c = 1.35 in the complex ξ B plane. Following Ref. [27], we choose N max so that the value of N max /V is approximately same when V changes. It means that we fix number of baryons that can exist in the system per volume. Then, we find that LYZs fall onto the almost same curve if V is sufficiently large as 16 3 or 20 3 .   There are some LYZs on the unit circle at small-N max . However, they go away from the unit circle when N max becomes larger: thus these LYZs are an artifact due to the small N max . As can be seen in Fig. 5, when we increase N max , the right edge of LYZs goes to zero, where the right edge is determined by a position of the LYZ for Re [ξ B ] > 0 and min (Im [ξ B ]). Therefore, there is no phase transition for real µ B /T at T /T c = 1.35.
On the contrary, an important observation is that the LYZ nearest to ξ B = −1 approaches this point as V increases and the LYZ is stable with respect to changing of N max ; See Fig. 4. This tells us that there is a RW phase transition at T /T c = 1.35.
We calculated LYZs also for N max ≥ 125. But the results are polluted by negative values of Z C (3n, T, V ) appearing from 3n = 312. As N max increases beyond N max = 125, the right edge of LYZs suddenly turns away from zero in spite of approaching zero until N max = 100, see Fig. 5. The distributions of the LYZs are sensitive to the inclusion of the negative Z C .
Next, we investigate the N max dependence of LYZs in the confinement phase. In Fig. 6, the result at T /T c = 0.93 is presented. Simulations are carried out up to N s = 16 because we find that the N s = 16 lattice is sufficiently large in Fig. 3. In the following figures, we only show results for N s = 16 and focus the discussion on the N max dependence. In Fig. 6, as N max increases, a right edge of LYZs approaches to the real positive axis. We extrapolate the right edges of LYZs to the real axis by linear or quadratic functions and a phase transition point can be roughly estimated as µ B /T ∼ 5-6 at T /T c = 0.93. At T /T c = 0.99, we also estimate a phase transition point µ B /T ∼ 3-3.5 with the same determination. We should point out that the determination of a phase transition strongly depends on a fitting function. Thus, we do not mention an estimation at T /T c = 0.84 because right edges of LYZs at T /T c = 0.84 are too close to zero: min(|ξ B |) = 0.0051.
The N max dependence of LYZs at T /T c = 1.08 in the complex ξ B plane is shown in Fig. 7. At T /T c = 1.08, the Z C (3n, T, V ) for all n are positive. Because LYZs near ξ B = −1 do not exist for N max larger than or equals to 125, there is no RW phase transition at T /T c = 1.08. The right edge of LYZs approaches the real positive axis as N max increases. However, the right edge is still far from zero even for N max = 300; This suggests that there might be a stable point around µ B /T ∼ (1.03,0.08).

Lee-Yang zeros in the complex µ B /T plane
Let us study the distributions of LYZs in the complex µ B /T plane. This includes the same information as in the complex fugacity plane, but allows us to see LYZ structure from different perspective. Fig. 8 shows the results of the LYZs in the complex µ B /T plane at T /T c = 1.35. The LYZs shown in Fig. 3 correspond to ones in the second quadrant in Fig. 8. Since LYZs have the properties discussed above, LYZs also exist in the first quadrant. We calculated the LYZs also at N max = 125, 150, and 200 but they suffered from the negative Z C .
We find that the LYZs far from the imaginary axis fall onto approximately a straight line for each N max . In other words, the LYZs in the upper unit circle in Fig. 3 locate on a curved line, with a parameter t and two constants p and q.

Summary
We studied LYZs, i.e., zeros of the grand canonical partition function, which we calculated from the number density using the integration method. The number density at the pure imaginary chemical potential were evaluated numerically in two-flavor full QCD. We used the Iwasaki gauge action and the clover improved Wilson fermion action. Simulations were carried out on N t = 4 and N s = 8, 16, and 20 lattices at m π /m ρ = 0.80 for temperatures in the range T /T c = 0.84 -1. 35.
We obtained the canonical partition functions, Z C (n, T, V ), after fitting the number density with Fourier series with N sin = 1, 1, 3, and 7 for T /T c = 0.84, 0.93, 0.99, and 1.08, and with polynomials with N poly = 5 and 2 for T /T c = 1.20 and 1.35.
We obtained LYZs with the cut Baum-Kuchen algorithm for each temperature at few N max values. The Roberge-Weiss phase transition clearly appears in the distribution of LYZs at T /T c = 1.35 and 1.20. Additionally, in the complex µ B /T plane, the LYZs far from the imaginary axis fall onto approximately straight line for each N max at T /T c = 1.35 and 1.20.
We estimated phase transitions in the confinement phase as µ B /T ∼ 3-3.5 at T /T c = 0.99 and µ B /T ∼ 5-6 at T /T c = 0.93. More precise values at these temperature and an estimation at T /T c = 0.84 need calculations with larger N max .
Moreover, in the deconfinement phase near T c (T /T c = 1.08), it is likely that there is a stable point around µ B /T ∼ (1.03,0.08) with respect to changing of N max . To confirm this, we need study more extensively the relation not only between the stable point and N max , but also between the point and N sin below T /T c = 1.08.
Above T /T c = 1.20, the negative Z C (n, T, V ) appear for large-n, which should not be used in the search of phase transitions from LYZs. Therefore, our searches of the N max dependence are limited in some N max . To avoid the negative Z C , we need more statistics and extract fitting coefficients a 2k−1 of a polynomial more accurate.
acknowledgments One of the authors (M.W.) wishes to thank M. Sekiguchi for discussions and encouragement. This work was completed due to support by the Russian Science Foundation Grant under Contract No. 15-12-20008. Work done by A. Nakamura on the theoretical formulation of the Z n was supported by JSPS KAKENHI Grant Numbers 26610072 and 15H03663. Computer simulations were performed on the GPU cluster Vostok-1 of Far Eastern Federal University. This research used computational resources of Cybermedia Center of Osaka University through the HPCI System (ID:hp170197). This work was supported by "Joint Usage/Research Center for Interdisciplinary Large-scale Information Infrastructures" in Japan (ID: EX18705).