Behavior and finite-size effects of the sixth order cumulant in the three-dimensional Ising universality class

The high-order cumulants of conserved charges are suggested to be sensitive observables to search for the critical point of Quantum Chromodynamics (QCD). This has been calculated to the sixth order in experiments. Corresponding theoretical studies on the sixth order cumulant are necessary. Based on the universality of the critical behavior, we study the temperature dependence of the sixth order cumulant of the order parameter using the parametric representation of the three-dimensional Ising model, which is expected to be in the same universality class as QCD. The density plot of the sign of the sixth order cumulant is shown on the temperature and external magnetic field plane. We found that at non-zero external magnetic field, when the critical point is approached from the crossover side, the sixth order cumulant has a negative valley. The width of the negative valley narrows with decreasing external field. Qualitatively, the trend is similar to the result of Monte Carlo simulation on a finite-size system. Quantitatively, the temperature of the sign change is different. Through Monte Carlo simulation of the Ising model, we calculated the sixth order cumulant of different sizes of systems. We discuss the finite-size effects on the temperature at which the cumulant changes sign.


Introduction
Quantum Chromodynamics (QCD) is generally believed to be the theory of describing the strong interaction in terms of quarks and gluons. It is expected that the QCD system undergoes a first order phase transition from hadronic matter to quark-gluon plasma (QGP) at high baryon density and low temperature [1][2][3][4][5][6]. The first order phase transition ends at a critical point, which is a unique character of the QCD phase diagram on the baryon chemical potential-temperature (µ B -T ) plane. Beyond the critical point, it is expected to be a rapid crossover. This has been proved at vanishing baryon chemical potential by lattice QCD [7]. One of the main goals of heavy-ion collision experiments is to locate the QCD critical point and determine the energy scale of QCD phase transitions.
In the thermodynamic limit, the correlation length (ξ) is infinite at the critical point. Susceptibility should diverge as ξ 2 . But in a finite system, ξ does not diverge and the susceptibility is rounded to a finite peak [8,9].
As a system created in heavy-ion collisions, the magnitude of the correlation length is also limited by the system size but mainly limited by the finite-size effects due to critical slowing down [10,11]. A more sensitive observable than the susceptibility is needed in order to search for the critical point. In recent years, the high-order cumulants (this is a mathematical term, corresponding to generalized susceptibilities in physics) of conserved charges are suggested to be sensitive observables to search for the QCD critical point in relativistic heavy-ion collisions, e.g., the net baryon number, net electric charge, and net strangeness [12][13][14][15][16][17]. They are particularly sensitive to the correlation length (ξ). For example, the third order cumulant diverges with ξ 4.5 , and the fourth order cumu-lant diverges with ξ 7 [14].
There are already many works on the high-order cumulants of conserved charges from both experiments and theory [15][16][17][18][19][20][21][22][23]. The theoretical results show there exists non-monotonic or sign change behavior in the vicinity of the critical point in the high-order cumulants. In Ref. [24], the authors found that the behavior of the first four orders of cumulants of net-proton number at the Relative Heavy-Ion Collider (RHIC) calculated by the STAR Collaboration [23] agrees well with the hadron resonance gas (HRG) model and lattice QCD calculations. Based on the three-dimensional O(4) scaling function, the sixth order cumulant of baryon number deviates considerably from the HRG model. It remains negative at the chiral transition temperature [17]. The preliminary results of the sixth order cumulant in STAR at RHIC have been shown in Ref. [25]. It is necessary to study the behavior of the sixth order cumulant from theory.
If the QCD critical point exists, it belongs to the same universality class as the three-dimensional Ising model [26][27][28][29]. They have the same values for the critical exponents. The behavior of the corresponding thermodynamic quantities is similar. That is to say the sixth order cumulant of the order parameter in the threedimensional Ising model can reflect the properties of the sixth order cumulant of the net-baryon.
In Ref. [30], through the Monte Carlo simulation of the Ising model, we have studied the temperature dependence of the sixth order cumulant of the order parameter at three fixed external magnetic fields in a finite system. We found that when approaching the critical point from the crossover side, the sixth order cumulant is negative.
In this paper, using the parametric representation, we further study the distribution of the sign of the sixth order cumulant of the order parameter in the temperatureexternal magnetic field (t-H) plane, and its temperature dependence at a fixed external field in the thermodynamic limit. The qualitative behavior is compared to the results from the Monte Carlo simulations.
As the magnitude of ξ is limited in the system created in heavy ion collisions, the finite-size effects cannot be neglected. In a system of finite-size, we can just get the pseudo-critical point. Within the linear sigma model with constituent quarks and a two-flavor quarkmeson model, the authors found that the pseudo-critical point decreases in the T -µ plane as the system size decreases [31][32][33][34]. In this paper, we study the finite-size effects on the temperature at which the sixth order cumulant changes sign on different sizes of lattice by Monte Carlo simulations.
The paper is organized as follows. In Section 2, the Ising model and its parametric representation are introduced. The parametric expression of the sixth order cumulant is derived. In Section 3, the distribution of the sign of the sixth order cumulant is presented. Its temperature dependence at a fixed external magnetic field is compared with the simulation results of the Monte Carlo method from the three-dimensional Ising model. In Section 4, the finite-size effects on the sixth order cumulant are discussed. The results of the sixth order cumulant are compared with the fourth order one. Finally, the conclusions and summary are given in Section 5.

The Ising model and parametric expression of the sixth order cumulant
The three-dimensional Ising model is defined as follows, where H is the Hamiltonian, and s i is the spin at site i on a cubic lattice which can take only two values ±1. J is the interaction energy between nearest-neighbor spins i, j , and H is the external magnetic field. The magnetization M (the order parameter) is where s = i s i and V = L d denote the total spin and volume of the lattice, respectively, where d = 3 is the dimension of the lattice. The magnetization is a function of the external magnetic field H and the reduced temperature t = (T −T c )/T c , where T c is the critical temperature. The reduced temperature describes the distance away from the critical point. At t > 0, it is the crossover side. At t < 0, it is the first order phase transition side. The high-order cumulants of the order parameter can be obtained from the derivatives of magnetization with respect to the external magnetic field H at fixed t, In particular, the sixth order cumulant is as follows, where δs = s − s .
To study the behavior of the high-order cumulants, one way is to simulate the Ising model by the Monte Carlo method. Another way is to use a parametric representation of the Ising model. In this representation, the magnetization M and reduced temperature t can be parameterized by two variables R and θ [35,36], The equation of state of the three-dimensional Ising model can be given by the parametric representation in terms of R and θ as where β and δ are critical exponents of the three-dimensional Ising universality class with values 0.3267(10) and 4.786 (14), respectively [37]. If M , t and h are analytic functions of θ, the analytic properties of the equation of state are satisfied [38]. The analytic expression of the high-order cumulants can be derived in the parametric representation.
In fact, the function h(θ) is analytic and an odd function of θ because the magnetization is an odd function of the magnetic field M (−H) = −M (H). h(θ) = 0 at θ = 0 corresponds to the crossover line H = 0, T > T c . It vanishes for another θ = θ 0 corresponds to the coexistence curve H = 0, T < T c (the first order phase transition line).
One simple function of h(θ) obeying all the demands is as follows, This is an exact representation of the equation of state of the Ising model to order ε 2 , where ε is a parameter related to the number of dimensions of space. ε-expansion is one of the techniques to explore the critical phenomena. It is enough for our purpose, although the parametric representation can be exact up to order ε 3 [36].
On the other hand, there is an excellent agreement of the scaling magnetization data from Monte Carlo simulations and the equation of state in parametric representations [39]. Using Eq. (4) to Eq. (7), the sixth order cumulant can be expressed in the parametric representation. For our purposes, it would be enough to take the approximate values 1/3 and 5 for the critical exponents β and δ, respectively. Then the sixth order cumulant can be expressed as follows, where a 2n are the expansion coefficients. Their values are listed in Table 1. The main difference in the two ways is that the results of the high-order cumulants are from finite systems by the Monte Carlo simulation, while in the parametric representation, the results are under the condition of the thermodynamic limit.

Behavior of the sixth order cumulant in the parametric representation
In Eq. (8), the sixth order cumulant has two zero values at θ = 0.5723 and θ = 0.183, respectively. Combining Eqs. (5), (6) and Eq. (7), the reduced temperature t can be expressed as a function of θ and the external magnetic field H, Then we can get the sign distribution of κ 6 on the t-H plane. We show it as a density plot in Fig. 1(a). In Fig. 1(a), the yellow point at H = 0 and t = 0 indicates the critical point of the three-dimensional Ising model. On the black curves, the values of the sixth order cumulant are zeros. The middle two black curves 093104-3 correspond to θ = 0.183. The other two correspond to θ = 0.5723. In the red regions, the sixth order cumulant takes negative values. It is positive in the purple region. Specially, the red region is separated into two parts by the purple region, which is different from the sign distribution of the fourth order cumulant given in Ref. [16].
Although the values of the sixth order cumulant are positive in the smaller purple region between the two red regions, it approximately equals to zero.
Along the green dashed line in Fig. 1(a), we study the temperature dependence of the sixth order cumulant at fixed external magnetic field, as shown in Fig. 1(b). The blue line shows where the sixth order cumulant is positive and the red line shows where it is negative. It is clear that the sixth order cumulant is positive but approximately zero at the higher temperature side, which corresponds to the green dashed line above the red region in Fig. 1(a). Then as the temperature decreases and approaches the critical point, the sixth order cumulant has a negative valley.
At fixed external magnetic field, the qualitative temperature dependence of the sixth order cumulant is similar to the parametric representation and Monte Carlo simulations [30]. Quantitatively, in different sizes of systems, due to the finite-size effects, the temperature of the sign change of the sixth order cumulant will be different.

Finite-size effects on the sixth order cumulant
Using Monte Carlo simulations, we study the finitesize effects on the temperature of the sign change of the sixth order cumulant. The simulations are performed on three different sizes of lattices with L = 12, 14, 16 at 11 different external magnetic fields between H = 0.005 and H = 0.04 by the Wolff cluster algorithm [40]. The helical boundary conditions are used. Because of the symmetry of the spin inversion, with negative external magnetic fields the values of the temperatures of the sign change are just the same as those of the corresponding positive ones.
On the t-H plane, we show the positions of the first sign changes at lower temperatures of the sixth order cumulant (e.g., the first zero at lower temperature in Fig. 1(b)) by the points, as shown in Fig. 2. The positions of the first sign changes are closer to the critical point and easier to be distinguished in the Monte Carlo simulations. The black points, blue squares and red stars present the results on lattices of size L = 12, 14, 16, respectively. Specially, pink triangles present the results from the parametric expression of the sixth order cumulant in the parametric representation.
In Fig. 2, for each size, below the points, the values of the sixth order cumulant are positive; above the points, the values are negative. It is clear that, at a fixed external magnetic field, the temperature of sign change increases as the size increases. The amplitude of increase decreases with increasing H. The bigger the external magnetic field, the smaller of the size of the system to get saturation. That is to say, the features of the thermodynamical quantities do not depend on the system size when it is bigger than some value. In fact, when it deviates from the critical point, the correlation length decreases. The system size needed to get saturation decreases, too. According to the results from the parametric expression of the sixth order cumulant, i.e. the pink triangles in Fig. 2, as the external field decreases to zero, the temperature of the sign change reduces to the critical temperature. At the critical point, the sixth order cumulant is divergent. In a finite system of size L = 12, as the external field decreases, the temperature of the sign change decreases. According to the trend, it can be inferred that the temperature of the sign change is smaller than the critical temperature at vanishing external field. As the system size increases from L = 12 to L = 16, the pseudocritical temperature (here we mean the temperature of the sign change at H = 0 in a finite system) increases. When reaching the thermodynamical limit, the pseudocritical temperature will approach the real critical temperature. The qualitative result is consistent with the linear sigma model with constituent quarks and the two flavor quark-meson model [31][32][33][34], which inferred that the QCD pseudo-critical point should shift to the higher temperature side as the size of the system increases.
Comparing the behavior of the sixth order cumulant with the fourth order one, e.g., Fig. 1(b) in this paper and Fig. 1(b) in Ref. [16], qualitatively, at a fixed H, their trends varying with temperature are similar. When the temperature is approached from the higher temperature side, the sixth and fourth order cumulants are both negative. Quantitatively, the negative valley is more obvious in the sixth order cumulant than the fourth order one. For a constant H, the ratio of the maximum to the minimum of the sixth order cumulant is a constant. According to Eq. (8), it approximates −6. For the fourth order cumulant, the ratio is approximately −28. From this result, the benefit of using the sixth order cumulant to search for the critical point is obvious.
On the other hand, a high price is needed to calculate the sixth order cumulant. First, the higher the order of the cumulant, the more statistics are needed. Second, in the vicinity of the critical point, the n-th order cumulant depends on the system size by L n(d−β/ν) [41]. The higher the order of the cumulant, the more serious the finite-size effects.

Summary
Taking advantage of the parametric representation of the three-dimensional Ising model, we have studied the sign distribution of the sixth order cumulant of the order parameter in the t-H plane. At a fixed external magnetic field, we calculated the temperature dependence of the sixth order cumulant. We found that at nonzero external magnetic field, when the critical point is approached from the higher temperature side, the sixth order cumulant has a negative valley. The width of the negative valley narrows with decreasing external field.
Through Monte Carlo simulations, the finite-size effects on the temperature of the sign change of the sixth order cumulant have been studied. At the same external magnetic field, the temperature of the sign change of the sixth order cumulant increases with the increasing size of the system in the Ising universality class. This indicates that the pseudo-critical point should shift towards higher temperatures with increasing size. The trend is consistent with that of the linear sigma model with constituent quarks and the two flavor quark-meson model.
Results of the sixth order cumulant are compared with the fourth order one. At fixed external magnetic field, qualitatively, their temperature dependence is similar. Quantitatively, the negative valley in the sixth order cumulant is more obvious. At the same time, the finitesize effects are also more serious in the higher order of cumulant.