# Connecting the Unstable Region of the Entropy to the Pattern of the Fisher’s
Zeros Map

J. C. S. Rocha  Departamento de Física, ICEB, Universidade Federal de Ouro
Preto, CEP 35402-136, Ouro Preto, Minas Gerais, Brazil  B.V. Costa  Retired
Professor, Laboratório de Simulação, Departamento de Física, ICEx
Universidade Federal de Minas Gerais, 31720-901 Belo Horizonte, Minas Gerais,
Brazil

###### Abstract

Phase transitions are one of the most interesting natural phenomena. For
finite systems, one of the concerns in the topic is how to classify a specific
transition as being of first, second, or even of a higher order, according to
the Ehrenfest classification. The partition function provides all the
thermodynamic information about the physical systems, and a phase transition
can be identified by the complex temperature where it is equal to zero. In
addition, the pattern of the zeros on the complex temperature plane can
provide evidence of the order of the transition. In this manuscript, we
present an analytical and simulational study connecting the microcanonical
analysis of the unstable region of the entropy to the canonical partition
function zeros. We show that, for the first-order transition, the zeros
accumulate uniformly in a vertical line on the complex inverse temperature
plane as discussed in previous works. We illustrate our calculation using the
$147$ particles Lennard-Jones cluster.

##  1 Introduction

The transitions between different states of matter observed in macroscopic
systems, such as the solidification of water into ice or the demagnetization
of a magnetic rod, are well-described phenomena. In the case of water
freezing, this transition involves the coexistence of both the liquid and
solid phases, whereas the demagnetization process does not exhibit a distinct
boundary between the ferromagnetic and paramagnetic phases. According to P.
Ehrenfest [1], these phase transitions are classified as being of first and
second order, respectively. This classification is based on appointing the
lowest discontinuous derivative of the free energy at the transition point,
thereby implying a non-analytical point in this function.

In the realm of equilibrium statistical mechanics, the key entity for
investigating a system is its partition function. The canonical partition
function, denoted as $Z(\beta,V,N)$, holds a connection to the Helmholtz free
energy through the limit
$F(\beta,v)=-\beta^{-1}\lim_{N\to\infty}\ln{Z(\beta,V,N)}$, while $v=V/N$
tends to a constant. Here, $\beta=1/k_{B}T$ represents the canonical inverse
temperature, $k_{B}$ signifies the Boltzmann constant, and $v$ denotes the
specific volume. Consequently, the non-analytic behavior of $F(\beta,v)$
occurs at points where $\lim_{N\to\infty}Z(\beta,V,N)=0$.

The groundbreaking work of Lee and Yang [2] and its extension by Fisher [3]
established that studying the zeros of the partition function offers a
rigorous framework for understanding phase transitions [4, 5, 6]. Although the
partition function is composed exclusively of positive terms, implying the
absence of real positive roots for any finite system, inspecting the zeros of
small systems allows for the revelation of some properties of the
thermodynamic system, for instance, the transition temperature.

Consider the analytical continuation of the partition function
$Z=Z(\mathcal{B})$ with $\mathcal{B}=\beta+i\tau$, where $\beta$ represents
the inverse temperature and $\tau$ is an imaginary parameter. In the
thermodynamic limit, a phase transition occurs at
$\lim_{N\rightarrow\infty}Z(\mathcal{B}_{k})=0$ if $\tau_{k}=0$. The manner in
which the zeros approach this limiting point is indicative of the transition’s
order.

In the late 1960s, S. Grossmann and W. Rosenhauer [7, 8] showed that the
phenomenologically known types of phase transition can be characterized by the
way that the density of zeros, which is the thermodynamic limit of the
distribution of zeros (DOZ), behave toward the transition point. They proposed
a general Finite-Size Scale (FSS) method for the DOZ which accumulates in
lines that tend to cut the real axis under a certain slope,
$\gamma=(\beta-\beta_{c})/\tau$, whereas the density function can be described
by a simple power law $\phi(\tau)\approx\tau^{\alpha}$. After that, S.
Grossmann and V. Lehmann [9] provided some results of this method for
realistic physical models.

Towards the end of the twentieth century, P. Borrmann et al [10] proposed a
classification scheme for phase transitions in finite systems based on the
method introduced by S. Grossmann and W. Rosenhauer. Similarly, by analyzing
the DOZ they classified the type of the transition by both: the angle of the
zeros lines toward the real axis and the distance between the zeros in this
line. For a pseudo-first-order phase transition, this line is perpendicular to
the real axis and, concomitantly, the zeros are evenly spaced, see Fig. 3.

More recently, M.P.Taylor et al[11] empirically demonstrated the connection
between the curvature properties of entropy, denoted as $S$, and the DOZ. In
the microcanonical analysis, a convex behavior of $S$, i.e. an unstable
region, is related to a first-order transition [12]. The double-touching
tangent line construction on this convex intruder can define both the energy
range of the non-stable region and the transition temperature. The authors
calculated the zeros of $Z$ by considering this truncated energy range and
$x=e^{-\mathcal{B}E}$ as a variable and showed that it leads to a circle on
the complex $x$ plane map. Solving it for $\mathcal{B}$, this circle leads to
a vertical line on the complex $\mathcal{B}$ plane map, which corroborates
with P. Bormann and collaborators’ results. Additionally, they observed
another pattern of zeros that pinch the real axis, which they attributed to a
higher-order transition.

In the present work, we propose an alternative analytical argument to
establish a connection between the unstable region of the entropy and the
vertical line pattern observed in the Distribution of Zeros (DOZ), as
empirically shown by M.P. Taylor and collaborators. This paper is organized as
follows: in section 2.1, we present the microcanonical analysis of phase
transition. Subsequently, in section 2.2 we introduce the Fisher zeros and in
the section 2.3 the classification scheme proposed by P. Borrmann et al. In
section 3.1 we outline the analytical arguments that a first-order transition
leads to a vertical line pattern of the zeros on the complex $\mathcal{B}$
plane map. Our results are then compared with a Monte Carlo simulation of the
$147$ particles Lennard-Jones cluster in section 3.2. Finally, in section 4 we
present our conclusions.

##  2 Methodology

###  2.1 Microcanonical Ensemble

In the microcanonical approach to statistical mechanics, entropy carries all
the information necessary to describe the system. The first probabilistic
statement for entropy was made for the ideal gas in 1872 by L. Boltzmann [13].
In 1901, M. Planck stated his famous formula,

$S(E)=k_{B}\ln{\Omega(E)},$ (1)

as the expression for the entropy of black bodies [14], with $\Omega(E)$
standing for the number of ways in which the system can be can be realized
with energy $E$. For simplicity, in this work we measure $S$ in units of
$k_{B}$.

In microcanonical statistics, the equilibrium state of a thermodynamic system
is defined by derivatives of $S$. For instance, the inverse microcanonical
temperature is given by the derivative of entropy $S$ with respect to energy
$E$, while keeping the set of independent extensive quantities $\\{X\\}$ (such
as volume $V$, number of particles $N$, magnetization $M$, etc.) constant,
i.e.

$\bar{\beta}(E)=\bar{T}^{-1}=\left(\frac{\partial S}{\partial
E}\right)_{\\{X\\}}.$ (2)

The overbar is used to emphasize that the quantity is a microcanonical
parameter. It is important to note that $\bar{\beta}/k_{B}$ recovers the usual
canonical $\beta$ in the thermodynamic limit.

Considering an energy region without any transition, the function $S(E)$ is
strictly monotonically increasing, concave, and positive. Consequently,
$\bar{\beta}$ is a monotonically decreasing, convex, and positive function.
Higher-order derivatives of entropy, denoted by

$\displaystyle\bar{\gamma}(E)=\left(\frac{\partial^{2}S}{\partial
E^{2}}\right)_{\\{X\\}}$ and
$\displaystyle\bar{\delta}(E)=\left(\frac{\partial^{3}S}{\partial
E^{3}}\right)_{\\{X\\}},$ (3)

are increasing concave negative and positive decreasing convex positive
functions, respectively, and so on.

A convex behavior of the entropy indicates a non-stable region, so that, a
change in the concavity of $S(E)$ corresponds to a first-order phase
transition. The touching points of the double-tangent line across the convex
region define the latent heat and the energy range of the transition,
$[E^{\prime},E^{\prime\prime}]$. Additionally, the slope of this line defines
the transition temperature, see Fig. 2.

###  2.2 Fisher’s Zeros

From the point of view of the Canonical Ensemble, the partition function can
be seen as the Laplace transform of $\Omega(E)$. For a system with a
continuous energy domain, it can be written as:

$Z(\mathcal{B},V,N)=\int\mathrm{d}E\ \Omega(E)e^{-\mathcal{B}E},$ (4)

where, $\mathcal{B}=\beta+i\tau$ represents the complex inverse temperature
and $E$ depends on $V$ and $N$. One can introduce a discretization approach by
considering an energy bin of size $\varepsilon$. Consequently, the energy of
the $k$-th bin can be expressed as:

$E_{k}=E_{0}+k\varepsilon,$ (5)

where $E_{0}$ stands for the ground state energy. In this approach,
$\Omega(E_{k})$ represents the number of states with energy between $E_{k}$
and $E_{k}+\varepsilon$. By considering a discrete version of eq. (4) and
incorporating the energy given by eq. (5), we can express the partition
function as:

$Z(\mathcal{B},V,N)=e^{-\mathcal{B}E_{0}}\sum_{k=0}^{n}\Omega_{k}e^{-\mathcal{B}k\varepsilon},$
(6)

where $\Omega_{k}\equiv\Omega(E_{k})$ and $n$ is the number of energy bins.
Following Fisher we define a new variable

$x\equiv
e^{-\varepsilon\mathcal{B}}=e^{-\varepsilon\beta}e^{-i\varepsilon\tau},$ (7)

which allows us to rewrite the partition function as a polynomial:

$Z(\mathcal{B},V,N)=e^{-\mathcal{B}E_{0}}\sum_{k=0}^{n}\Omega_{k}x^{k}=e^{-\mathcal{B}E_{0}}\prod_{k=1}^{n}\left(x-x_{k}\right),$
(8)

where $\Omega_{k}^{\prime}s$ are identified as the coefficients of the
polynomial and $x_{k}$ represents the $k$-th zero. It is worth mentioning that
a multiplicative constant in the polynomial does not alter its roots.
Consequently, instead of dealing directly with the number of states, sometimes
it can be preferable to work with the density of states (DOS), defined as
$g(E)=\Omega(E)/\sum_{E}\Omega(E)$.

According to the fundamental theorem of algebra, an $n$-th-order polynomial
has exactly $n$ zeros, including multiplicities. Besides that, the roots of
the polynomial occur in complex conjugated pairs
($x_{k_{\pm}}=e^{-\varepsilon\beta_{k}}e^{\pm i\varepsilon\tau_{k}}$). Since
all coefficients in the polynomial are real and positive, any real zeros must
be negative, at least for a finite-order polynomial. If $Z$ possesses real
positive roots, the corresponding $F$ becomes singular at those points,
indicating the presence of phase transitions in the system. Implying that a
real positive zero is only possible at the thermodynamic limit.

All thermodynamic functions can be derived from the zeros, for instance, the
specific heat at constant volume,

$\displaystyle c_{V}$ $\displaystyle=$
$\displaystyle\frac{k_{B}\beta^{2}}{N}\left(\frac{\partial^{2}\ln{Z}}{\partial\beta^{2}}\right)$
(9) $\displaystyle=$
$\displaystyle\frac{k_{B}x(\ln{|x|})^{2}}{N}\sum_{k=1}^{n}\left(\frac{-x_{k}}{(x-x_{k})^{2}}\right).$

In this work, $c_{V}$ is measured in units of $k_{B}$. It is observed that a
singular behavior of the specific heat may emerge for $x=x_{k}$ and
$\tau_{k}\ll 1$. While true phase transitions are not possible in finite
systems, it is expected that a particular zero, known as the dominant or
leading zero, approaches the real positive axis, indicating a pseudo-phase
transition.

###  2.3 Classification of the Order of the Phase Transition

P. Borrmann et al [10] proposed a discretized version of the phase transition
classification scheme of S. Grossmann and W. Rosenhauer [7, 8]. In this
section, we provide a brief outline of their main results. They considered the
zeros that are close to the real axis to lie approximately on a straight line
making an angle $\delta=\arctan{(\gamma)}$ with the imaginary axis, where

$\gamma=\frac{\beta_{2}-\beta_{1}}{\tau_{2}-\tau_{1}},$ (10)

as show in Fig. 1. It is worth mentioning that the indexes start from $1$ and
increase with $\tau$, the zero labeled $1$ is the leading zero. The crossing
point of the line with the real axis is
$\beta_{cut}=\beta_{1}-\gamma\tau_{1}$.

Figure 1: (Color online) Reproduction of the scheme of the DOZ toward the real
axis from P. Borrmann et al [10]

A discrete density of zeros, $\phi(\tau_{k})$, is defined as the average of
the distances between the first near zeros as

$\phi(\tau_{k})=\frac{1}{2}\left(\frac{1}{\|\mathcal{B}_{k}-\mathcal{B}_{k-1}\|}+\frac{1}{\|\mathcal{B}_{k+1}-\mathcal{B}_{k}\|}\right),$
(11)

with $k=2,3,4\cdots$. Since zeros with small imaginary parts contribute more
to the specific heat at the transition (or any other thermodynamic functions
that is singular at this point) they supposed that $\phi$ can be approximated
by a simple power law, i.e. $\phi(\tau)\sim\tau^{\alpha}$. An estimate of the
exponent $\alpha$ can be done using two zeros as

$\alpha=\frac{\ln{\phi(\tau_{3})}-\ln{\phi(\tau_{2})}}{\ln{\tau_{3}}-\ln{\tau_{2}}}.$
(12)

A first-order phase transition is defined by $\alpha=0$ and $\gamma=0$, i.e. a
vertical line of evenly spaced zeros.

##  3 Results

###  3.1 Fisher’s zeros for a first-order phase transition

In this section, we present an alternative demonstration that for a pseudo-
first-order transition, the zeros of the partition function exhibit a vertical
line pattern in the complex inverse temperature plane. We divide the domain of
the partition function, given by eq. (6), into three parts:
$Z(\mathcal{B},V,N)=Z_{<}+Z^{\prime}+Z_{>}$. The first part, $Z_{<}$, includes
energies $E<E^{\prime}$, $Z^{\prime}$ covers the energy range of the non-
stable region [$E^{\prime}$, $E^{\prime\prime}$] (as discussed in section
2.1), and $Z{>}$ accounts for energies $E>E^{\prime\prime}$.

It can be claimed that $Z^{\prime}(\mathcal{B}=\mathcal{B}_{j})\approx 0$
since approaches that truncate the energy range, such as the zeros of the
density of states [15, 16, 17], can capture indications of phase transitions.
Thus, we have:

$Z^{\prime}(\mathcal{B}_{j},V,N)=\sum_{E=E^{\prime}}^{E^{\prime\prime}}\Omega(E)e^{-\mathcal{B}_{j}E}\approx
0.$ (13)

To deal with the convexity of the entropy, we expand $S$ in a Taylor series
around the midpoint $E_{in}=(E^{\prime}+E^{\prime\prime})/2$ and collect terms
up to the third order:

$\displaystyle S(E)$ $\displaystyle\approx$ $\displaystyle
S_{in}+\bar{\beta}_{in}(E-E_{in})+\frac{\bar{\gamma}_{in}}{2}(E-E_{in})^{2}$
(14) $\displaystyle+$
$\displaystyle\frac{\bar{\delta}_{in}}{6}(E-E_{in})^{3},$

where $S_{in}=S(E_{in})$ and $\bar{\beta}_{in}=\bar{\beta}(E_{in})$,
$\bar{\gamma}_{in}=\bar{\gamma}(E_{in})$, and
$\bar{\delta}_{in}=\bar{\delta}(E_{in})$ are the derivatives of $S$ as defined
in section 2.1. In the considered energy range, $E=E^{\prime}+k\varepsilon$.
Defining $\Delta E=E^{\prime\prime}-E^{\prime}$ so that
$E^{\prime}=E_{in}-\Delta E/2$, we can write:

$S(E)\approx S^{\prime}+\bar{\beta^{\prime}}\varepsilon
k+\frac{\bar{\gamma^{\prime}}}{2}\varepsilon^{2}k^{2}+\frac{\bar{\delta}_{in}}{6}\varepsilon^{3}k^{3},$
(15)

where,

$S^{\prime}=S_{in}-\frac{\bar{\beta}_{in}}{2}\Delta
E+\frac{\bar{\gamma}_{in}}{8}\Delta E^{2}-\frac{\bar{\delta}_{in}}{48}\Delta
E^{3},$ (16)
$\bar{\beta^{\prime}}=\bar{\beta}_{in}-\frac{\bar{\gamma}_{in}}{2}\Delta
E+\frac{\bar{\delta}_{in}}{8}\Delta E^{2}=-\frac{\partial S^{\prime}}{\partial
E^{\prime}},$ (17)

and

$\bar{\gamma^{\prime}}=\bar{\gamma}_{in}-\frac{\bar{\delta}_{in}}{2}\Delta
E=-\frac{\partial\bar{\beta^{\prime}}}{\partial
E^{\prime}}=\frac{\partial^{2}S^{\prime}}{\partial E^{\prime 2}}.$ (18)

Inserting eq. (15) into eq. (1) and solving for $\Omega(E)$, eq. (13) can be
rewritten as:

$Z^{\prime}_{n}(\mathcal{B}_{j})\approx
e^{-\mathcal{B}_{j}F^{\prime}}\sum_{k=0}^{n^{\prime}}x^{k}y^{k^{2}}z^{k^{3}},$

where $n^{\prime}$ is the number of energy levels in the energy range of the
non-stable region, $F^{\prime}=E^{\prime}-S^{\prime}/(k_{B}\mathcal{B}_{j})$,

$\displaystyle x$ $\displaystyle=$
$\displaystyle\exp{\left[-\left(\mathcal{B}_{j}-\frac{\bar{\beta^{\prime}}}{k_{B}}\right)\varepsilon\right]}$
(19) $\displaystyle=$
$\displaystyle\exp{\left[-\left(\beta_{j}-\frac{\bar{\beta^{\prime}}}{k_{B}}\right)\varepsilon\right]}\exp{\Big{[}-i\tau_{j}\varepsilon\Big{]}},$
$y=\exp{\left(\frac{\bar{\gamma^{\prime}}}{2k_{B}}\varepsilon^{2}\right)},$

and

$z=\exp{\left(\frac{\bar{\delta}_{in}}{6k_{B}}\varepsilon^{3}\right)}.$

Usually, $\varepsilon$, $\bar{\gamma}_{in}$ and $\bar{\delta}_{in}$ are small
quantities, so $y\approx z\approx 1$ giving:

$Z^{\prime}\approx
e^{-\mathcal{B}_{j}F^{\prime}}\sum_{k=0}^{n^{\prime}}x^{k}=e^{-\mathcal{B}_{j}F^{\prime}}\frac{1-x^{n^{\prime}+1}}{1-x}.$
(20)

By collecting terms up to first order, i.e. considering a linear behavior of
the entropy, it leads to the same relation for $Z^{\prime}$. Hence, one can
say that the double-touching tangent line construction is a good approach even
for finite systems. By inspecting eqs. (20) and (19), we get $Z^{\prime}=0$ if

$\beta_{j}=\frac{\bar{\beta^{\prime}}}{k_{B}},$ (21)

and

$\tau_{j}=\frac{2\pi j}{\varepsilon(n^{\prime}+1)}=\frac{2\pi}{\Delta E}\ j$
(22)

where $j=1,2,\cdots,n^{\prime}$. It is worth mentioning that $j\neq
0,(n^{\prime}+1)$, since the denominator in the last term of eq. (20) requires
that $x\neq 1$, hence $\mathcal{B}_{j}$ can not be a positive real number.
Furthermore, any other $j$ will lead to multiplicities and can be neglected.
Since $\bar{\beta^{\prime}}$ is a constant, given by eq. (17), plotting the
ordered pairs ($\beta_{j}$, $\tau_{j}$) leads to a vertical line of evenly
spaced points as claimed before. Besides that, by inserting eq. (17) into eq.
(21), we obtain:

$k_{B}\beta_{j}=\bar{\beta}_{in}-\frac{\bar{\gamma}_{in}}{2}\Delta
E+\frac{\bar{\delta}_{in}}{8}\Delta E^{2}.$ (23)

###  3.2 Zeros Map for the Lennard-Jones Cluster

In this section, we illustrate the discussion above by using the example of
the Lennard-Jones (LJ) cluster with $N=147$ particles, which is considered a
prototype of a pseudo-first-order phase transition. The LJ cluster consists of
particles interacting through the pairwise LJ potential:

$U_{LJ}(r_{ij})=4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right],$
(24)

where $r_{ij}=|\mathbf{r}_{j}-\mathbf{r}_{i}|$ is the distance between
particles identified by the indices $i$ and $j$, and $\mathbf{r}_{i}$ and
$\mathbf{r}_{j}$ denote the respective positions of these particles. The
reduced parameters are set such that the minimum of the potential is at
$r_{ij}=r_{0}=1$, and energy is measured in units of $\epsilon$
($\sigma=2^{-1/6}$ and $\epsilon=1$). The particles are confined to a sphere
of radius $r_{c}=4\sigma$ to reproduce the transition temperature
($T_{tr}\approx 0.36$) presented by P.A. Frantsuzov and V.A. Mandelshtam[18].
The results presented in this section are averages of five independent
simulations, and errors are given by the standard deviation, except for Fig.
3, where the zeros map of each individual simulation is shown. The specific
details of the simulations can be found in Appendix A.

Fig. 2 shows the specific entropy, $s=S/N$, as a function of the energy
density, $e=E/N$, for the $147$-LJ cluster. One can observe the convex
intruder inside the dotted green rectangle, which is zoomed in the inset. The
blue dashed line is the double-touching tangent line construction, which leads
to a slope $\bar{\beta}_{tan}=2.751(9)$, and the energy density range of the
unstable region is $[e^{\prime}=-5.2286(9),e^{\prime\prime}=-4.861(1)]$. The
specific latent heat then is calculated to be $q_{L}=2.78(1)$.

Figure 2: (Color online) Estimation of the specific entropy for the $147$
particles Lennard-Jones Cluster. The error bars are in the same order as the
line width. The dotted green rectangle demarcates the unstable region. The
inset is a zoom in this region where the convex intruder can be perceived. The
dashed blue line is the double-touching tangent line construction. The small
dashed purple vertical line marks $e_{in}=(e^{\prime\prime}+e^{\prime})/2$.

Fig. 3 displays the Fisher zeros distribution map for the $147$-LJ cluster.
Each symbol indicates the result of an independent simulation. The leading
zero is found at ($\beta_{1}=2.761(2)$, $\tau_{1}=0.0609(6)$), our maps also
show that although the zeros are sensitive to statistical fluctuations, the
zeros in the transition region are quite stable [19]. In Fig. 4 we show an
adaptation of the scaling analysis proposed by Borrmann et al, discussed in
section 2.3. We propose a linear fit in
$\ln{(1/\|\mathcal{B}_{k}-\mathcal{B}_{k-1}\|)}\times\ln{(\tau_{k})}$, for
$k=2,3,4,$ and $5$. We found the coeficient $\alpha=0.058(7)$, which is
coherent to the approach proposed by eq. (12), $\alpha=0.041(5)$. In the inset
of this figure we show the linear fit of the dominant zeros where we found the
slope $\gamma=-0.004(3)$ which leads to an angle $\delta=0.2(2)^{\circ}$, and
the crossing point $\beta_{cut}=2.7601(9)$. Those parameters are consistent
with the first-order phase transition. Besides that, they are also consistent
with the approach values proposed by Borrmann et all, $\gamma=-0.021(1)$, and
$\beta_{cut}=2.762(2)$. The average of the distances between the dominant
zeros is $0.110(2)$. From eq. (22) one can see that this distance is
$\Delta\tau=2\pi/[N(e^{\prime\prime}-e^{\prime})]=0.1162(4)$, corroborating
for the validity of the demonstration.

Figure 3: (Color online) The Fisher zeros distribution map for the $147$
particles Lennard-Jones cluster. Each symbol indicates the results of an
independent simulation. $\delta$ is the angle between the fitted line of the
dominant zeros with the vertical axis. See the inset of Fig. 4 for a zoom in
this region.  Figure 4: (Color online) $\log\times\log$ graph of the inverse
of the absolute value of the difference between the complex inverse
temperature of adjacent dominant zeros versus the complex part of the inverse
temperature, i.e.
$-\ln{\|\mathcal{B}_{k}-\mathcal{B}_{k-1}\|}\times\ln{(\tau_{k})}$, for
$k=2,3,4$, and $5$. In the inset we show the real part versus the imaginary
part of the dominant zeros.

The reliability of the zeros maps and their relationship with other quantities
is further discussed. We have chosen the MPSolve [20, 21] routine as the zeros
finder for this study. Besides the roots of polynomials, this routine’s output
can also return error bars. In this examination, the error bars are the order
of $10^{-12}$. Upheld by obtaining $\sum_{i}\tau_{i}\approx 0$, since the
zeros come in complex conjugated pairs, we can endorse the precision of the
routine in this case. To prove accuracy, one can calculate a given
thermodynamic function by the Fisher’s zeros and compare it with one obtained
via DOS. As a check, we compare the specific heat at constant volume obtained
by eq. (9) and by the standard canonical average:

$c_{V}=\frac{k_{B}\beta^{2}}{N}\left(\left\langle
E^{2}\right\rangle-\left\langle E\right\rangle^{2}\right),$ (25)

where

$\left\langle E^{k}\right\rangle=\sum_{E}E^{k}P(E,\beta),$ (26)

and

$P(E,\beta)=\frac{g(E)e^{-\beta E}}{Z},$ (27)

is the Boltzmann probability density. We then define the relative difference,

$\Delta c_{V}=\left\|1-\frac{c_{V}(z)}{c_{V}(g)}\right\|,$ (28)

where $c_{V}(g)$ is obtained from the DOS and $c_{V}(z)$ is obtained from the
zeros, as comparative metric. This inspection is shown in Fig. 5, where we can
state that the numerical imprecision provided by the zeros finder is
negligible in this case. Thus, we have high confidence in the legitimacy of
the zeros map. In addition, one can recognize that the $\beta_{1}$, indicated
by the dotted-dashed green line, is close to the temperature of the peak
position of the $c_{V}$.

Figure 5: (Color online) Specific heat at constant volume for the $147$
Lennard-Jones Cluster ($V=4^{4}\pi\sigma^{3}/3$). The black circles stands for
$c_{V}$ evaluated via the DOS, eq. (25). The red square stands for $c_{V}$
calculated via the Fisher zeros, eq. (9). The inset shows the relative
difference between the two values, see eq. (28). The dotted-dashed green line
indicates $\beta_{1}$ from the zeros maps, the dashed blue line indicates
$\bar{\beta}_{tan}$ from the tangent line of double-touching tangent line
construction, and the small dashed purple line indicates $\bar{\beta}_{in}$
and the double-dotted-dashed magenta line indicates $\bar{\beta}_{tr}$, which
will be discussed later.

Due to the coexistence of phases, the Boltzmann probability density presents
two peaks in a first-order transition, each related to a phase. At the
transition temperature, one expected that those peaks have the same height.
Since one can rewrite eq. (27) as $P(E,\beta)=\exp{(-\beta F)}/Z$, this
analysis is similar to the minimum condition of the Helmholtz free energy. In
Fig. 6 we show the Boltzmann probability density for four temperatures:
$\beta_{1}$ and $\bar{\beta}_{tan}$, already estabilished, and
$\bar{\beta}_{in}$ and $\bar{\beta}_{tr}$, discussed in the next paragraph.
One can see that the Fisher zeros analysis is coherent with the equal
probability condition, and the double-touching tangent construction slightly
deviates from it.

Figure 6: (Color online) The Boltzmann probability density of the $147$
Lennard-Jones Cluster. The dotted-double-dashed red vertical line marks the
microcanonical transition point. The unstable region is demarcated by the
dotted green line.

Finally, we show the microcanonical analysis of least-sensitive inflection
points for the $147$-LJ cluster. The change in the curvature of $S(E)$ causes
an inflection point, called the inflection point of least sensitivity if the
derivative changes least on variation in energy and provides a signal of the
transition at this energy, $E_{tr}$ [12]. Let
$\bar{\beta}_{tr}=\bar{\beta}(E_{tr})$,
$\bar{\gamma}_{tr}=\bar{\gamma}(E_{tr})$, and
$\bar{\delta}_{tr}=\bar{\delta}(E_{tr})$ the higher order derivatives of $S$
evaluated in $E_{tr}$. According to the least-sensitive inflection point
microcanonical analysis, for a pseudo-first-order transition,
$\bar{\gamma}_{tr}$ is a maximum positive value. In Fig. 7 we show the
microcanonical inverse temperature, $\bar{\beta}$, just for the unstable
region, i.e. the derivative of the entropy shown in the inset of Fig. 2. In
conformity, the dashed blue line is the derivative of the double-touching
tangent line construction. For comparison purposes, we show $k_{B}\beta_{1}$
in the dotted-dashed green line. We also show
$\bar{\beta}_{in}=\bar{\beta}(e_{in})$ in the small dashed purple line.
$\bar{\beta}_{in}>k_{B}\beta_{1}$ as predicted by eq. (23). $k_{B}\beta_{1}$
line is in accordance with the Maxwell’s equal area construction, since
$A_{1}\approx A_{2}$. In the inset, we show $\bar{\gamma}$, measured in units
of $k_{B}/\epsilon^{3}$, where the peak position defines the microcanonical
transition point, $e_{tr}$. The double-dotted-dashed magenta line indicates
the microcanonical transition temperature, i.e.
$\bar{\beta}_{tr}=\bar{\beta}(e_{tr})$.

Figure 7: (Color online) The microcanonical inverse temperature in the
unstable region. The dashed blue line indicates $\bar{\beta}_{tan}$ from the
double-touching tangent line construction, the dotted-dashed green line
indicates $k_{B}\beta_{1}$ from the zeros maps, the small dashed purple line
indicates $\bar{\beta}_{in}=\bar{\beta}(e_{in})$, and the double-dotted-dashed
magenta line indicates $\bar{\beta}_{tr}=\bar{\beta}(e_{tr})$ from the
microcanonical analysis. The error bars are the same size as the symbols. The
hued areas $\mathrm{A}_{1}$ and $\mathrm{A}_{2}$ are consistent with the
Maxwell’s equal area construction. The inset shows $\bar{\gamma}$. The dotted-
double-dashed red line marks the peak position of $\bar{\gamma}$, i.e. the
microcanonical transition point.

It is worth mentioning that, although the Fisher zeros analysis corroborates
with the Maxwell’s equal area construction, equal probability condition, and
provides a transition temperature close to the temperature of the peak
position of the specific heat. It is well known that, for finite systems,
different quantities provides different transitions temperatures [19],
converging to the transition value as the thermodynamic limit is approached.
Thus, this specific study is inconclusive about the accuracy of distinct
methods, a statement in this regard requires extensive work, and this is not
the purpose of this manuscript.

##  4 Conclusion

In this work, a mathematical argument to connect the nonstable region of the
entropy to the pattern of the Fisher’s zeros map was presented. The first term
of Taylor’s series approach of the entropy leads to the vertically lined
equally spaced zeros of the partition function on the complex inverse
temperature plane for the first-order phase transition. This behavior was
illustrated via the Lennard-Jones cluster. For this specific model, the
transition temperature defined by the leading zero corroborates with the peak
position of the specific heat, also with the Maxwell’s equal area
construction, and with the equal probability condition of phases at the first-
order transition.

## Acknowledgments

We would like to acknowledge helpful conversations with Dr. Michael Bachmann.
This work received public financial support from Fundação de Amparo à Pesquisa
do Estado de Minas Gerais (FAPEMIG), Brazil, under the grant RED-00458-16.

## Declarations

The authors have no competing interests to declare that are relevant to the
content of this article.

##  Appendix A Details of the Simulations

In this appendix, we present a detailed description of the Monte Carlo
simulation method used for studying the Lennard-Jones cluster. The Monte Carlo
method is a class of statistical algorithms that sample a limited but
representative number of states to infer some properties of the system under
study. Those states can be chosen in a Markov chain, i.e. the probability to
sample each state depends only on the previous state. Mathematically, this
condition can be stated by the detailed balance,

$P_{i}W_{i\to j}=P_{j}W_{j\to i},$ (29)

where $W_{i\to j}$ is the transition probability from state $i$ to state $j$,
and $P_{i}$ is the equilibrium probability of being in state $i$ [22]. The
Metropolis prescription to satisfy this condition is

$\text{W}_{i\to j}=\text{min}\left\\{1,\frac{P_{j}}{P_{i}}\right\\}.$ (30)

We want a Monte Carlo scheme to estimate the entropy, it can be done by one of
the flat histogram methods, here we choose the Wang-Landau Sampling [23]. To
understand this method, let us look at the Boltzmann distribution for
$\beta=0$. In this situation, eq. (27) can be written as $P(E)=\Omega(E)/Z$.
So, the probability of randomly tossing a state with energy $E_{i}$ is
proportional to $\Omega(E_{i})$. If we accept the tossed state to our sampled
set with probability $P_{i}=1/\Omega(E_{i})$, all energies will be equally
sampled. Of course, we are unaware of $\Omega(E)$, but we can use this equally
sampled energies fact to estimate it as follows: We create a histogram to
count how many states with a given energy were sampled, $h(E)$. Since
$\Omega(E)$ can assume very large numbers, let us work with the entropy. So we
first guess an initial value to $S(E)/k_{B}$, for instance,
$\ln(\Omega(E))=1$, and define an initial state, $i$. Hereinafter, we randomly
guess a new state, $j$, and compare the states $i$ and $j$ by the Metropolis
prescription. Considering the proposed probability, it can be written as

$\text{W}_{i\to
j}=\text{min}\left\\{1,\frac{\Omega(E_{i})}{\Omega(E_{j})}\right\\}.$ (31)

If the trial state is accepted we set $j$ as the current one, $i\leftarrow j$.
At every trial move $\Omega(E_{i})$ is updated by a multiplicative factor $f$,
i.e., $\ln{\Omega(E_{i})}\leftarrow\ln{\Omega(E_{i})}+\ln{(f)}$.
Simultaneously, the histogram is also updated, $h(E_{i})\leftarrow
h(E_{i})+1$. When $h(E)$ is flat we can say that we approach $\ln{\Omega(E)}$
with precision equal to $\ln(f)$. We considered the energy ranging from
$0.95E_{\text{min}}$ to $E_{\text{max}}=0$. Where $E_{\text{min}}$ is the
ground state given by J.A. Northby [24]. We also consider one trial move the
attempt to change the position of a single particle. The new position is
chosen inside a small sphere of radius $r_{t}$ centered in the original
position of the particle. The value of $r_{t}$ is chosen so that the
acceptance ratio is close to $60\%$. To quickly sample the entire
configuration space, a large initial value for $f=f_{0}$ is required, the
original recommendation states that $\ln{(f_{0})}=1$. To improve the precision
$f$ must be decreased and the scheme repeated.

The histogram flatness is tested after $10^{6}$ Monte Carlo sweeps (MCS). One
MCS is counted after a sequential attempt to change all particles of the
system once. If the histogram is flat, it is reset, $h(E)=0$ and $f$ is
decreased. The histogram is considered flat when the ratio of its lowest value
by the mean value is greater than $p$, in this work $p=0.70$. At first, any
function can be used to decrease $f$, we also used the original suggestion,
i.e. $\ln{(f_{i+1})}=\ln{(f_{i})}/2$. The scheme is repeated until the desired
precision is reached, in this work we cease the process when
$\ln{(f)}=\varepsilon=10^{-9}$. Regrettably, the capacity to diminish the
inaccuracy of $\Omega(E)$ asymptotically halts as the modification factor $f$
decreases. This phenomenon is recognized as the saturation of the error
between the calculated and the exact $\Omega(E)$, a concept first elucidated
in Ref. [25]. Various improvements to the WL-method have been proposed, for
instance, the $1/t$-Wang-Landau approach [26, 27, 28] and the optimal
modification factor [29]. Given the analogous behavior observed in the
convergence patterns of simple sampling Monte Carlo and $1/t$-Wang-Landau
[30], we compare our results with those obtained by the regular Metropolis
algorithm [31], as illustrated in Fig. 8. The discussion of this comparison is
presented in the final paragraph of this appendix.

Moreover, the standard WL method is very time-consuming, so we opted for a
parallelization procedure, called Replica Exchange Wang-Landau (REWL) method
[32]. The idea is to divide the energy range into several smaller pieces,
called windows. In this work, we divided the energy range in $N_{W}=10$ equal
size windows of $10^{4}$ energy bins. One or more WL samples, called walkers,
are performed in parallel at each window. Here we consider $N_{R}=4$ walkers
in each window. In addition, an attempt to exchange configurations of walkers
between adjacent windows is proposed after $10^{3}$ MCS. An exchange between
conformations $X$ and $Y$, respectively located at neighboring windows $i$ and
$j$, is proposed with the probability

$\text{P}_{\text{acc}}=\text{min}\left\\{\frac{\Omega_{i}(E[X])}{\Omega_{i}(E[Y])}\frac{\Omega_{j}(E[Y])}{\Omega_{j}(E[X])},1\right\\}.$
(32)

This exchange allows the walkers to efficiently sample different parts of the
configuration space, this procedure is as crucial as dividing the windows to
improve the simulation time. The acceptance ratio of the replica exchange is
tied to the overlap between the windows, in this work we set an overlap of
$75\%$. When the final precision is reached, the pieces are combined to form
the entire entropy. We concatenate the pieces at the point of the smallest
difference of the inverse temperature between the adjacent windows. There are
$N_{R}^{N_{W}}$ possible combinations of concatenations of the pieces of the
entropy, we randomly chose $10^{3}$ of them and the final result is the
average value of those combinations via Jackknife resampling. The average
procedure presented here holds significant importance as it effectively
mitigates the potential introduction of spurious non-analytical points in the
entropy, arising from the interconnection of distinct entropy fragments. These
artifacts, if not properly handled, might erroneously suggest the presence of
non-existent phase transitions. Since each pair of walkers is connected at a
different point, the implementation of the Jackknife resampling technique
induces a larger weight in the continuous region of the pieces of the entropy,
which ultimately improves the accuracy of the mean value.

On the question of the convergence issue of the Wang-Landau method we check
the Boltzmann distribution obtained by the REWL method with the one obtained
by the regular Metropolis Algorithm [31], see Fig. 8. We calculate the
$P(E,\beta)$ for two temperatures, one above the transition temperature
($\beta=2$) and another below ($\beta=3$). Those temperatures are far away
from the transition to avoid the Metropolis algorithm being stuck in a meta-
stable state [33]. For the Metropolis Algorithm, we excluded the first
$10^{5}$ MCS for thermalization purposes and, after that, performed $10^{7}$
MCS to obtain $P(E)$. The result presented here is also an average of 5
independent simulations. Besides that, the trial move is similar to that used
for the WL method. The relative differences between the two methods are of the
order of the error bars, see the inset in Fig. 8, demonstrating the
reliability of the REWL procedure.

Figure 8: (Color online) Comparison of the Boltzmann probability density of
the $147$ Lennard-Jones Cluster for $\beta=2.0$ and $\beta=3.0$ obtained by
the Replica-Exchange-Wang-Landau method and by the Metropolis algorithm.

## References

  * [1] Tilman Sauer.  Statistical theory of equations of state and phase transitions. ii. lattice gas and ising model.  The European Physical Journal Special Topics, 226(4):539–549, 2017.
  * [2] T. D. Lee and C. N. Yang.  Statistical theory of equations of state and phase transitions. ii. lattice gas and ising model.  Phys. Rev., 87:410–419, Aug 1952.
  * [3] M.E. Fisher.  The nature of critical points.  In W. E. Brittin, editor, Lectures in Theoretical Physics, Volume VII C - Statistical Physics, Weak Interactions, Field Theory, Lectures Delivered at the Summer Institute for Theoretical Physics. University of Colorado Press, Boulder, 1965.
  * [4] David Ruelle.  Some remarks on the location of zeroes of the partition function for lattice systems.  Communications in Mathematical Physics, 31(4):265–277, 1973.
  * [5] Bo-Bo Wei and Ren-Bao Liu.  Lee-yang zeros and critical times in decoherence of a probe spin coupled to a bath.  Phys. Rev. Lett., 109:185701, Oct 2012.
  * [6] Xinhua Peng, Hui Zhou, Bo-Bo Wei, Jiangyu Cui, Jiangfeng Du, and Ren-Bao Liu.  Experimental observation of lee-yang zeros.  Phys. Rev. Lett., 114:010601, Jan 2015.
  * [7] S. Grossmann and W. Rosenhauer.  Temperature dependence near phase transitions in classical and quant. mech. canonical statistics.  Zeitschrift für Physik, 207(2):138–152, 1967.
  * [8] S. Grossmann and W. Rosenhauer.  Phase transitions and the distribution of temperature zeros of the partition function.  Zeitschrift für Physik A Hadrons and nuclei, 218(5):437–448, 1969.
  * [9] S. Grossmann and V. Lehmann.  Phase transitions and the distribution of temperature zeros of the partition function.  Zeitschrift für Physik A Hadrons and nuclei, 218(5):449–459, 1969.
  * [10] Peter Borrmann, Oliver Mülken, and Jens Harting.  Classification of phase transitions in small systems.  Phys. Rev. Lett., 84:3511–3514, Apr 2000.
  * [11] Mark P. Taylor, Pyie Phyo Aung, and Wolfgang Paul.  Partition function zeros and phase transitions for a square-well polymer chain.  Phys. Rev. E, 88:012604, Jul 2013.
  * [12] Kai Qi and Michael Bachmann.  Classification of phase transitions by microcanonical inflection-point analysis.  Phys. Rev. Lett., 120:180601, Apr 2018.
  * [13] S.G. Brush and N.S. Hall.  The Kinetic Theory of Gases: An Anthology of Classic Papers with Historical Commentary.  History of modern physical sciences. Imperial College Press, 2003.
  * [14] M. Planck.  Phase transitions and the distribution of temperature zeros of the partition function.  Annalen der Physik, 4:553–562, 1901.
  * [15] B. V. Costa, L. A. S. Mól, and J. C. S. Rocha.  Energy probability distribution zeros: A route to study phase transitions.  Computer Physics Communications, 216:77–83, 2017.
  * [16] J. J. Carvalho and A. L. Mota.  Finding the dominant zero of the energy probability distribution.  International Journal of Modern Physics C, 32(12):2150155, 2021.
  * [17] 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.
  * [18] Pavel A. Frantsuzov and Vladimir A. Mandelshtam.  Size-temperature phase diagram for small lennard-jones clusters.  Phys. Rev. E, 72:037102, Sep 2005.
  * [19] Julio C. S. Rocha, Stefan Schnabel, David P. Landau, and Michael Bachmann.  Leading fisher partition function zeros as indicators of structural transitions in macromolecules.  Physics Procedia, 57:94–98, 2014.  Proceedings of the 27th Workshop on Computer Simulation Studies in Condensed Matter Physics (CSP2014).
  * [20] Dario Andrea Bini and Giuseppe Fiorentino.  Design, analysis, and implementation of a multiprecision polynomial rootfinder.  Numerical Algorithms, 23(2):127–173, 2000.
  * [21] Dario A. Bini and Leonardo Robol.  Solving secular and polynomial equations: A multiprecision algorithm.  Journal of Computational and Applied Mathematics, 272:276–292, 2014.
  * [22] D. Landau and K. Binder.  A Guide to Monte Carlo Simulations in Statistical Physics.  Cambridge University Press, 2021.
  * [23] Fugao Wang and D. P. Landau.  Efficient, multiple-range random walk algorithm to calculate the density of states.  Phys. Rev. Lett., 86:2050–2053, Mar 2001.
  * [24] J. A. Northby.  Structure and binding of lennard‐jones clusters: $13\leq n\leq 147$.  The Journal of Chemical Physics, 87(10):6166–6177, 1987.
  * [25] Qiliang Yan and Juan J. de Pablo.  Fast calculation of the density of states of a fluid by monte carlo simulations.  Phys. Rev. Lett., 90:035701, Jan 2003.
  * [26] R. E. Belardinelli and V. D. Pereyra.  Fast algorithm to calculate density of states.  Phys. Rev. E, 75:046701, Apr 2007.
  * [27] A. A. Caparica and A. G. Cunha-Netto.  Wang-landau sampling: Improving accuracy.  Phys. Rev. E, 85:046702, Apr 2012.
  * [28] B.B. Rodrigues, J.C.S. Rocha, and B.V. Costa.  Phase diagram of flexible polymers with quenched disordered charged monomers.  Physica A: Statistical Mechanics and its Applications, 604:127787, 2022.
  * [29] Chenggang Zhou and Jia Su.  Optimal modification factor and convergence of the wang-landau algorithm.  Phys. Rev. E, 78:046705, Oct 2008.
  * [30] R. E. Belardinelli, S. Manzi, and V. D. Pereyra.  Analysis of the convergence of the 1/t and wang-landau algorithms in the calculation of multidimensional integrals.  Phys. Rev. E, 78:067701, Dec 2008.
  * [31] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics, 21(6):1087–1092, 1953.
  * [32] Thomas Vogel, Ying Wai Li, Thomas Wüst, and David P. Landau.  Generic, hierarchical framework for massively parallel wang-landau sampling.  Phys. Rev. Lett., 110:210603, May 2013.
  * [33] Bernd A. Berg and Thomas Neuhaus.  Multicanonical ensemble: A new approach to simulate first-order phase transitions.  Phys. Rev. Lett., 68:9–12, Jan 1992.