Net baryon number fluctuations across the chiral phase transition at finite density in the strong coupling lattice QCD

We investigate the net-baryon number fluctuations across the chiral phase transition at finite density in the strong coupling and chiral limit. Mesonic field fluctuations are taken into account by using the auxiliary field Monte-Carlo method. We find that the higher-order cumulant ratios, $S\sigma$ and $\kappa\sigma^2$, show oscillatory behavior around the phase boundary at $\mu/T\gtrsim 0.2$, and there exists the region where the higher-order cumulant ratios are negative. The negative region of $\kappa\sigma^2$ is found to shrink with increasing lattice size. This behavior agrees with the expectations from the scaling analysis.


Introduction
Fluctuations of conserved charges are promising observables in relativistic heavy-ion collisions in search of the QCD phase transition [1][2][3][4][5]. In the first phase of the Beam Energy Scan (BES I) program at the Relativistic Heavy Ion Collider (RHIC), net-proton [6,7], as a proxy for net-baryon [8,9], and net-electric charge event-by-event fluctuations [10] have been measured in Au+Au collisions for a broad energy range from √ s N N = 7.7 to 200 GeV. The success of the statistical thermal model [11] indicates that the event-by-event fluctuations of the multiplicity of conserved charges can be regarded as particle number fluctuations in the grand canonical ensemble specified by temperature, volume, and chemical potential. 1 The particle number fluctuations of the conserved charge, i.e., susceptibilities or cumulants, reflect the properties of the phase in QCD [13], which is expected to exhibit a rich structure [14]. In particular, if QCD has a critical point (CP) [15], the susceptibilities and higher-order cumulants diverge at the critical point owing to the divergent correlation length [16,17]; thus, one expects that anomalously large fluctuations would be observed if the system passes around CP, as shown in Fig. 1. In heavy-ion collisions, the correlation length cannot diverge due to the finite size of the system, but these fluctuations are expected to remain sensitive to the remnant criticality of the system around Schematic QCD phase diagram. In the chiral limit, we expect that there exist the second-order phase transition at low μ and the first-order phase transition at high μ. The two phase transition lines are connected at the tricritical point (TCP). Around TCP and the second-order phase transition line, anomalous κσ 2 could be realized in the chiral limit. At finite mass, the transition at low μ becomes a crossover and the first-order phase transition at high μ ends at the Z(2) critical point (CP). An anomalous κσ 2 region at finite mass around Z(2) CP is also expected.
CP [1,2]. Specifically, such criticality can be more pronounced in higher-order cumulants, whose accurate measurements require a high-precision tail of the event-by-event multiplicity distribution. In Refs. [18][19][20], it was shown that the tail of the multiplicity distribution must be properly included to obtain the correct value of the higher-order cumulants.
The measured higher-order cumulant ratios of net-proton number, Sσ = χ (3) p /χ (2) p and κσ 2 = χ (4) p /χ (2) p , show non-monotonic behavior as a function of the incident energy, where σ 2 , S, and κ are referred to as the variance, skewness, and kurtosis, respectively. At around √ s N N = 20 GeV, the data show κσ 2 < 1 below the expected value in the Skellam distribution [6,7]. The decrease of κσ 2 in the net-proton number can be the signal of critical behavior. According to theoretical arguments [21,22], κσ 2 of the net-baryon number can be reduced by critical behavior from universality. In QCD, the expected universality class depends on the quark masses. For twoflavor massless quarks with axial anomaly, the QCD phase transition at finite T (μ = 0) belongs to the universality class of the 3D O(4) symmetric spin model [23][24][25][26], and the second-order transition at low μ may turn into the first-order transition at the tricritical point (TCP). At physical quark masses, a finite-T phase transition would be a crossover [27,28], and the fluctuations are governed by the approximate chiral symmetry. The pseudo-critical line at low μ may be connected with the first-order transition line at CP, whose criticality is expected to belong to the Z(2) universality class [15,16,29]. A schematic QCD phase diagram is shown in Fig. 1.
While the universality argument dictates the singular behavior of the thermodynamic quantities close to the phase transition, the actual magnitude of the fluctuations is smeared by the finite-volume effect in addition to the finite quark mass. Finite volume makes the transition smoother. We need to perform a finite-size scaling analysis to observe the precise critical behavior [27,28]. The critical behavior depends on the relative strength of the singular part to the regular (non-singular) part of the free energy density (or its derivatives). The regular part would also mask the expected singular behavior from criticality [21]. Thus, explicit calculations are desirable to see how much criticality can be present in the fluctuation observables, and to pin down the origin of the observed decrease of κσ 2 : Z(2) critical behavior around CP [17], the remnant of the O(4) chiral phase transition [21], or other mechanisms. Lattice QCD (LQCD) is the most powerful non-perturbative approach to QCD based on Monte Carlo (MC) simulations, and is successful at vanishing or low baryon chemical potential. For instance, higher-order cumulants have been studied at zero chemical potential [30][31][32][33][34][35][36][37] and at finite chemical potential [37][38][39][40][41][42][43]. At large baryon chemical potential, however, LQCD faces the notorious sign problem, which makes it difficult to carry out MC simulations owing to the complex fermion determinant. There have been many attempts to circumvent the sign problem [25,37,[44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59]. Most of the methods of evading the sign problem are reliable only in limited circumstances such as μ/T 1, small volume, or heavier quark masses than the physical one, depending on the method. Conclusive results at the physical point have not yet been obtained.
One possible alternative way to attack the finite density region is the strong-coupling approach of LQCD with fermions . The strong-coupling approach with fermions could reduce the severity of the sign problem and is applied to investigate the chiral phase transition. Recently, theoretical frameworks including fluctuation effects in the strong-coupling limit (SCL) have been developed: the auxiliary field Monte Carlo (AFMC) method [69] and the monomer-dimer-polymer (MDP) simulation [70][71][72][73][74][75][76]. These two methods give consistent phase boundaries of the chiral phase transition in the chiral limit. Although SCL is the coarse lattice limit, we might expect that the observables related to the phase transition are not so sensitive to the coarse lattice spacing, since long-wave dynamics dominates the property of the phase transition, including the influences on higher-order cumulants [63]. Furthermore, one can take the chiral limit, in which the chiral transition becomes second order. Thus, the singular behavior associated with the phase transition is smeared by the finite-size effect only. As a result, the divergent part of the higher-order cumulants, which can be described by the relevant scaling function of the singular part of the free energy density [21], could be replaced by sign changes across the transition [21,22,60,61,66].
In this article, we focus on the fluctuations of net-baryon number and discuss the higher-order cumulant ratios in the strong-coupling limit of LQCD by utilizing the AFMC method to take the fluctuation effects into account. We here consider the LQCD action with one species of unrooted staggered fermion in the chiral limit. It should be noted that one species of unrooted staggered fermion at finite lattice spacing has O(2) symmetry, while it corresponds to four flavors in the continuum limit. We expect that we could extract qualitative behaviors of higher-order cumulants around TCP and CP according to the O(2) symmetry at finite lattice spacing, as discussed in Sects. 2.2 and 3.4, since the sign of the critical exponent α is the same for O (2) and O(4) symmetries. We demonstrate that the higher-order cumulant ratios, Sσ and κσ 2 , show oscillatory behavior, and that there exists a region where higher-order cumulants are negative in the strong-coupling and chiral limits on not-too-small lattices at μ/T 0.2. We also discuss the lattice-size dependence of the negative region.
This paper is organized as follows. We briefly provide the formalism in Sect. 2. Our main results on the higher-order fluctuations of the net-baryon number are presented in Sect. 3. We summarize our paper in Sect. 4. Some formulae for the net-baryon number cumulants used in Sect. 3 are found in the appendix.

Lattice QCD in the strong-coupling limit with fluctuations
We consider a lattice QCD action with one species of unrooted staggered fermion in the d(= 3) + 1-dimensional anisotropic Euclidean spacetime with N τ temporal and L spatial lattice sizes. This LQCD action has remnant chiral symmetry U (1) R × U (1) L in the chiral limit (m 0 → 0). We set the number of the color N c = 3 and the lattice unit a = 1 throughout this paper. In the following, we briefly summarize the formalism developed in Ref. [69].

Effective action in the strong-coupling limit with fluctuations
In the strong-coupling limit (SCL) g → ∞, we can ignore the plaquette terms, which are proportional to 1/g 2 , so the partition function and the action of the lattice QCD become where χ x and U ν,x represent the staggered quark field and the link variable, respectively, and V ± x and M x are mesonic composites. The staggered sign factor η j,x = (−1) x 0 +···+x j−1 is related to the gamma matrices in the continuum limit. The quark mass m 0 is taken to be zero throughout this paper. The quark chemical potential μ is introduced together with the physical lattice spacing ratio , where γ is the anisotropy parameter. We adopt f (γ ) = γ 2 following the arguments in Refs. [69,74,[87][88][89].
We obtain the effective action of the auxiliary fields, S AF eff , after the following three steps. First, by integrating out spatial link variables [77][78][79][80][81][82][83][84][85][86][87][88][89][90][95][96][97] in the leading order of the strong-coupling and 1/d (large-dimensional) expansions [83], one finds a convenient expression for the effective action: It should be noted that spatial baryon hopping terms are not included in Eq. (4) as they are higher order in the 1/d expansion. By comparison, we exactly integrate out temporal link variables later, then the temporal baryon hopping effects are taken into account.
We can now perform the Monte Carlo integral over the auxiliary fields (σ k,τ , π k,τ ) using the effective action S AF eff . We generate Monte Carlo configurations based on the phase-quenched action at finite μ and T and calculate observables by the reweighting method: where θ = −Im(S AF eff ) and · · · pq denotes the phase-quenched average. Through this AFMC method, we can numerically take account of fluctuation effects.
It should be noted that we have a sign problem that stems from the bosonization procedure, but it is not severe and we can investigate the QCD phase diagram as done in Ref. [69]. The reason for the milder sign problem may be understood as follows. First, the auxiliary field effective action is T. Ichihara et al. obtained by integrating out the link variables, then it contains only the color singlet field. Color singlet states are considered to be closer to energy eigenstates than colored states are, and we expect smaller phases in the path integral. Next, there is no sign problem in the mean-field approximation, where σ x is assumed to be constant and the π x fields are neglected; the sign problem is then not severe with small |k|. In the case where long-wave physics dominates, only small |k| auxiliary fields are relevant, and thus the π field can be regarded as almost constant. As seen in Eq. (8), the complex phase of one site has an opposite sign to that of the nearest-neighbor sites. Therefore, we could expect that the complex phase on one site coming from the bosonization is canceled out by the nearest-neighbor site contributions as long as we study the long-wave phenomena. As a result, we have a milder sign problem and can directly generate MC configurations at finite μ and T as long as the lattice size is not very large. In actual calculations, high |k| modes are not negligible, and the average phase factor exp(iθ) pq is suppressed; 2 exp(iθ) pq 0.9 and 0.4 for μ/T 0.8 on 4 3 × 4 and 8 3 × 8 lattices around the phase boundary, respectively [69]. The average phase factor is related to the free energy density as exp(iθ) pq = e − f ·L 3 /T , where the free energy density difference between full and phase-quenched simulations f = f − f pq characterizes the severity of the sign problem [74]. In our calculations, the free energy density difference takes almost the same values on 6 3 × 4, 6 4 , and 8 4 lattices. Supposing that ( f ) max 10 −3 as obtained at T 0.9 and μ/T = 0.8 on an 8 4 lattice [69], the average phase factor is evaluated to be about 0.329, 0.147, 0.0106 on 10 3 , 12 3 , and 16 3 spatial lattices. It is not impossible to perform reliable calculations with these average phase factors.
In Fig. 2, we show the chiral condensate (σ = τ σ k=0,τ /N τ ) and the logarithm of the baryon number susceptibility (χ ) in the chiral limit on a 6 3 × 6 lattice as a function of T /T c , where T c ( 1.468 62) denotes critical temperature at μ = 0 on a 6 3 × 6 lattice defined by the peak position of chiral susceptibility (χ σ = ∂ 2 log Z /∂m 2 0 /(L 3 N τ )) [69]. We can see the decrease of the chiral condensate with increasing T , indicating chiral phase transition. The decrease of the chiral condensate becomes steeper with increasing μ/T . The baryon number susceptibility has a sharp peak at high μ/T . These behaviors suggest the influence of the singularity at the tricritical point. It should be noted that the baryon number susceptibility at zero chemical potential decreases at high T . This trend is not consistent with standard lattice QCD results and could be an artifact of the truncation of the action in Eq. (4), so we only focus on the vicinity of the phase transition.

Some remarks on the universality class
The staggered fermion has a remnant chiral symmetry in the chiral limit and has O(2) symmetry as long as the lattice spacing is coarse. We introduce the auxiliary fields to maintain the original symmetry, O(2). The auxiliary field action S AF eff in Eq. (10) is invariant under the chiral transformation, which corresponds to the transformation of the staggered fermion fields: PTEP 2015, 113D01 T. Ichihara et al. In the framework of the MDP in the strong-coupling limit, the critical exponent was studied, the values of which are consistent with O(2) [71]. It should be noted that we need careful discussions about the difference in symmetry since the critical exponents are different in O (2), O(4), and Z(2). We will mention the relation between the O(N ) scaling function and critical behavior in Sect. 3.4.

Net-baryon number cumulants
In this section, we present results on the higher-order cumulant ratios of the net-baryon number.
The nth-order cumulant of the net-baryon number in the grand canonical ensemble is given by To discuss the effects of the phase transition on the net-baryon number fluctuations, it is convenient to take the ratio of a higher-order cumulant to the second-order one [102][103][104], since the volume factor, V , in Eq. (15) cancels. We show the results of the following cumulant ratios, the normalized skewness, and kurtosis: The skewness S and kurtosis κ probe the asymmetry with respect to the mean and peakedness of the underlying probability distribution, respectively. The normalization by χ (2) μ implies that one sets a reference, since κσ 2 = 1 for the Skellam distribution, which describes the distribution of the difference n 1 − n 2 , where n 1 and n 2 follow the Poisson distribution. Then it corresponds to a free hadron gas with the Boltzmann approximation.
The detailed method for the calculation of the cumulants is summarized in Appendix A. As discussed in Appendix A, the observables have an imaginary part due to the high |k| modes of π k,τ , so we take the real part of the observables and error bars in this article. The absolute mean values of the maximum imaginary parts for the normalized kurtosis and skewness are about 70 and 2 for μ/T = 0.8 and 0.5 and 0.04 for μ/T = 0.2, respectively, which are smaller than the peak heights and valley depths of the real part, as seen in Figs. 4 and 5. It should be noted that the high-temperature behavior of the cumulants should reveal the artifacts of the present treatment, presumably due to the absence of spatial baryon hopping effects. Figures 4 and 5 indicate that Sσ stays negative in the large chemical potential region and κσ 2 can also take small negative values depending on μ/T at high T .  For a free ideal baryon gas, we expect Sσ → 6μ/ 3μ 2 + π 2 and κσ 2 → 6/ 3μ 2 + π 2 at high temperature [37,102,103,105]. These differences at high temperature may be due to the truncation in the 1/d expansion. The spatial baryon hopping term is missing in the leading order in the 1/d expansion; the fermion momentum integral then generates only a volume factor. As a result, pressure is proportional to T rather than T 4 in the Stefan-Boltzmann behavior at high temperature. Since we are interested in the transition region, where we expect long-wave mesonic fluctuations to dominate, we focus on the cumulant ratios around the phase boundary in our subsequent discussion.

T and μ dependence
In Fig. 3, we show the results of Sσ and κσ 2 on a 6 4 lattice at several μ/T as a function of T /T c . One sees the strong dependence of these higher-order cumulants on both temperature and chemical potential. In the low chemical potential region, both Sσ and κσ 2 are positive and have a broad peak.
As μ/T increases, characteristic structures emerge. Sσ increases with T , then exhibits a strong positive peak followed by sudden decrease to a large negative value. Further increase of temperature leads to a constant small value. The negative region of the skewness, indicating the critical behavior [66], starts to appear between μ/T = 0.3 and μ/T = 0.5. κσ 2 shows a similarly moderate increase with temperature and a sharp positive peak. Then it becomes negative with a sharp valley but becomes positive again with a somewhat milder peak. Such an oscillatory behavior appears in μ/T 0.2. Both cumulant ratios have one negative valley, and the skewness (kurtosis) has one (two) positive peak(s) at high chemical potential. Comparing these values with the phase boundary [69], one finds these behaviors appear around the phase boundary (see below). In the following subsections, we discuss the behaviors in more detail by looking at lattice size dependence.

Lattice size dependence
In Figs  chiral condensate vanishes. We will use these mean-field values in Sect. 3.5 to remove the artifact, as discussed in Sect. 3.1.
The lattice size dependence of Sσ is moderate at low μ/T and prominent at large μ/T . From Fig. 4, one sees that the qualitative structure of the temperature dependence of Sσ does not strongly depend on the lattice size. For μ/T = 0.2, the skewness is positive around the phase-transition region and has a peak. While the peak height does not show a statistically significant dependence on the lattice size, one finds that the temperature dependence becomes slightly stronger in the largest one, 6 4 . The lattice size dependence becomes prominent at large chemical potential μ/T = 0.8. The skewness has one positive peak and one negative valley at μ/T = 0.8. The widths of the peak and the valley become narrower on a 6 4 lattice.
A similar lattice size dependence is found for κσ 2 . In the low chemical potential region (μ/T = 0.2), the smallest lattice size 4 4 does not exhibit the characteristic structure seen in larger lattice cases, but one sees a moderate decrease around T /T c 0.9 following a peak at lower T . The broad negative valley appears on the 6 3 × 4 and 6 4 lattices, implying that one should not use too small a lattice size to see the influence of the critical fluctuations. We find that the oscillatory behavior with two positive peaks and one negative valley and its lattice size dependence becomes prominent for μ/T = 0.8. Again, as lattice size increases, the peak height tends to have a larger value and the width of the oscillatory region becomes narrower.
For a nonzero μ, the first divergent cumulant is the third-order one; thus, Sσ diverges in the thermodynamic limit for the 3D O(2) and O(4) universality classes. Once the singular part is smeared by the explicit breaking (nonzero current quark mass) or finite-volume effects, 3 the leading singular contribution is suppressed by the small multiplicative factor −(2κ q )(μ/T ) n , where κ q denotes the curvature of the phase boundary near μ = 0 in the chiral limit [21]. Then, whether one sees the effects of critical fluctuations in the higher-order cumulants or not depends on the value of μ/T and the magnitude of the regular part of the free energy density and its derivatives. For instance, Sσ positively (negatively) diverges in the thermodynamic limit when approached from below (above) the phase transition. The negative divergent contribution is replaced by a large negative value in the presence of explicit breaking or finite-volume effects. Then, the appearance of the negative Sσ depends on whether the singular part overcomes the regular part. A similar argument applies to κσ 2 , and the absence of the negative region in κσ 2 in the smallest lattice (Fig. 5) is a consequence of the finitevolume effects. Since κσ 2 exhibits an oscillation around the phase transition along constant μ/T , including the narrow negative region between the two positive peaks, one expects κσ 2 to positively diverge both from below and above the phase transition and the negative region to disappear in the thermodynamic limit. The lattice size dependence of the negative κσ 2 region meets this expectation. This is in contrast to the case of the finite quark mass. With a finite quark mass, one finds the negative kurtosis region even in the thermodynamic limit as a remnant of the divergence in the chiral limit [21,60,61].

Kurtosis in the QCD phase diagram in the strong-coupling limit
In Fig. 6, we show the negative kurtosis region in the chiral limit (m 0 → 0) on a 6 4 lattice. We also show the phase boundary determined by the peak position of the chiral susceptibility (χ σ = ∂ 2 log Z /∂m 2 0 /(L 3 N τ )) for μ/T ≤ 0.8 and by the cross point of the effective potential defined by F eff = S eff /(N τ L 3 ) for μ/T ≥ 1.0. The transition would be of second order for μ/T ≤ 0.8. While the numerical finite-size scaling analysis cannot rule out crossover due to the statistics [69], we expect a second-order transition from the symmetry arguments [23]. We guess that the tricritical point exists at μ/T > 0.8.
We here define the negative kurtosis region, where κσ 2 is smaller than the mean-field results at high T , in order to reduce the artifact at high T , as discussed in Sect. 3.1. It should be noted that 3 Recalling that the behavior of the order parameters becomes moderate not only by finite-mass but also finite-size effects, we could guess that the singular part of the cumulants is also masked by the finite-size effects. By using models, we can find that the characteristic behavior of physical quantities becomes smoother when the system is finite [114][115][116]. In the framework of the 3D Ising model, finite-size results are studied in Ref. [117] and the peak heights of Sσ and κσ 2 increase with larger volume due to the correlation length ξ . The relations between higher-order net-baryon number cumulants and the correlation length in the 3D O(4) spin model are pointed out in Ref. [118]. we do not take account of error bars to define the region here. The negative kurtosis region appears from about μ/T = 0.2, and the region is the largest at μ/T = 0.3. 4 By comparison, the kurtosis is positive at μ/T < 0.2. This might be due to the suppression of the singular contribution by the factor (μ/T ) n and the regular part dominates over the singular part [21].
The negative κσ 2 region almost coincides with the phase boundary. As discussed in Sects. 3.3 and 3.4, the region will shrink to the phase boundary and finally vanish in the thermodynamic limit. Although the dependence of the scaling function is different between the finite-size effect and the symmetry-breaking term, one could expect that the negative region also appears with finite mass in the strong-coupling limit.
We conclude that we find negative skewness or kurtosis regions in the AFMC method in the chiral and strong-coupling limits due to the finite-volume effects as for the kurtosis, as a consequence of the critical fluctuations around the phase boundary incorporated through the AFMC method.

Summary
We have investigated the higher-order cumulant ratios of the net-baryon number, Sσ = χ (3) μ /χ (2) μ , κσ 2 = χ (4) μ /χ (2) μ , in the strong-coupling (g → ∞) and chiral limits (m 0 → 0) of QCD in the leading order of large-dimensional expansions. Mesonic fluctuation effects are taken into account by making use of the auxiliary field Monte Carlo (AFMC) method. We find negative skewness and kurtosis regions around the boundary of the chiral phase transition. The skewness and kurtosis exhibit characteristic temperature dependences influenced by the critical fluctuations. The skewness is found to be negative near the phase boundary for large μ/T . The oscillatory behavior around the phase boundary seems to be consistent with the expectations from the finitevolume effect such that the positive (negative) peak at lower (higher) temperature diverges in the thermodynamic limit [21]. Similarly, the kurtosis has one negative valley between two positive peaks for large μ/T . With increasing lattice size, the negative valley is found to shrink, as anticipated for the finite-volume effect. Thus, we expect the two positive peaks to diverge and the negative region to disappear in the thermodynamic limit [21]. It should be noted that we need careful treatments to discuss the connection between strong-coupling lattice QCD studies on finite-size lattices and experiments. In Ref. [113], the relations between finite-size effects in standard LQCD calculations and momentum cut-off effects in experiments are addressed.
One of the important next steps could be to study the effect of finite mass. This can be carried out by applying the AFMC method. Another important step is to see the finite-size effects [27,28,[114][115][116][119][120][121][122] and determine the negative kurtosis area with or without the bare quark mass. Whether one could have such a region in lattice QCD is a mandatory question. As discussed above, it will depend on whether the smeared singular contribution by finite-volume and/or explicit symmetry-breaking terms can overwhelm the regular contribution [21,[114][115][116], which corresponds to a hadron resonance gas in finite-temperature QCD below the phase transition. Since our present formulation ignores the spatial baryon hopping, we would expect that the regular contribution might be smaller than the realistic case; thus, our results could serve as a lower limit for the value of the chemical potential where κσ 2 < 0.

Appendix A. Higher-order derivatives of the baryon chemical potential
In this appendix, we show the higher-order derivatives with respect to the dimensionless chemical potentialμ(= N c μ/T ). In the following, the action S and the partition function Z are denoted as S AF eff and Z AF , respectively. Then, higher-order derivatives are given as whereμ = N c μ/T . The derivatives of the action with respect toμ are given as where S = k,τ, f >0 In this analysis, we have a complex phase coming from the fermion determinant, so the observables also take complex values even if the observables are real in essence. Invoking the cancellation of the imaginary part in the case of higher statistics in principle, we can take the real part of the observables. We here use the jackknife method in order to evaluate the statistical errors (see, e.g., Ref. [123]). Since the observables have an imaginary part, we evaluate the error with complex phase and take the correlation between complex observables and the complex phase into account. Then, we take the real part of the obtained error bars.