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

###### Abstract

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, and , show oscillatory behavior around the phase boundary at , and there exists the region where the higher-order cumulant ratios are negative. The negative region of is found to shrink with increasing lattice size. This behavior agrees with the expectations from the scaling analysis.

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

Terukazu Ichihara, Kenji Morita and Akira Ohnishi

Department of Physics, Kyoto University, Kyoto 606-8502, Japan

Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

Institute of Theoretical Physics, University of Wroclaw, PL-50204 Wrocław, Poland

E-mail: t-ichi@ruby.scphys.kyoto-u.ac.jp, †ohnishi@yukawa.kyoto-u.ac.jp, ‡kmorita@yukawa.kyoto-u.ac.jp

\@hangfrom Subject Index B64, D28, D30, D34

## 1 Introduction

Fluctuations of conserved charges are promising observables in
relativistic heavy-ion collisions in search for QCD phase transition
SRS1999 [], AHM-JK2000 [], Koch_review [].
In the first phase of the Beam Energy Scan (BES I) program at Relativistic
Heavy Ion Collider (RHIC), net-proton BESI-proton [],
as a proxy of net-baryon Hatta-Stephanov [], KA_proton-baryon [],
and net-electric charge event-by-event fluctuations BES1-charge []
have been measured in Au+Au collisions for a broad energy range from
to 200 GeV.
The success of the statistical thermal model Thermalmodel []
indicates 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}^{1}1This argument
holds for a limited acceptance. If one covers entire phase space in heavy
ion collisions, those fluctuations is of course absent. One can invoke
more informations by changing the rapidity acceptance beyond the equilibrium
regime Kitazawa_PLB728 [].

The particle number fluctuations of the conserved charge, i.e., susceptibilities or cumulants, reflect the property of the phase in QCD SFR_susceptibility [], which is expected to exhibit a rich structure Fukushima-Hatsuda []. Especially, if QCD has a critical point (CP) QCDCP [], the susceptibilities and higher-order cumulants diverge at the critical point owing to the divergent correlation length Hatta-Ikeda [], Stephanov_higherorder [], thus one expects anomalously large fluctuations would be observed if the system passes around CP. In heavy-ion collisions, the correlation length cannot diverge due to finite size of the system, but these fluctuations are expected to remain sensitive to the remnant criticality of the system around CP SRS1999 []. Since this sensitivity implies the tail of the event-by-event multiplicity distribution has importance in the behavior of the cumulants Morita_prob [], the analysis becomes statistically demanded.

The measured higher-order cumulant ratios of net-proton number, and , show non-monotonic behavior as a function of the incident energy, where , and are referred to as the variance, skewness and kurtosis, respectively. At around , the data show below the expected value in the Skellam distribution BESI-proton [].

The decrease of in the net-proton number can be the signal of the critical behavior. According to the theoretical arguments SF-4-CM [], Stephanov_negativekappa [], 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 two-flavor massless quarks with axial anomaly, the QCD phase transition at finite () belongs to the universality class of 3d O(4) symmetric spin model PW84 [], MagneticEOS [], and the second order transition at low may turn into the first order transition at the tricritical point (TCP). At physical quark masses, finite phase transition would be crossover FSS [], 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 Z(2) universality class QCDCP [], Hatta-Ikeda [], Fujii-Ohtani [].

While the universality argument dictates the singular behavior of the thermodynamic quantities close to the phase transition, the actual magnitude of the fluctuations are smeared by the finite volume effect in addition to the finite quark mass. Finite volume makes the transition smoother. We need to perform finite size scaling analysis to observe precise critical behavior FSS []. 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 SF-4-CM []. 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 : Z(2) critical behavior around CP Stephanov_higherorder [], the remnant of the O(4) chiral phase transition SF-4-CM [], 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 HM-Lattice-zero-chemi-pote [], Jin:2014hea [], Allton:2005gk [] and at finite chemical potential HM-lattice [], Jin:2013wta [], Allton:2005gk [], Takeda:2014nta [], Gattringer:2014hra []. At large baryon chemical potential, however, LQCD faces the notorious sign problem which makes it difficult to carry out MC simulations owing to complex fermion determinant. There are many attempts to circumvent the sign problem Reweighting [], Taylor [], Histogram [], ImagMu [], Langevin [], LThimble [], Fugacity [], Canonical []. Most of the methods evading the sign problem are reliable only in limited circumstances such as , small volume, or heavier quark masses than physical one, depending on the method. Conclusive results at physical point have not been obtained yet.

Owing to the sign problem in LQCD, most of fluctuation studies at nonzero density have been carried out by using chiral models CModel-3rd [], CModel-HM [], FRG-HM [], SF-4-CM [], Morita_prob [], Stokic [], Morita_ptep []. It has been pointed out that taking into account fluctuations around the phase transition is essential to correctly describe the behavior of the fluctuation of conserved charges FRG-HM [], SF-4-CM [], Morita_prob []. In chiral models, this has been done by implementing functional renormalization group (FRG) FRG [], Schaefer-Wambach []. This is a natural consequence of the fact that the critical behavior of the conserved charge fluctuations is governed by the critical exponent of the specific heat which requires beyond mean-field treatment FRG-HM [], SF-4-CM [], Morita_prob [], Stokic [], Morita_ptep [].

One of the possible alternative ways to attack the finite density region is the strong coupling approach of LQCD with fermions SCL [], large_d [], Faldt [], BilicDeme [], MF-SCL [], Bilic [], Jolicoeur [], NLOchiral [], Aoki:1987us [], NLOPD [], NNLO [], KarschMutter [], MDP [], MDP-CE [], SC-Rewei [], SC-Pol [], Ichihara:2014ova [], PhDFromm []. 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: auxiliary field Monte-Carlo (AFMC) method Ichihara:2014ova [] and the monomer-dimer-polymer (MDP) simulation MDP [], MDP-CE [], SC-Rewei [], KarschMutter [], PhDFromm []. 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 is not so sensitive to the coarse lattice spacing since long-wave dynamics dominates the property of the phase transition, including the influences on the higher-order cumulants Morita_ptep []. 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 relevant scaling function of the singular part of the free energy density SF-4-CM [], could be replaced by sign changes across the transition CModel-3rd [], FRG-HM [], SF-4-CM [], Stephanov_negativekappa [].

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. We demonstrate that the higher-order cumulant ratios, and , show oscillatory behavior, and there exists the region where higher-order cumulants are negative in the strong coupling and chiral limit on not too small lattices at . We also discuss the lattice size dependence of the negative region.

This paper is organized as follows. We briefly provide the formalism in Sec. 2. Our main results on the higher-order fluctuations of the net-baryon number will be presented in Sec. 3. We summarize our paper in Sec. 4. Some formulae for the net-baryon number cumulants used in Sec. 3 are found in the appendix.

## 2 Lattice QCD in the strong coupling limit with fluctuations

We consider a lattice QCD action with one species of unrooted staggered fermion in the dimensional anisotropic Euclidean spacetime with temporal and spatial lattice sizes. This LQCD action has remnant chiral symmetry in the chiral limit (). We set the number of the color and the lattice unit throughout this paper. In the following, we briefly summarize the formalism developed in Ref. Ichihara:2014ova [].

### 2.1 Effective action in the strong coupling limit with fluctuations

In the strong coupling limit (SCL) , we could ignore the plaquette terms, which are proportional to , so the partition function and the action of the lattice QCD become

(1) | ||||

(2) | ||||

(3) |

where and represent the staggered quark field and the link variable, respectively, and and are mesonic composites. The staggered sign factor is related to the gamma matrices in the continuum limit. The quark mass is taken to be zero throughout this paper. Quark chemical potential is introduced together with the physical lattice spacing ratio , where is the anisotropy parameter. We adopt following the arguments in Refs. Bilic [], BilicDeme [], PhDFromm [], Ichihara:2014ova []

We obtain the effective action of the auxiliary fields, , after the following three steps. First, by integrating out spatial link variables SCL [], large_d [], MF-SCL [], Faldt [], BilicDeme [], Bilic [], NLOPD [], NNLO [], Jolicoeur [] in the leading order of the strong coupling and (large-dimensional) expansion large_d [], one finds a convenient expression for the effective action,

(4) |

It should be noted that spatial baryonic hopping terms are not included in Eq. (2.1) as they are higher order in the expansion. By comparison, we exactly integrate out temporal link variables later, then the temporal baryonic hopping effects are taken into account.

In the second step, we transform the effective action to the fermion-bilinear form by using the extended Hubbard-Stratonovich (eHS) transformation NLOPD [], NNLO [],

(5) |

where and . The four-Fermi interaction terms in Eq. (2.1) are diagonal in the momentum space, and are separated into two parts based on the momentum regions: the positive modes () and the negative () modes.

(6) |

where , , and . The effective action after the eHS transformation of Eq. (6) reads

(7) | ||||

(8) |

where and . Note that and , since terms in Eq. (6) are bosonized at a time. At small , they are also regarded as the chiral () fields and the Nambu-Goldstone () fields, respectively, since the staggered fermion identifies the spin and flavor by specifying space-time position. The sign factor, , corresponds to in the spinor-taste space in the continuum limit.

In the third step, by integrating over the Grassmann and temporal link () variables Faldt [], BilicDeme [], Bilic [], the partition function and the effective action are reduced to

(9) | ||||

(10) | ||||

(11) |

where and . In the last line, we use a recursion formula to obtain Faldt [], Bilic [], BilicDeme []. In the cases where is static, we find . It should be noted that the action represents the confinement phase, since the baryonic chemical potential appears as , which can be understood as baryonic contribution SC-Pol [], similarly to the PNJL model with vanishing Polyakov loop Pol [].

We can now perform the Monte-Carlo integral over the auxiliary fields using the effective action . We generate Monte-Carlo configurations based on the phase quenched action at finite and and calculate observables by the reweighting method;

(12) |

where and 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. Ichihara:2014ova [].
The reason for the milder sign problem may be understood as follows.
First, the auxiliary field effective action is obtained by
integrating out the link variables,
then it contains only the color singlet field.
Color singlet states are, in general, closer to the energy eigenstates
than colored states, and we expect smaller phases in the path integral.
Next, there is no sign problem in the mean field approximation,
where is assumed to be constant and fields are neglected,
then the sign problem is not severe with small .
In the case where
the long wave physics dominates,
only small auxiliary fields are relevant,
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
as long as the lattice size is not very large.
In actual calculations, high modes are not
negligible, and the average phase factor
is suppressed;
^{2}^{2}2The average phase factor indicates the severity of weight cancellation.
We have no weight cancellation when .
and 0.4
for
on the and lattices around the phase boundary, respectively Ichihara:2014ova [].

In Fig. 1, we show the chiral condensate () and the logarithm of the baryon number susceptibility () in the chiral limit on a lattice as a function of , where denotes critical temperature at on a lattice defined by the peak position of chiral susceptibility () Ichihara:2014ova []. We can see the chiral phase transition behavior from the decrease of the chiral condensate with increasing . The decrease of the chiral condensate becomes steeper with increasing . The baryon number susceptibility has a sharp peak at high . 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 . This trend is not consistent with standard lattice QCD results and would be an artifact of the truncation of the action in Eq. (2.1), so we only focus on the vicinity of the phase transition.

### 2.2 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 in Eq. (10) is invariant under the chiral transformation,

(13) |

which corresponds to the transformation of the staggered fermion fields

(14) |

In the framework of the MDP in the strong coupling limit, the critical exponent was studied and that values are consistent with O(2) MDP-CE []. It should be noted that we need careful discussions about the difference in symmetry since the critical exponents are different in , , and Z(2). We will mention the relation between the O() scaling function and critical behavior in Sec. 3.4.

## 3 Net-baryon number fluctuations in the strong coupling limit

### 3.1 Net-baryon number cumulants

In this section, we present results on the higher-order cumulant ratios of the net-baryon number. The -th order cumulant of the net-baryon number in the grand canonical ensemble is given by

(15) |

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 CumuRatio [], CumuRatio-1 [] since the volume factor, , in Eq. (15) cancels. We show the results of the following cumulant ratios, the normalized skewness and kurtosis,

(16) |

The skewness and kurtosis probe the asymmetry with respect to the mean and peakedness of the underlying probability distribution, respectively. The normalization by implies that one sets a reference, since for the Skellam distribution which describes the distribution of the difference , where and follow the Poisson distribution. Then it corresponds to the free hadron gas with Boltzmann approximation.

The detailed method for the calculation of the cumulants are summarized in Sec. A. As discussed in Sec. A, observables have an imaginary part due to the high modes of , so we take the real part of the observables and error bars in this article. The absolute mean value of maximum imaginary part for the normalized kurtosis and skewness are about 70 and 2 for and 0.5 and 0.04 for , respectively, which are smaller than the peak heights and valley depth of the real part as seen in Figs. 3 and 4. It should be noted that the high temperature behavior of the cumulants would bare the artifacts of the present treatment. Figures 3 and 4 indicate that stays negative in the large chemical potential region and can also take small negative values depending on at high . For a free ideal baryon gas we expect and at high temperature CumuRatio [], Stokic:2008jh [], Allton:2005gk []. These differences at high temperature may be due to the truncation in the expansion. Spatial baryon hopping term is missing in the leading order in the expansion, then the fermion momentum integral generates only a volume factor. As a result, pressure is proportional to rather than in the Stefan-Boltzmann behavior at high temperature. Since we are interested in the transition region, where we expect that the long wave mesonic fluctuations dominate, we focus on the cumulant ratios around the phase boundary in the later discussion.

### 3.2 and dependence

In Fig. 2, we show the results of and on a lattice at several as a function of . One sees the strong dependence of these higher-order cumulants on both temperature and chemical potential. In the low chemical potential region, both and are positive and have a broad peak.

As increases, characteristic structures emerge. increases with , then exhibits a strong positive peak followed by sudden decrease to large negative value. Further increase of temperature leads to a constant small value. The negative region of the skewness, indicating the critical behavior CModel-3rd [], starts to appear between and . shows similarly moderate increase with temperature and a sharp positive peak. Then it turns to be negative with a sharp valley but becomes positive again with a somewhat milder peak. Such an oscillatory behavior appears in . Both cumulant ratios have one negative valley, and the skewness (kurtosis) has one (two) positive peak(s) at high chemical potential. Comparing with the phase boundary Ichihara:2014ova [], one finds these behaviors appear around the phase boundary (see below). In the following subsections, we discuss the behaviors in more details by looking at lattice size dependence.

### 3.3 Lattice size dependence

In Figs. 3 and 4, we show the results of the normalized skewness and the normalized kurtosis , respectively, on the fixed lines, (left panel) and (right panel), on and lattices. Horizontal lines show the mean field values at high temperature, where the chiral condensate vanishes. We will use these mean field values in Sec. 3.5 to remove the artifact as discussed in Sec. 3.1

Lattice size dependence of is moderate at low and prominent at large . From Fig. 3, one sees that the qualitative structure of the temperature dependence of does not strongly depend on the lattice size. For , the skewness is positive around the phase transition region and has a peak. While the peak height does not show statistically significant dependence on the lattice size, one finds the temperature dependence becomes slightly stronger in the largest one, . The lattice size dependence becomes prominent at large chemical potential . The skewness has one positive peak and one negative valley at . The widths of the peak and the valley become narrower on a lattice.

Similar lattice size dependence is found for . In the low chemical potential region (), the smallest lattice size does not exhibit the characteristic structure seen in larger lattice cases, but one sees a moderate decrease around following a peak at lower . The broad negative valley appears on and lattices, implying that one should not take too small 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 . Again, as lattice size becomes larger, the peak height tends to be larger value and the width of the oscillatory region becomes narrower.

### 3.4 Relation to O() scaling

The behavior of the higher-order net-baryon number cumulants around chiral phase transition are expected to be described by the property of the O() scaling function O2-CE-alpha [], O4-CE-alpha [], SF-4 [], O2-O4-study [], SF-4-CM []. According to Ref. SF-4-CM [], the first divergent cumulant at in the thermodynamic limit is the sixth order one for and . This result comes from the fact that the critical exponent of the specific heat is negative in these cases, in O(2) O2-CE-alpha [] and in O(4) O4-CE-alpha []. For positive , the divergence appears at the fourth order cumulants while it has a cusp for FRG-HM [], Morita_prob [].

For a nonzero , the first divergent cumulant is the third order
one thus diverges in the thermodynamic limit for 3d O(2) and O(4)
universality class. Once the singular part is smeared by the explicit
breaking (nonzero current quark mass) or finite volume effects
^{3}^{3}3
Recalling 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 the characteristic behavior of physical quantities becomes smoother when the system is finite ModelStudy-Volume-effect [].
In the framework of 3d Ising model, finite size results are studied in Ref. Chen:2010ej []
and the peak heights of and increase with larger volume
due to the correlation length . The relations between higher-order net-baryon number cumulants and the correlation length in 3d O(4) spin model are pointed out in Ref. chiral-model-and-scaling [].
,
the leading singular contribution is suppressed by
the small multiplicative
factor where denotes the curvature of
the phase boundary near in the chiral limit SF-4-CM [].
Then, whether one sees the effects of critical fluctuations in the
higher-order cumulants or not depends on the value of and
magnitude of the regular part of the free energy density and its derivatives.
For instance, 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 the explicit breaking or finite volume effects.
Then, the appearance of the negative depends on whether the
singular part overcomes the regular part.
Similar argument applies to , and the
absence of the negative
region in
in the smallest lattice (Fig. 4)
is a consequence of the finite volume effects.
Since
exhibits the oscillation
around the phase transition along constant including the narrow negative region between the two positive peaks,
one expects
positively diverges both from below and above the phase transition
and negative region would disappear in the thermodynamic limit.
The lattice size dependence of the negative region
meets this expectation.
This is in contrast to the case of the finite quark mass.
With finite quark mass,
one finds the negative kurtosis region even in the thermodynamic limit as a remnant of the divergence
in the chiral limit SF-4-CM [], FRG-HM [].

### 3.5 Kurtosis in the QCD phase diagram in the strong coupling limit

In Fig. 5, we show the negative kurtosis region in the chiral limit on a lattice. We also show the phase boundary determined by the peak position of the chiral susceptibility () for and by the cross point of the effective potential defined by for . The transition would be the second order for . While the numerical finite size scaling analysis cannot rule out crossover due to the statistics Ichihara:2014ova [], we expect the second order transition from symmetry arguments PW84 []. We guess that the tricritical point exists at .

We here define the negative kurtosis region where is smaller
than the mean field results at high ,
in order to reduce the artifact at high as discussed in Sec. 3.1.
It should be noted that we do not take account of error bars to define the region here.
The negative kurtosis region appears from about ,
and the region is the largest at
^{4}^{4}4
Along the line at high ,
the mean field value and AFMC results are so close to each other
and
it seems that the AFMC results oscillate around the MF result.
Thus we define the boundary of the negative kurtosis region as the first intersection of the AFMC and MF results.
.
By comparison, the kurtosis is positive at .
This might be due to
the suppression of the singular contribution by the factor and the regular part dominates
over the singular part SF-4-CM [].

The negative region almost coincides with the phase boundary. As discussed in Secs. 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 the negative skewness or kurtosis region in the AFMC method in the chiral and strong coupling limit 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.

## 4 Summary

We have investigated the higher-order cumulant ratios of the net-baryon number, , in the strong coupling and the chiral limit () 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 the negative skewness and the kurtosis region 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 . The oscillatory behavior around the phase boundary seems to be consistent with the expectations from the finite volume effect such that the positive (negative) peak at lower (higher) temperature diverges in the thermodynamic limit SF-4-CM []. Similarly, the kurtosis has one negative valley between two positive peaks for large . With increasing lattice size, the negative valley is found to shrink, as anticipated for the finite volume effect. Thus, we expect two positive peaks diverge and the negative region disappears in the thermodynamic limit SF-4-CM [].

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 FSS [], ModelStudy-Volume-effect [], Vol-effect-study [] and determine the negative kurtosis area with or without bare quark mass. Whether one could have such 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 term can overwhelm the regular contribution SF-4-CM [], ModelStudy-Volume-effect [], which corresponds to 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 chemical potential where .

## Acknowledgment

The authors would like to thank Sinya Aoki, Frithjof Karsch, Swagato Mukherjee, Hiroshi Ohno, Philippe de Forcrand, Hideaki Iida, Yu Maezawa, Keitaro Nagata, Shuntaro Sakai, Takahiro Sasaki, and Wolfgang Unger for useful discussions. Part of numerical computation in this work was carried out at the Yukawa Institute Computer Facility. TI is supported by the Grants-in-Aid for JSPS Fellows (No.25-2059) and would like to thank for RIKEN-BNL Brain Circulation Program. KM was supported by Polish Science Foundation (NCN), under Maestro grant 2013/10/A/ST2/00106. This work is supported in part by the Grants-in-Aid for Scientific Research from JSPS (Nos. 23340067, 24340054, 24540271, 15K05079), by the Grants-in-Aid for Scientific Research on Innovative Areas from MEXT (No. 2404: 24105001, 24105008), by the Yukawa International Program for Quark-hadron Sciences.

## A Higher order derivatives of the baryon chemical potential

In this appendix, we show the higher order derivatives with respect to dimensionless chemical potential . In the following, the action and the partition function are denoted as and , respectively. Then, higher order derivatives are given as

(A1) | ||||

(A2) | ||||

(A3) | ||||

(A4) |

where . The derivatives of the action with respect to are given as

(A5) | ||||

(A6) | ||||

(A7) | ||||

(A8) |

where and .

In this analysis, we have a complex phase coming from the fermion determinant, so the observables also take a complex value 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 JK []. 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.

When we evaluate the fourth derivative of with respect to , , (the function form is like ), we generate jackknife samples of each expectation value (). Then, we calculate the jackknife samples of the observable (). Finally we calculate expectation value and error bars of the observable by using . When we calculate the cumulant ratios (), we use the same technique as well.

## References

- [1] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998); Phys. Rev. D 60, 114028 (1999).
- [2] M. Asakawa, U. W. Heinz and B. Müller, Phys. Rev. Lett. 85, 2072 (2000); S. Jeon and V. Koch, Phys. Rev. Lett. 85, 2076 (2000).
- [3] V. Koch, arXiv: 0810.2520.
- [4] M. M. Aggarwal et al. (STAR Collaboration), Phys. Rev. Lett. 105, 022302 (2010); L. Adamczyk, et al. (STAR Collaboration), Phys. Rev. Lett. 112, 032302 (2014).
- [5] Y. Hatta and M. A. Stephanov, Phys. Rev. Lett. 91, 102003 (2003).
- [6] M. Kitazawa and M. Asakawa, Phys. Rev. C 85, 021901 (2012).
- [7] L. Adamczyk et al. (STAR Collaboration), Phys. Rev. Lett. 113, 092301 (2014).
- [8] P. Braun-Munzinger, K. Redlich, and J. Stachel, in Quark-GluonPlasma 3, edited by R. C. Hwa and X. N. Wang (World Scientific, 2004) p. 491.
- [9] M. Kitazawa, M. Asakawa, and H. Ono, Phys. Lett. B 728, 386 (2014).
- [10] C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. D 75, 074013 (2007).
- [11] K. Fukushima and T. Hatsuda, Rep. Prog. Phys. 74, 014001 (2011).
- [12] M. Asakawa and K. Yazaki, Nucl. Phys. A 504, 668 (1989).
- [13] Y. Hatta and T. Ikeda, Phys. Rev. D 67, 014028 (2003).
- [14] M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
- [15] K. Morita, B. Friman, K. Redlich, and V. Skokov, Phys. Rev. C 88, 034903 (2013); K. Morita, V. Skokov, B. Friman, and K. Redlich, Eur. Phys. J. C 74, 2706 (2014); K. Morita, B. Friman, and K. Redlich, Phys. Lett. B 741, 178 (2015).
- [16] B. Friman, F. Karsch, K. Redlich, and V. Skokov, Eur. Phys. J. C 71, 1694 (2011).
- [17] M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
- [18] R. D. Pisarski and F. Wilczek, Phys. Rev. D 29, 338 (1984).
- [19] S. Ejiri, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner, W. Unger, Phys. Rev. D 80, 094505 (2009); O. Kaczmarek, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner, and W. Unger, Phys. Rev. D 83, 014504 (2011); A. Bazavov, et. al., (HotQCD Collaboration), Phys. Rev. D 85, 054503 (2012).
- [20] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature 443, 675 (2006); G. Endrodi, Z. Fodor, S. D. Katz, K. K. Szabo, JHEP 1104, 001 (2011).
- [21] H. Fujii and M. Ohtani, Phys. Rev. D 70, 014016 (2004).
- [22] A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa and S. Mukherjee et al., Phys. Rev. Lett. 111, 082301 (2013) [arXiv:1304.7220 [hep-lat]]; R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz and C. Ratti, Phys. Rev. Lett. 111, 202302 (2013) [arXiv:1305.6297 [hep-lat]]; A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa and S. Mukherjee et al., Phys. Lett. B 737, 210 (2014) [arXiv:1404.4043 [hep-lat]]; A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa and S. Mukherjee et al., Phys. Rev. Lett. 113, 072001 (2014) [arXiv:1404.6511 [hep-lat]]; S. Sharma [Bielefeld-BNL-CCNU Collaboration], Nucl. Phys. A 931, 1199 (2014) [arXiv:1408.2332 [hep-lat]]; H. T. Ding, Nucl. Phys. A 931, 52 (2014) [arXiv:1408.5236 [hep-lat]].
- [23] X. Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda and A. Ukawa, Phys. Rev. D 91, 014508 (2015) [arXiv:1411.7461 [hep-lat]].
- [24] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, Phys. Rev. D 71, 054508 (2005) [hep-lat/0501030];
- [25] C. Schmidt, Prog. Theor. Phys. Suppl. 186, 563 (2010) [arXiv:1007.5164 [hep-lat]]; R. V. Gavai and S. Gupta, Phys. Lett. B 696, 459 (2011) [arXiv:1001.3796 [hep-lat]]; A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, S. Mukherjee and P. Petreczky et al., Phys. Rev. Lett. 109, 192302 (2012) [arXiv:1208.1220 [hep-lat]].
- [26] S. Takeda, X. Y. Jin, Y. Kuramashi, Y. Nakamura and A. Ukawa, PoS LATTICE 2014, 196 (2014) [arXiv:1411.1148 [hep-lat]].
- [27] X. Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda and A. Ukawa, Phys. Rev. D 88, 094508 (2013) [arXiv:1307.7205 [hep-lat]].
- [28] C. Gattringer and H. P. Schadler, Phys. Rev. D 91, 074511 (2015) [arXiv:1411.5133 [hep-lat]].
- [29] Z. Fodor, S. D. Katz, JHEP 0203, 014 (2002); JHEP 0404 , 050 (2004).
- [30] C. R. Allton, M. Döring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, and K. Redlich, Phys. Rev. D 71, 054508 (2005); R. V. Gavai and S. Gupta, Phys. Rev. D 71, 114014 (2005) [hep-lat/0412035]; O. Kaczmarek, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt, W. Soeldner, and W. Unger, Phys. Rev. D 83, 014504 (2011).
- [31] H. Saito et al. [WHOT-QCD Collaboration], Phys. Rev. D 84, 054502 (2011) [Erratum-ibid. D 85, 079902 (2012)].
- [32] P. de Forcrand and O. Philipsen, Nucl. Phys. B 642, 290 (2002); JHEP 0701, 077 (2007); M. D’Elia and M. -P. Lombardo, Phys. Rev. D 67, 014505 (2003).
- [33] G. Aarts, E. Seiler and I. -O. Stamatescu, Phys. Rev. D 81, 054508 (2010); E. Seiler, D. Sexty, I. -O. Stamatescu, Phys. Lett. B 723, 213 (2013).
- [34] M. Cristoforetti, F. D. Renzo, and L. Scorzato [AuroraScience Collaboration], Phys. Rev. D 86, 074506 (2012); H. Fujii, D. Honda, M. Kato, Y. Kikukawa, S. Komatsu and T. Sano, JHEP 1310, 147 (2013).
- [35] K. Nagata and A. Nakamura, Phys. Rev. D 82, 094027 (2010).
- [36] A. Hasenfratz and D. Toussaint, Nucl. Phys. B 371, 539 (1992); S. Ejiri, Phys. Rev. D 78, 074507 (2008); A. Li, A. Alexandru and K. -F. Liu, Phys. Rev. D 84, 071503 (2011); K. Nagata, S. Motoki, Y. Nakagawa, A. Nakamura, and T. Saito (XQCD-J Collaboration), Prog. Theor. Exp. Phys. 2012, 01A103 (2012).
- [37] V. Skokov, B. Stokić, B. Friman, and K. Redlich, Phys. Rev. C 82, 015206 (2010); V. Skokov, B. Friman, and K. Redlich, Phys. Rev. C 83, 054904 (2011).
- [38] B. Stokić, B. Friman, and K. Redlich, Eur. Phys. J. C 67, 425 (2010).
- [39] K. Morita and K. Redlich, Prog. Theor. Exp. Phys. 2015, 043D03 (2015).
- [40] B. J. Schaefer and M. Wagner, Phys. Rev. D 85, 034027 (2012) [arXiv:1111.6871 [hep-ph]]; A. Bhattacharyya, S. K. Ghosh, A. Lahiri, S. Majumder, S. Raha and R. Ray, Phys. Rev. C 89, 064905 (2014) [arXiv:1212.6134 [hep-ph]].
- [41] M. Asakawa, S. Ejiri, and M. Kitazawa, Phys. Rev. Lett. 103, 262301 (2009),
- [42] B. -J. Schaefer and J. Wambach, Nucl. Phys. A 757, 479 (2005).
- [43] C. Wetterich, Phys. Lett. B 301, 90 (1993).
- [44] T. Ichihara, A. Ohnishi, and T. Z. Nakano, PTEP 2014 12, 123D02 (2014).
- [45] F. Karsch and K. H. Mutter, Nucl. Phys. B 313, 541 (1989).
- [46] G. Boyd, J. Fingberg, F. Karsch, L. Kärkkäinen, and B. Petersson, Nucl. Phys. B 376, 199 (1992).
- [47] P. de Forcrand and M. Fromm, Phys. Rev. Lett. 104, 112005 (2010) [arXiv:0907.1915 [hep-lat]]; W. Unger and P. de Forcrand, J. Phys. G 38, 124190 (2011) [arXiv:1107.1553 [hep-lat]].
- [48] M. Fromm, Ph.D. thesis, ETH-19297, Eidgenssische Technische Hochschule ETH Zrich, 2010.
- [49] P. de Forcrand, J. Langelage, O. Philipsen, and W. Unger, PoS LATTICE 2013, 142 (2013) [arXiv:1312.0589 [hep-lat]]; W. Unger, Acta Phys. Polon. Supp. 7 No. 1, 127 (2014).
- [50] N. Kawamoto and J. Smit, Nucl. Phys. B 192, 100 (1981); H. Kluberg-Stern, A. Morel, B. Petersson, Phys. Lett. B 114, 152 (1982); J. Hoek, N. Kawamoto, and J. Smit, Nucl. Phys. B 199, 495 (1982); P. H. Damgaard, N. Kawamoto and K. Shigemoto, Phys. Rev. Lett. 53, 2211 (1984); P. H. Damgaard, D. Hochberg, and N. Kawamoto, Phys. Lett. B 158, 239 (1985); P. H. Damgaard, N. Kawamoto, and K. Shigemoto, Nucl. Phys. B 264, 1 (1986).
- [51] H. Kluberg-Stern, A. Morel, B. Petersson, Nucl. Phys. B 215 [FS7], 527 (1983);
- [52] K. Fukushima, Prog. Theor. Phys. Suppl. 153, 204 (2004); Y. Nishida, Phys. Rev. D 69, 094501 (2004).
- [53] G. Faldt and B. Petersson, Nucl. Phys. B 265, 197 (1986).
- [54] N. Bilic, K. Demeterfi and B. Petersson, Nucl. Phys. B 377, 651 (1992).
- [55] N. Bilic, F. Karsch, and K. Redlich, Phys. Rev. D 45, 3228 (1992); N. Bilic and J. Cleymans, Phys. Lett. B 355, 266 (1995).
- [56] T. Jolicoeur, H. Kluberg-Stern, M. Lev, A. Morel, and B. Petersson, Nucl. Phys. B 235, 455 (1984).
- [57] I. Ichinose, Phys. Lett. B 135, 148 (1984); ibid. B 147, 449 (1984); I. Ichinose, Nucl. Phys. B 249, 715 (1985).
- [58] S. Aoki, Nucl. Phys. B 314, 79 (1989).
- [59] K. Miura, T. Z. Nakano and A. Ohnishi, Prog. Theor. Phys. 122, 1045 (2009); K. Miura, T. Z. Nakano, A. Ohnishi and N. Kawamoto, Phys. Rev. D 80, 074034 (2009).
- [60] T. Z. Nakano, K. Miura and A. Ohnishi, Prog. Theor. Phys. 123, 825 (2010).
- [61] T. Z. Nakano, K. Miura, and A. Ohnishi, Phys. Rev. D 83 (2011) 016014.
- [62] K. Fukushima, Phys. Lett. B 591, 277 (2004) [hep-ph/0310121]; C. Ratti, M. A. Thaler and W. Weise, Phys. Rev. D 73, 014019 (2006) [hep-ph/0506234]. K. Morita, V. Skokov, B. Friman and K. Redlich, Phys. Rev. D 84, 076009 (2011) [arXiv:1107.2273 [hep-ph]].
- [63] F. Karsch, S. Ejiri and K. Redlich, Nucl. Phys. A 774, 619 (2006) [hep-ph/0510126]; S. Ejiri, F. Karsch, and K. Redlich, Phys. Lett. B 633, 275 (2006).
- [64] F. Karsch and K. Redlich, Phys. Lett. B 695, 136 (2011).
- [65] B. Stokic, B. Friman and K. Redlich, Phys. Lett. B 673, 192 (2009) [arXiv:0809.3129 [hep-ph]].
- [66] M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi and E. Vicari, Phys. Rev. B 63, 214503 (2001) [cond-mat/0010360].
- [67] J. Engels and F. Karsch, Phys. Rev. D 85, 094506 (2012).
- [68] J. Engels and F. Karsch, Phys. Rev. D 90, 014501 (2014).
- [69] J. Engels and T. Mendes, Nucl. Phys. B 572, 289 (2000); J. Engels, S. Holtmann, T. Mendes and T. Schulze, Phys. Lett. B 492, 219 (2000); J. Engels, S. Holtmann, T. Mendes and T. Schulze, Phys. Lett. B 514, 299 (2001); A. Cucchieri, J. Engels, S. Holtmann, T. Mendes and T. Schulze, J. Phys. A 35, 6517 (2002);
- [70] R. A. Tripolt, J. Braun, B. Klein and B. J. Schaefer, Phys. Rev. D 90, 054012 (2014); Phys. Rev. C 91, 041901 (2015); A. Bhattacharyya, R. Ray and S. Sur, Phys. Rev. D 91, 051501 (2015).
- [71] X. Pan, L. Chen, X. S. Chen and Y. Wu, Nucl. Phys. A 913, 206 (2013).
- [72] B. Friman, Nucl. Phys. A 928, 198 (2014) [arXiv:1404.7471 [nucl-th]].
- [73] J. Braun and B. Klein, Eur. Phys. J. C 63, 443 (2009) [arXiv:0810.0857 [hep-ph]]; J. Braun, B. Klein and P. Piasecki, Eur. Phys. J. C 71, 1576 (2011) [arXiv:1008.2155 [hep-ph]]; J. Braun, B. Klein and B. J. Schaefer, Phys. Lett. B 713, 216 (2012) [arXiv:1110.0849 [hep-ph]]; P. Springer and B. Klein, arXiv:1506.00909 [hep-ph].
- [74] For example, I. Montvay and G. Munster, Quantum Fields on a Lattice (Cambridge University Press, Cambridge, 1994)