# Cumulants and ordering of their ratios in 2D Potts models: Lessons for QCD?

Rajiv V. Gavai1, Bedangadas Mohanty2,3, Jaydev Singh Rao1 and Swati Saha2,3
1Department of Physics, Indian Institute of Science Education and Research
Bhopal, Bhopal 462066, Madhya Pradesh, India<EMAIL_ADDRESS>2School of
Physical Sciences, National Institute of Science Education and Research, Jatni
752050, Odisha, India<EMAIL_ADDRESS>3Homi Bhabha National Institute,
Training School Complex, Anushaktinagar, Mumbai 400094, Maharastra, India

###### Abstract

Theoretical considerations suggest an ordering of the ratios of net-baryon
number fluctuations in the vicinity of the transition from the low-temperature
hadronic phase to the high temperature quark-gluon plasma phase at small
values of the baryon chemical potential, $\mu_{B}$, in the QCD phase diagram.
The ordering hierarchy is
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$,
where $\chi_{n}$ is the $n^{\mathrm{th}}$ order cumulant of net-baryon number
fluctuation. The STAR experiment observed this hierarchy in the ordering of
cumulant ratios of net-proton number (a proxy of net-baryon number) for a
range of colliding energies. These inequalities can be tested in spin models
by taking the corresponding order parameters in the model as an analog of
baryon density. We employed two different models: the two-state and three-
state Potts models in two dimensions, which undergo a transition from an
ordered phase to a disordered phase at their respective critical temperature.
Simulations were performed on square lattices of different sizes using the
Wolff algorithm. The cumulants of total magnetization are obtained up to the
sixth order in both of these models in a temperature range near their
corresponding critical temperatures. With increasing lattice size, height
(trough) of the peaks (dips) of the higher-order cumulants appears to increase
with the increase in the order of the cumulants. Except in a narrow range
above the critical temperature of the three-state Potts model, the complete
inequality or its complete reverse is not satisfied in the temperature ranges
simulated.

##  I Introduction

Understanding the phase structure and the nature of phase transitions in the
temperature-baryon density (or equivalently the baryon chemical potential,
$\mu_{B}$) plane for the strongly interacting matter, known as quantum
chromodynamics (QCD) phase diagram, is crucial for unraveling the fundamental
properties of matter [1, 2]. While the QCD phase diagram is not well-explored
for larger $\mu_{B}$, the region below $\mu_{B}$ = 450 MeV, spanning down to
$\mu_{B}=0$, has been extensively studied, both in theory [2, 3] and in
experiments [4]. For vanishing and small values of $\mu_{B}$, the transition
from the low-temperature hadronic region to the quark-gluon plasma at high
temperature is well-known to be a smooth transition (cross over) [5]. At
larger values of $\mu_{B}$, however, it is generally expected that a first-
order phase transition line exists, which ends in a second-order critical
point belonging to the three-dimensional Ising model universality class [6, 7,
8]. Comprehensive studies of net-proton fluctuations have been carried out
with the Beam Energy Scan (BES) program of the Relativistic Heavy Ion collider
(RHIC) in the STAR experiment to search for the QCD critical point [9, 10].
Recently, an ordering of the cumulant ratios of net-proton number,
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$,
has been reported by the STAR experiment across the energy range 7.7 to 200
GeV [11] where $\chi_{n}$ is the $n^{\mathrm{th}}$ order cumulant of event-by-
event net-proton number distribution. This ordering was suggested previously
by the Lattice-QCD (LQCD) calculations for the cumulants of net-baryon number
fluctuations at small values of $\mu_{\mathrm{B}}$ [12, 13]. Functional
renormalization group (FRG) calculations were also found to follow this
hierarchy of the inequalities [14]. In contrast, predictions from the Hadron
Resonance Gas (HRG) model, utilizing an ideal gas equation of state within a
grand canonical ensemble framework, consistently yield positive unity for all
cumulant ratios, lacking the observed ordering [15]. One may wonder whether
these inequalities are generically valid in the vicinity of any phase
transition or even in wider regions surrounding it. In the QCD phase diagram,
either near the cross over or near the QCD critical point, if it exists,
whether such inequalities exist is an interesting question to ask.

With the advent of lattice gauge theories, the Z(N) spin models have gained
significant importance in the study of critical behavior in SU(N) gauge
theories at finite temperature. Here, the variable N in general represents the
number of degrees of freedom: within SU(N) gauge theories, it denotes the
number of colors comprising the gauge group, while in Z(N) spin models, it
represents the number of distinct spin states assignable to each lattice site.
Among the Z(N) spin models, the three-state Potts model (Z(3)) in three
dimensions (3D) holds particular significance as it provides a valuable
framework for understanding the phase structures of QCD based on the SU(3)
gauge symmetry group. The Polyakov loop, which serves as an order parameter of
$\text{SU}(3)$ gauge theories, behaves similarly as spins in the low
temperature phase of the 3D three-state Potts model [16, 17]. The line of
first order phase transition in QCD phase diagram [18, 19] is closely related
to the first order phase transition of 3D three-state Potts model in the
absence of an external magnetic field [20, 17]. Addition of a small magnetic
field ($h$) leads to a critical point in the Ising universality class for the
three-state Potts model in 3D at $h_{c}$ [21], making it a nice toy model to
explore the hierarchy of inequalities of the ratios of net-baryon number
cumulants. As a first step towards such an investigation, we studied the
higher-order cumulants of magnetization in the two-state (Z($2$)) and the
three-state (Z($3$)) Potts model on a two-dimensional lattice in this paper.
Both are known to have critical points albeit of another universality class.
The finite-size scaling of the higher-order cumulants of magnetization were
investigated near the critical region in the respective models and the
ordering of their ratios were tested.

The paper is structured as follows: in Sec. II, we provide an overview of the
$q$-state Potts models, detailing their formalism and the computation of
cumulants of magnetization up to $6^{\mathrm{th}}$ order in these models. We
also discuss the application of finite size scaling relation, which enables
the extraction of exponents associated with the scaling of higher-order
cumulants. Section III presents the methodology employed for simulating the
models and outlines the simulation validation process. In Sec. IV, we present
the temperature-dependent behavior of cumulants and their ratios. The
exponents involved with the finite-size scaling of the higher-order cumulants
are estimated in terms of known critical exponents of the model and compared
with the numerical results. Finally, Sec. V summarizes the key findings of our
study.

##  II  Model and Observables

The Hamiltonian for the the $q$-state Potts model [22] is given by

$\mathcal{H}=-\mathcal{J}\sum_{\langle
i,j\rangle}\delta(\sigma_{i},\sigma_{j})-h\sum_{i}\delta(\sigma_{i},0)\,,$ (1)

where $\mathcal{J}$ is a coupling constant in units of energy, determining the
interaction strength among neighbouring spins, $h$ is the external magnetic
field apllied to the model and $\delta$ is the Kronecker delta function. Here,
$\sigma_{k}$ can take any integer values in the set {0, 1, 2, …, ($q$ \- 1)}
and for each $\sigma_{k}$, the spin orientations in the Potts model are
represented by $s_{k}=\text{exp}\left(\frac{i2\pi\sigma_{k}}{q}\right)$. The
sum is over nearest neighbors $\langle i,j\rangle$. In $d$-dimensions there
are 2$d$ such pairs. We focus on $d=2$ in this work. In the absence of an
external magnetic field ($h$ = 0), the model exhibits an exact symmetry
represented by the group of permutations $S_{q}$, which encompasses all
possible global permutations of the $q$ spin values. This symmetry implies
that each spin state is equally likely, and there is no preferred direction
for alignment. However, below a critical temperature, the symmetry is
spontaneously broken resulting in the alignment of spins along a specific
direction and the establishment of an ordered phase. To measure the extent of
symmetry breaking in the system, a quantity called the order parameter is
used, that provides a measure of the overall alignment of spins in the model.
The order parameter for the $q$-state Potts model ($m_{\text{Potts}}$) that we
use in this study is given by

$m_{\text{Potts}}=\frac{q}{q-1}\left(\frac{\text{max}(N_{0},N_{1},...,N_{q-1})}{N}-\frac{1}{q}\right)\,,$
(2)

where $N_{\alpha}=\sum_{j=1}^{L^{d}}\delta(\sigma_{i},\alpha)$ with
$\alpha\in\\{0,1,2,...,q-1\\}$ represents the number of sites in the lattice
with $\sigma_{i}=\alpha$ and $N=L^{d}=\sum_{\alpha}N_{\alpha}$ is the total
number of spin sites in the lattice.

In the following analysis and results, we will be dealing Potts model with
both $q$ = 2 and $q$ = 3 states. The two-state Potts model is equivalent to
the well-known Ising model for which the exact solution by Onsager gives the
critical point at
$T_{C}=\frac{2}{\text{ln}(1+\sqrt{2})}\,\mathcal{J}/k_{\text{B}}\approx
2.2692\,\mathcal{J}/k_{\text{B}}$ [23]. The three-state Potts model that
belongs to the same universality class as the two-dimensional Ising model,
undergoes a second order phase transition at
$T_{C}=\frac{1}{\text{ln}(1+\sqrt{3})}\,\mathcal{J}/k_{\text{B}}\approx
0.9950\,\mathcal{J}/k_{\text{B}}$ [22].

###  II.1  Cumulants

Figure 1: (a) The second order susceptibility of magnetization for the two-
dimensional Ising model (equivalent to $q$ = 2 state Potts model) shown as a
function of temperature. The results shown by different colored markers
represent different lattice sizes ($L$). (b) The reduced temperature,
$t^{\prime}=T(\chi_{2}^{\mathrm{max}})-T_{C}$ is plotted as a function of
$1/L$, and fitted with the function $y=Ax^{1/\nu}$ to extract the value of
$\nu$.

In this paper, the main objects of interest are the higher-order cumulants
(also referred as susceptibilities in the paper), which reflect the higher-
order correlations amongst the spins. In spin models, the $p^{\text{th}}$
order susceptibility normalized by the total number of spin sites in the
lattice is defined [24] by

$\chi_{p}(T)=\frac{1}{L^{d}\beta^{p}}\left(\frac{\partial^{p}\text{ln}Z}{\partial
h^{p}}\right)_{h\to 0}\,,$ (3)

where $\beta=1/(k_{\text{B}}T)$, $T$ is the temperature, $k_{\text{B}}$ is the
Boltzmann constant. Here, $Z$ is the partition function of the model given by
its corresponding $\mathcal{H}$ as

$Z=\sum_{a}\text{exp}(-\beta\mathcal{H})\,.$ (4)

The sum over $a$ in Eq. 4 means the summation over all possible spin
configurations on the lattice. The higher-order spin-spin susceptibilities can
be associated with the cumulants of total magnetization ($M$). This is because
$Z$ is proportional to the moment generating function [25] of $M$ and therfore
one can easily deduce,

$\left(\frac{\partial^{p}\mathrm{ln}Z}{\partial h^{p}}\right)_{h\to
0}=\beta^{p}K_{p}(M)\,.$ (5)

Here, $M$ is referred as magnetization which is equal to the product of the
order parameter and the volume of the lattice, i.e.,
$M=L^{d}m_{\mathrm{Potts}}$. $K_{p}(M)$ is the $p^{\text{th}}$ order cumulant
of $M$. The expressions for $\chi_{p}$ (obtained by using the relation of
cumulants in terms of central moments [26]) up to $p=6$ is listed below,

$\displaystyle\chi_{1}$ $\displaystyle=\langle M\rangle/L^{d}$ (6)
$\displaystyle\chi_{2}$ $\displaystyle=\mu_{2}/L^{d}$ (7)
$\displaystyle\chi_{3}$ $\displaystyle=\mu_{3}/L^{d}$ (8)
$\displaystyle\chi_{4}$ $\displaystyle=(\mu_{4}-3\mu_{2}^{2})/L^{d}$ (9)
$\displaystyle\chi_{5}$ $\displaystyle=(\mu_{5}-10\mu_{3}\mu_{2})/L^{d}$ (10)
$\displaystyle\chi_{6}$
$\displaystyle=(\mu_{6}-15\mu_{4}\mu_{2}-10\mu_{3}^{2}+30\mu_{2}^{3})/L^{d}$
(11)

where $\mu_{p}=\langle(M-\langle M\rangle)^{p}\rangle$ is the $p^{\text{th}}$
order central moment of the magnetization and the angular brackets,
$\langle..\rangle$ represents average over all possible spin configurations on
the lattice. We will be studying these $\chi_{p}$’s as a function of
temperature in the critical region.

###  II.2  Finite size scaling

In this section, we obtain the generalized relation for finite size scaling of
the higher-order susceptibilities of magnetization in the Potts model. We
define, $t=(T-T_{C})/T_{C}$ and $h_{r}=\beta h$. For a second order
(continuous) phase transition, there are at least two relevant scaling fields
for the free energy: temperature like and magnetic field like. The dependence
of free energy on these scaling fields near the critical point can be
established by a bulk Renormalization Group (RG) transformation. RG
transformation involves length scale re-scaling, under which, the scaling
fields generally look like $C_{1}L^{1/\nu}t$ and $C_{2}L^{\Delta/\nu}h_{r}$
respectively [27, 28]. Here, the exponent $\nu$ is defined such that $\xi\sim
t^{-\nu}$ (for a finite lattice $\xi\sim L$ at the critical point) and
$\Delta\equiv\beta+\gamma$ is such that $\chi_{1}\sim t^{\beta}$ and
$\chi_{2}\sim t^{\gamma}$ in the zero field limit, i.e., $h\rightarrow 0$.

In the obtained scaling relation for the singular part of the free energy
after RG transformations, irrelevant scaling fields and non-linear
contributions from relevant scaling fields can be safely removed by setting
them to zero. This simplification is valid under the condition that $d<d_{>}$,
where $d_{>}$ is the upper critical dimension. Consequently, only the linear
order of temperature-like ($t$) and magnetic field-like ($h_{r}$) scaling
fields remains, serving as the sole parameters that determine the distribution
of the system [27]. The asymptotic finite size scaling relation of the
singular part of reduced free energy density, i.e $f_{r}^{(s)}=L^{-d}\ln Z$,
can therefore be written as (for $t\rightarrow 0,h_{r}\rightarrow 0$)

$f_{r}^{(s)}(t,h;L)\approx
L^{-d}Y(C_{1}L^{1/\nu}t,C_{2}L^{\Delta/\nu}h_{r})\,.$ (12)

Here, the function $Y(x,y)$ is same for every system in the universality
class, i.e. it is universal, while the coefficients $C_{1}$ and $C_{2}$ are
generally not universal. The $n^{th}$ derivative of Eq. 12 with respect to
$h_{r}$ gives the scaling relation for the $n^{th}$ order finite-size
susceptibility,

$\chi_{n}^{(s)}(t;L)\approx-
C_{2}^{n}L^{\frac{n\Delta}{\nu}-d}\frac{\partial^{n}Y}{\partial
h_{r}^{n}}\bigg{|}_{h_{r}=0}\,.$ (13)

From Eq. 13, we can extract the scaling relation for the singular part of
$\chi_{n}$, i.e., $\chi_{n}^{(s)}\sim L^{\frac{n\Delta}{\nu}-d}$. To simplify
the exponents, we utilize the relations among the critical exponents:
$\gamma=2\Delta-(2-\alpha)$ and $2-\alpha=d\nu$. Both of these relations are
valid for $d=2$, but in general, the latter is only valid for any $d<d_{>}$.
Thus, we can write the scaling exponent of $n^{\mathrm{th}}$ order
susceptibility as,

$e_{n}\equiv\frac{n\Delta}{\nu}-d=\frac{1}{2}\frac{n\gamma}{\nu}+(n-2)d\,.$
(14)

These exponents for $n$ = 2 to 6 can now be calculated using the respective
values of $\gamma$ and $\nu$ for the model under consideration. For the $q=2$
state Potts model on a two-dimensional lattice, we have $\gamma=\frac{7}{4}$
and $\nu=1$ for $q=2$ [29], while $\gamma=\frac{13}{9}$ and $\nu=\frac{5}{6}$
for $q=3$ [22].

Figure 2: (a) The second order susceptibility of magnetization for the two-
dimensional $q$ = 3 state Potts model shown as a function of temperature. The
results shown by different colored markers represent different lattice sizes
($L$). (b) The reduced temperature, $t^{\prime}=T(C_{2}^{\mathrm{max}})-T_{C}$
is plotted as a function of $1/L$, and fitted with the function $y=Ax^{1/\nu}$
to extract the value of $\nu$.

##  III  Numerical simulation

Both $q=2$ and $q=3$ state Potts model were simulated on two-dimensional
lattices using the Wolff Cluster algorithm [30] at temperatures close to their
corresponding critical temperature ($T_{C}$). In this algorithm, clusters of
similarly oriented spin sites are determined randomly which are then flipped
together to generate new states. This algorithm reduces the effect of critical
slowing down near $T_{C}$ with increase in lattice size, i.e. it has a smaller
dynamical exponent resulting in more efficient simulations in this region
[31]. For $q=2$, square lattices of linear size $L$ = 50, 60, 70, 80 and 90
were used in the simulations and during the simulation $10^{4}$ independent
configurations were generated after ensuring thermalization by throwing away
sufficient initial configurations. Similarly, for $q=3$, simulations were
performed on square lattices with $L$ = 48, 64, 80, 96 and 128 to generate
$2\times 10^{4}$ independent configurations after approprite thermalization.
In order to obtain independent spin configurations on the lattice,
configurations were recorded at the interval of $2\tau$ where $\tau$ is the
auto-correlation time that is obtained from the auto-correlation function (of
$M$) defined in Eq. 15 using an exponential decay ansatz
$\phi(t)\sim\mathrm{exp}\left({-\frac{t}{\tau}}\right)$.

$\phi(t)=\frac{\langle M(0)M(t)\rangle-\langle M\rangle^{2}}{\langle
M^{2}\rangle-\langle M\rangle^{2}}$ (15)

The magnetization is calculated for each spin configuration in the simulation
and the higher-order susceptibilities of magnetization are calculated using
Eqs. 6 \- 11 by averaging over all possible spin configurations on the
lattice.

Figure 3: The third and fourth order susceptibilities of magnetization,
$\chi_{3}$ and $\chi_{4}$, are shown for different lattice sizes as a function
of temperature scaled by critical temperature of the corresponding lattice,
$T/T_{C,L}$. The left side plot belongs to two-dimensional $q=2$ state Potts
model whereas the right side plot is for $q=3$ state Potts model. The upper
panel represents $\chi_{3}$ and the bottom panel represents $\chi_{4}$ [where
$\chi_{3}$ and $\chi_{4}$ divided by a factor $10^{5}$ ($10^{6}$) and $10^{8}$
($10^{9}$) respectively for $q=2$ ($q=3$)].

##  IV  Results

The top panels of Fig. 1 and Fig. 2 display the second-order susceptibility,
denoted as $\chi_{2}$, for the $q=2$ and $q=3$ state Potts model respectively
with varying lattice sizes, $L$. In both cases, $\chi_{2}$ exhibits a
distinctive peak structure. As may be expected, the locations ($T_{C,L}$) and
magnitudes ($\chi_{2}(T_{C,L})$) of these peaks vary with changes in the
lattice size. Their scaling behavior can be utilized to extract the critical
exponents, specifically $\gamma$ and $\nu$, which will be further elaborated
upon in the subsequent section, IV.1.

Figure 4: The fifth and sixth order susceptibilities of magnetization,
$\chi_{5}$ and $\chi_{6}$, are shown for different lattice sizes as a function
of temperature scaled by critical temperature of the corresponding lattice,
$T/T_{C,L}$. The left side plot belongs to two-dimensional $q=2$ state Potts
model whereas the right side plot is for $q=3$ state Potts model. The upper
panel represents $\chi_{5}$ and the bottom panel represents $\chi_{6}$ [where
$\chi_{5}$ and $\chi_{6}$ divided by a factor $10^{12}$ ($10^{13}$) and
$10^{16}$ ($10^{17}$) respectively for $q=2$ ($q=3$)].

###  IV.1  Critical exponents

The critical exponents, $\gamma$ and $\nu$ can be obtained using the scaling
relations, $\chi_{2}(T_{C,L})\propto L^{-\gamma/\nu}$ and
$T_{C,L}-T_{C}\propto L^{-1/\nu}$, where $T_{C,L}$ is the psuedo-critical
temperature for a given lattice size, and $T_{C}$ represents the critical
temperature of the model in the limit $L\to\infty$. The pseudo-critical
temperature is an effective critical temperature observed in finite-sized
systems, reflecting behavior similar to the critical point in infinite-sized
systems. $T_{C,L}$ is defined as the location of the peak in the second order
susceptibility, $\chi_{2}$ distribution as shown in Figs. 1 and 2. By fitting
$\chi_{2}$ as a Gaussian function of the temperature, $T$, the peak location
($T_{C,L}$, $\chi_{2}(T_{C,L})$) is determined for each lattice size. The
obtained values of $T_{C,L}$ and $\chi_{2}(T_{C,L})$ as a function of $L$ are
tabulated in Tables 1 and 2 for $q$ = 2 and $q$ = 3 state Potts model
respectively. Also, we define the critical region around $T_{C,L}$ for each
lattice size by the width of the $\chi_{2}$ distribution peak, which is equal
to the $\sigma$ of the fitted Gaussian distribution function, represented here
as $\sigma_{C}$. Using the values of $T_{C,L}$, the difference
$t^{\prime}=T_{C,L}-T_{C}$ is plotted as a function of $1/L$, and fitted with
its scaling relations for both $q=2$ and $q=3$ state Potts model as shown in
the bottom panels of Figs. 1 and 2 respectively. Alternatively, $T_{C}$ can be
estimated by plotting directly $T_{C,L}$ against $L$ and fitting with
$y=A+Bx^{-1/\nu}$. The resulting estimates ($T_{C}=2.2780\pm 0.0089$ for $q=2$
and $T_{C}=0.9953\pm 0.0001$ for $q=3$) are found to be consistent with the
known values [23, 22] within uncertainties. For precise measurements of the
critical exponents, however we employ the known values of $T_{C}$ in the
former method and extracted the coefficient $\nu$. The peak position,
$\chi_{2}(T_{C,L})$ as a function of $L$, is also fitted with its scaling
relation. Therefore, the extracted values of $\gamma$ and $\nu$ from the fits
with the scaling relations are as follows: for $q=2$, $\gamma=1.758\pm 0.060$,
$\nu=1.019\pm 0.034$ whereas for $q=3$, $\gamma=1.435\pm 0.014$, $\nu=0.830\pm
0.008$. The excellent agreement observed between the calculated values of
$\gamma$ and $\nu$ and their theoretical counterparts affirms the validation
of our simulation.

Lattice size | $T_{C,L}$ | $\chi_{2}(T_{C,L})$ | $\sigma_{C}$
---|---|---|---
50 | 2.3151 $\pm$ 0.0007 | 105.656 | 0.0495 $\pm$ 0.0011
60 | 2.3076 $\pm$ 0.0006 | 143.913 | 0.0424 $\pm$ 0.0009
70 | 2.3020 $\pm$ 0.0007 | 188.497 | 0.0360 $\pm$ 0.0009
80 | 2.2974 $\pm$ 0.0005 | 238.766 | 0.0310 $\pm$ 0.0007
90 | 2.2956 $\pm$ 0.0005 | 289.739 | 0.0288 $\pm$ 0.0007
Table 1: The extracted values for position ($T_{C,L}$), height ($\chi_{2}(T_{C,L})$) and width ($\sigma_{C}$) of the peak of $\chi_{2}$ distribution as a function of temperature for $q$ = 2 state Potts model. The quoted values after ’$\pm$’ represent statistical uncertainties derived from the fit.  Lattice size | $T_{C,L}$ | $\chi_{2}(T_{C,L})$ | $\sigma_{C}$
---|---|---|---
48 | 1.0039 $\pm$ 0.0001 | 105.656 | 0.0102 $\pm$ 0.0002
64 | 1.0012 $\pm$ 0.0001 | 143.913 | 0.0071 $\pm$ 0.0002
80 | 0.9998 $\pm$ 0.0001 | 188.497 | 0.0056 $\pm$ 0.0001
96 | 0.9989 $\pm$ 0.0001 | 238.766 | 0.0046 $\pm$ 0.0001
128 | 0.9978 $\pm$ 0.0001 | 289.739 | 0.0033 $\pm$ 0.0001
Table 2: The extracted values for position ($T_{C,L}$), height
($\chi_{2}(T_{C,L})$) and width ($\sigma_{C}$) of the peak of $\chi_{2}$
distribution as a function of temperature for $q$ = 3 state Potts model. The
quoted values after ’$\pm$’ represent statistical uncertainties derived from
the fit.

We plot the higher-order susceptibilities, $\chi_{3}$ and $\chi_{4}$ for both
$q=2$ and $q=3$ state Potts model in Fig. 3 as a function of temperature
scaled by $T_{C,L}$. Since, these $\chi_{n}$’s are one order derivative of its
preceding $\chi_{n-1}$ (for $n>2$), the sign of peak structure gets reversed
in $\chi_{n}$ with respect to $\chi_{n-1}$. Similarly, $\chi_{5}$ and
$\chi_{6}$ are shown as a function of $T/T_{C,L}$ for both the models in Fig.
4. The height of the peaks or dips in these higher-order susceptibilities are
found to be lattice size dependent. They also follow a scaling as a function
of $L$ as discussed in section II.2. In order to determine the scaling
exponents of the higher-order susceptibilities, we fit the maximum or minimum
value of $\chi_{n}$, depending on which is prominent, as a function of $L$
against the scaling relation: $\chi_{n}^{\text{max/min}}\propto L^{b}$. The
scaling exponents of $\chi_{3}$, $\chi_{4}$, $\chi_{5}$ and $\chi_{6}$
obtained from the fit are referred as $e_{n}^{{}^{\prime\prime}}$. For $q=2$
state Potts model, the values of scaling exponents are given in Table 3,
whereas for $q=3$ state Potts model, they are provided in Table 4.
Additionally, we compare the scaling exponents (Eq. 14) obtained using exact
values of $\gamma$ and $\nu$ as well as by using the estimates of $\gamma$ and
$\nu$ (from scaling of $\chi_{2}$) obtained in this paper, and they are
referred as $e_{n}$ and $e_{n}^{\prime}$ respectively in Tables 3 and 4. The
numerically estimated scaling exponents, $e_{n}^{\prime\prime}$ for the
higher-order susceptibilities show a good agreement with the analytically
obtained values. We also note that determining the exponents from finite size
scaling of $\chi_{n}$’s instead of using the numerical estimate of $\gamma$
and $\nu$ from $\chi_{2}$ result in smaller numerical errors.

Quantity | $e_{n}$ | $e_{n}^{\prime}$ | $e_{n}^{\prime\prime}$
---|---|---|---
$\chi_{3}$ | 3.625 | $3.566\pm 0.097$ | $3.608\pm 0.026$
$\chi_{4}$ | 5.500 | $5.420\pm 0.130$ | $5.500\pm 0.074$
$\chi_{5}$ | 7.375 | $7.280\pm 0.160$ | $7.257\pm 0.067$
$\chi_{6}$ | 9.250 | $9.130\pm 0.190$ | $9.196\pm 0.066$
Table 3: Finite size scaling exponents of higher-order susceptibilities, $e_{n}^{{}^{\prime\prime}}$ obtained from fits to peak (or dip) values of $\chi_{n}$ as a function of lattice size in two-dimensional $q$ = 2 state Potts model. Here, $e_{n}$ is computed based on exact values of $\gamma$ and $\nu$ ($\gamma=7/4$ and $\nu=1$) in Eq. 14, while $e_{n}^{\prime}$ is derived using our estimates of $\gamma$ and $\nu$ discussed in Sec. IV.1. The quoted values after ’$\pm$’ represent statistical uncertainties derived from the fit.  Quantity | $e_{n}$ | $e_{n}^{\prime}$ | $e_{n}^{\prime\prime}$
---|---|---|---
$\chi_{3}$ | 3.599 | $3.593\pm 0.036$ | $3.619\pm 0.021$
$\chi_{4}$ | 5.466 | $5.458\pm 0.047$ | $5.501\pm 0.047$
$\chi_{5}$ | 7.333 | $7.322\pm 0.059$ | $7.301\pm 0.047$
$\chi_{6}$ | 9.199 | $9.187\pm 0.071$ | $9.197\pm 0.077$
Table 4: Finite size scaling exponents of higher-order susceptibilities,
$e_{n}^{{}^{\prime\prime}}$ obtained from fits to peak (or dip) values of
$\chi_{n}$ as a function of lattice size in two-dimensional $q$ = 3 state
Potts model. Here, $e_{n}$ is computed based on exact values of $\gamma$ and
$\nu$ ($\gamma=13/9$ and $\nu=5/6$) in Eq. 14, while $e_{n}^{\prime}$ is
derived using our estimates of $\gamma$ and $\nu$ discussed in Sec. IV.1. The
quoted values after ’$\pm$’ represent statistical uncertainties derived from
the fit.

Figure 5: Ratios of the susceptibilities, $\chi_{3}/\chi_{1}$,
$\chi_{4}/\chi_{2}$, $\chi_{5}/\chi_{1}$ and $\chi_{6}/\chi_{2}$ for two-
dimensional $q$ = 2 state Potts model with lattice sizes, L = 50 (left plot),
60 (middle plot), and 80 (right plot) are shown as a function of scaled
temperature, $T/T_{C,L}$. The yellow shaded region represents the critical
region, $T_{C,L}\pm\sigma_{C}$.

Figure 6: Ratios of the susceptibilities, $\chi_{3}/\chi_{1}$,
$\chi_{4}/\chi_{2}$, $\chi_{5}/\chi_{1}$ and $\chi_{6}/\chi_{2}$ for two-
dimensional $q$ = 3 state Potts model with lattice sizes, L = 50 (left plot),
60 (middle plot), and 80 (right plot) are shown as a function of scaled
temperature, $T/T_{C,L}$. The yellow shaded region represents the critical
region, $T_{C,L}\pm\sigma_{C}$.

###  IV.2  Ordering of Cumulant ratios

In this section, we delve into the discussion of whether the two-dimensional
$q=2$ and $q=3$ state Potts models exhibit ordering in the susceptibility
ratios,
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$.
Figure 5 presents a comparison of susceptibility ratios,
$\frac{\chi_{3}}{\chi_{1}}$, $\frac{\chi_{4}}{\chi_{2}}$,
$\frac{\chi_{5}}{\chi_{1}}$, and $\frac{\chi_{6}}{\chi_{2}}$, in the $q$ = 2
state Potts model for lattice sizes of 50, 60, and 80 as a function of
temperature scaled by $T_{C,L}$. The shaded yellow band corresponds to the
critical region, $T_{C,L}\pm\sigma_{C}$, where the values of $\sigma_{C}$ are
listed in Table 1. It is already clear from Figs. 3 and 4 that the scale of
the susceptibilities differ by factors larger than $\sim\mathcal{O}(10^{3})$.
Therefore, to investigate the inequalities among the susceptibility ratios (if
it exists) thoroughly, the ratios are displayed pairwise in three panels (for
easy inspection by eye). Furthermore, a meaningful discussion of the ordering
of the ratio of susceptibilities is feasible in the critical region around
$T_{C,L}$ as a result of the vast changes in the magnitudes of successively
higher orders. A clear ordering,
$\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$
is observed for $T>T_{a}$ in the critical region from the top and middle panel
plots of Fig. 5, where $T_{a}\sim 1.01T_{C,L}$ ($T_{a}\sim 1.005T_{C,L}$) for
$L=50$ ($L=80$). But the ordering,
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}$ is not observed for
$T>T_{a}$. Instead, a reverse ordering is observed, with
$\frac{\chi_{5}}{\chi_{1}}$ being smaller than $\frac{\chi_{6}}{\chi_{2}}$ for
the entire $T>T_{C,L}$. Therefore, we find that the expected complete ordering
as seen in [11] is not observed in any temperature range near the critical
region in $q=2$ state Potts model. Only a truncated ordering of the
susceptibility ratios up to $\frac{\chi_{5}}{\chi_{1}}$, i.e.,
$\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$
holds true for $T>1.01T_{C,L}$ across all lattice sizes shown in Fig. 5. Such
a truncated inequality is not observed for the entire temperature range,
$T<T_{a}$.

The susceptibility ratios for the three-state Potts model are similarly
examined for lattice sizes 48, 80 and 128 in Fig. 6 and similar trends are
observed for the ratios as seen in the $q$ = 2 state Potts model. An ordering
of the ratios,
$\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$
is observed for $T>T_{a}$ within the critical region where $T_{a}\sim
1.002T_{C,L}$ ($T_{a}\sim 1.001T_{C,L}$) for $L=48$ ($L=80$). For the largest
lattice size, $L=128$ in Fig. 6 (right plot), the ratios show distinct
differences only in a narrow region around $T/T_{C,L}=1$ between 0.984 to
1.008, making it possible to talk of any inequality in this region around the
transition point. In the top panel of the right plot,
$\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$ for $T/T_{C,L}$ in
(0.996, 1.003) and the inequality is reversed for (0.984, 0.996). Similarly,
in the middle panel, $\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}$ in
(0.99, 0.994) and (0.999, 1.004) while inequality is reversed in (0.994,
0.999). Finally, in the bottom panel,
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}$ in two regions of
$T/T_{C,L}$, (0.993, 0.997) and (1.001, 1.006), whereas the ratio is reversed
in the two complementary regions, (0.986, 0.983) and (0.997, 1.001).
Therefore, it can now be easily summarized that only in a narrow range of
temperature, (1.001$T_{C,L}$, 1.003$T_{C,L}$), a complete ordering of ratios,
i.e.,
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$
is observed. In no other region of temperature, such an inequality or its
complete reversed version is seen, including the entire region below
$T_{C,L}$. As one decreases the lattice size, a similar trend is observed
albeit in yet narrower interval for temperature greater than $T_{C,L}$.

A partial reversal of the ordering up to $\frac{\chi_{5}}{\chi_{1}}$, i.e.,
$\frac{\chi_{5}}{\chi_{1}}>\frac{\chi_{4}}{\chi_{2}}>\frac{\chi_{3}}{\chi_{1}}$
is seen in a narrow temperature range below the $T_{C,L}$, but a complete
reversal is not observed for both $q=2$ and $q=3$ state Potts model. An
important point to note is that regions of inequality exclusively lie in the
critical region, which itself shrinks as in the infinite volume limit. This
suggests that for infinite system size, we may not be able to establish any
inequalities for the susceptibility ratios of magnetization in the spin models
considered, and any observed orderings are a finite size effect. It would be
interesting to check if any such volume dependence of the produced medium is
seen in the relativistic heavy ion collisions.

##  V  Conclusions

LQCD predicts a specific hierarchy of net-baryon number susceptibility ratios
near the quark-hadron transition temperature for small $\mu_{B}$:
$\frac{\chi_{6}}{\chi_{2}}<\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$.
The STAR experiment tested these predictions for net-proton number cumulant
ratios in Au+Au collisions at various center-of-mass energies. Cumulant ratios
from 7.7 to 200 GeV exhibited the predicted hierarchy, while a reverse
ordering was observed at 3 GeV. If these inequalities are universally
applicable for any system undergoing phase transition, is what we tried to
understand.

In this study, we investigated higher-order susceptibilities of magnetization
and their ratios in the critical regions of two-dimensional two and three-
state Potts models. Utilizing the Wolff cluster algorithm to mitigate critical
slowing down effects, simulations were conducted. We employed finite size
scaling techniques to extract scaling exponents associated with higher-order
susceptibilities, extending up to the sixth order. Additionally, we estimated
these scaling exponents using the known critical exponents of the model,
specifically the $\nu$ and $\gamma$ exponents. The obtained values from both
methods exhibited a good agreement, indicating that the behavior of higher-
order susceptibilities aligns well with our understanding of second-order
phase transitions in the spin models.

We observed an absence of a consistent order in susceptibility ratios both
above and below the pseudo-critical point (pseudo-critical temperature). A
truncated ordering,
$\frac{\chi_{5}}{\chi_{1}}<\frac{\chi_{4}}{\chi_{2}}<\frac{\chi_{3}}{\chi_{1}}$
is observed in a small temperature range, greater than the pseudo-critical
temperature for a given lattice size in both two and three-state Potts models.
In a yet narrow temperature range just above the pseudo-critical temperature,
the complete ordering hierarchy in the susceptibility ratios is noted in the
three-state Potts models for $L=128$ (illustrated in Fig. 6). While a complete
reversal of this ordering is not observed in either of the models, we
acknowledge that a partial reversal of the ordering,
$\frac{\chi_{5}}{\chi_{1}}>\frac{\chi_{4}}{\chi_{2}}>\frac{\chi_{3}}{\chi_{1}}$
is seen in a small temperature range below the pseudo-critical temperature.
Furthermore, we found that the region suitable for making comparisons between
these ratios diminishes notably as the lattice size (or “lattice volume”)
increases, shrinking to a point in the infinite volume limit. Consequently,
the orderings observed in the susceptibility ratios can be attributed to
finite size effects in spin model.

The spin models considered here have only a second order phase transition. A
more realistic simulation which will be in better analogy with the
expectations for the QCD phase diagram entails employing $q=3$ state Potts
model on a three-dimensional lattice and with a small but non-vanishing
magnetic field. It will be interesting to see what changes does a first order
phase transition or the cross over region beyond the critical magnetic field
make to the above observations on the inequality structure.

###### Acknowledgements.

We acknowledge the support of the Department of Atomic Energy (DAE) and the
Department of Science and Technology (DST), Government of India. RVG expresses
gratitude for the financial support from DAE through the Raja Ramanna
Fellowship and the hospitality of the Department of Physics at IISER Bhopal.
BM and SS acknowledge the financial support of DAE and DST. We also
acknowledge the use of Garuda HPC facility at the School of Physical Sciences,
NISER. SS would like to thank V K S Kashyap for technical help in running
simulations, as well as Ashish Pandav and Debasish Mallick for insightful
discussions regarding STAR results.

## References

  * [1] Edward V. Shuryak.  Quantum Chromodynamics and the Theory of Superdense Matter.  Phys. Rept., 61:71–158, 1980.
  * [2] Heng-Tong Ding, Frithjof Karsch, and Swagato Mukherjee.  Thermodynamics of strong-interaction matter from Lattice QCD.  Int. J. Mod. Phys. E, 24(10):1530007, 2015.
  * [3] Szabocls Borsanyi, Zoltan Fodor, Christian Hoelbling, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo.  Full result for the QCD equation of state with 2+1 flavors.  Phys. Lett. B, 730:99–104, 2014.
  * [4] Anton Andronic, Peter Braun-Munzinger, Krzysztof Redlich, and Johanna Stachel.  Decoding the phase structure of QCD via particle production at high energy.  Nature, 561(7723):321–330, 2018.
  * [5] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo.  The Order of the quantum chromodynamics transition predicted by the standard model of particle physics.  Nature, 443:675–678, 2006.
  * [6] Juergen Berges and Krishna Rajagopal.  Color superconductivity and chiral symmetry restoration at nonzero baryon density and temperature.  Nucl. Phys. B, 538:215–232, 1999.
  * [7] M. A. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. M. Verbaarschot.  Phase diagram of qcd.  Phys. Rev. D, 58:096007, Sep 1998.
  * [8] M. Asakawa and K. Yazaki.  Chiral Restoration at Finite Density and Temperature.  Nucl. Phys. A, 504:668–684, 1989.
  * [9] A. Pandav, D. Mallick, and B. Mohanty.  Search for the QCD critical point in high energy nuclear collisions.  Prog. Part. Nucl. Phys., 125:103960, 2022.
  * [10] Adam Bzdak, Shinichi Esumi, Volker Koch, Jinfeng Liao, Mikhail Stephanov, and Nu Xu.  Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan.  Phys. Rept., 853:1–87, 2020.
  * [11] Bassam Aboona et al.  Beam Energy Dependence of Fifth and Sixth-Order Net-proton Number Fluctuations in Au+Au Collisions at RHIC.  Phys. Rev. Lett., 130(8):082301, 2023.
  * [12] A. Bazavov, D. Bollweg, H.-T. Ding, P. Enns, J. Goswami, P. Hegde, O. Kaczmarek, F. Karsch, R. Larsen, Swagato Mukherjee, H. Ohno, P. Petreczky, C. Schmidt, S. Sharma, and P. Steinbrecher.  Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice qcd confront high-statistics star data.  Phys. Rev. D, 101:074502, Apr 2020.
  * [13] A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Swagato Mukherjee, H. Ohno, P. Petreczky, E. Rinaldi, H. Sandmeyer, C. Schmidt, Chris Schroeder, S. Sharma, W. Soeldner, R. A. Soltz, P. Steinbrecher, and P. M. Vranas.  Skewness and kurtosis of net baryon-number distributions at small values of the baryon chemical potential.  Phys. Rev. D, 96:074510, Oct 2017.
  * [14] Wei-jie Fu, Xiaofeng Luo, Jan M. Pawlowski, Fabian Rennecke, Rui Wen, and Shi Yin.  Hyper-order baryon number fluctuations at finite temperature and density.  Phys. Rev. D, 104:094047, Nov 2021.
  * [15] Szabolcs Borsanyi, Zoltan Fodor, Jana N. Guenther, Sandor K. Katz, Kalman K. Szabo, Attila Pasztor, Israel Portillo, and Claudia Ratti.  Higher order fluctuations and correlations of conserved charges from lattice QCD.  JHEP, 10:205, 2018.
  * [16] Benjamin Svetitsky and Laurence G. Yaffe.  Critical Behavior at Finite Temperature Confinement Transitions.  Nucl. Phys. B, 210:423–447, 1982.
  * [17] R. V. Gavai, F. Karsch, and B. Petersson.  A Study of the Correlation Length Near a First Order Phase Transition: The 3-$D$ Three State Potts Model.  Nucl. Phys. B, 322:738–758, 1989.
  * [18] L. G. Yaffe and B. Svetitsky.  First-order phase transition in the su(3) gauge theory at finite temperature.  Phys. Rev. D, 26:963–965, Aug 1982.
  * [19] F Green and Frithjof Karsch.  Mean field analysis of su (n) deconfining transitions in the presence of dynamical quarks.  Nuclear Physics B, 238(2):297–306, 1984.
  * [20] H. W. J. Blöte and R. H. Swendsen.  First-order phase transitions and the three-state potts model.  Phys. Rev. Lett., 43:799–802, Sep 1979.
  * [21] Frithjof Karsch and Sven Stickan.  The Three-dimensional, three state Potts model in an external field.  Phys. Lett. B, 488:319–325, 2000.
  * [22] F. Y. Wu.  The Potts model.  Reviews of Modern Physics, 54(1):235–268, January 1982.
  * [23] Lars Onsager.  Crystal statistics. i. a two-dimensional model with an order-disorder transition.  Phys. Rev., 65:117–149, Feb 1944.
  * [24] R. K. Pathria and Paul D. Beale.  Statistical Mechanics.  Elsevier, 3rd edition, 2011.
  * [25] R. G. M. Rodrigues, B. V. Costa, and L. A. S. Mól.  Moment-generating function zeros in the study of phase transitions.  Phys. Rev. E, 104:064103, Dec 2021.
  * [26] Peter J Smith.  A recursive formulation of the old problem of obtaining moments from cumulants and vice versa.  The American Statistician, 49(2):217–218, 1995.
  * [27] V. Privman and M. E. Fisher.  Universal critical amplitudes in finite-size scaling.  Phys. Rev. B, 30:322–327, Jul 1984.
  * [28] R.K. Pathria and Paul D. Beale.  14 - phase transitions: The renormalization group approach.  In R.K. Pathria and Paul D. Beale, editors, Statistical Mechanics (Third Edition), pages 539–581. Academic Press, Boston, third edition edition, 2011.
  * [29] Andrea Pelissetto and Ettore Vicari.  Critical phenomena and renormalization group theory.  Phys. Rept., 368:549–727, 2002.
  * [30] Ulli Wolff.  Collective monte carlo updating for spin systems.  Phys. Rev. Lett., 62:361–364, Jan 1989.
  * [31] M. E. J. Newman and G. T. Barkema.  Monte Carlo methods in statistical physics.  Clarendon Press, Oxford, 1999.