# A novel Bayesian approach to the computation
of the configurational density of states

Felipe Moreno<EMAIL_ADDRESS>Sergio Davis<EMAIL_ADDRESS>Joaquín Peralta<EMAIL_ADDRESS>Departamento de Física, Facultad de
Ciencias Exactas, Universidad Andrés Bello, Sazié 2212, Santiago, Chile.
Research Center in the Intersection of Plasma Physics, Matter and Complexity
(P2mc), Comisión Chilena de Energía Nuclear, Casilla 188-D, Santiago 7600713,
Chile.

###### Abstract

In this work we develop and implement a novel Bayesian method for computing
the DOS of a system. This method is based on the use of a test function with
adjustable parameters and we use Bayes theorem to find the best parameters
given a certain number of measurements done on the system. This measurements
can be done in any ensemble defined by a distribution function. We found that
the algorithm can find the DOS in a reasonable amount of time, and that if the
test function is suitable enough, the DOS found by the algorithm is very close
to the true DOS.

###### keywords:

Density of States, Bayes Theorem, Algorithm

##  1 Introduction

The density of states (DOS), on its different flavors, such like energy,
vibrational, configurational, among others [1], is essentially the number of
states per unit of energy. The configurational density of states, defined for
a system with degrees of freedom $\bm{x}$ and potential energy $\Phi(\bm{x})$
is given by

$D(\phi)\mathrel{\mathop{:}}=\int d\bm{x}\delta(\Phi(\bm{x})-\phi).$ (1)

In physics and chemistry this DOS have an outstanding importance, because a
large amount of information can be obtained from it. In spite of the enormous
success in characterizing physical and chemical quantities based on the DOS
[2, 3, 4, 5, 6, 7], it is not a easy to obtain function, which reduce its use,
in spite of the their characteristics, on multiples research areas [8, 9, 10].

Different approaches has been implemented to find the DOS for a large variety
of systems [11, 12, 13, 14]. Nowadays one of the most successful techniques
correspond to determine the DOS is the one based in Wang Landau algorithm
[14], this technique is based on flat-histogram procedure and have been showed
that plenty of programming techniques and improvements has been implemented
during last decades [15, 16, 17]. However it is not easy to calculate nor
programming, and is even more difficult to apply on specific scenarios, such
like ab-initio.

While Wang-Landau has been very popular this time, there are different
techniques for calculating DOS that are still in use, such as parallel
tempering [18], metadynamics [19, 20, 21], and other techniques based on flat
histograms [9, 22]. In most of these cases, solution are focused on specific
problems, and not always cover the variety of cases where DOS could be used.
Motivated on this, and also in the incorporation of most general statistical
techniques, we worked in a new methodology to determine the DOS.

In this work we present a novel procedure to obtain the energetic density of
states (DOS) for solid materials, based on Bayesian method. Here we develop
this novel technique by: i) incorporate the Bayesian Theorem on the procedure,
ii) develop a modular implementation on a wide and popular programming
language as Python3, and iii) the speed-up of certain function using C
language.

In what follows we present the Algorithm description which connect the
Bayesian theorem and the DOS (sec. 2). On section 3 a brief description of the
implementation of the algorithm is presented using Python programming
language. Section 4 present cases of study for different scenarios, such like
classical solid structure using a Lennard-Jones interatomic potential, and the
Ising model. Finally on section 5 conclusions and scope of our work is
discussed.

##  2 Algorithm

Suppose the DOS of a system is completely determined by $m$ parameters
$\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{m})$

$D(\phi)=D(\phi;\boldsymbol{\theta}).$ (2)

The probability that an ensemble $S$ has a certain energy $\phi$ is

$P(\phi;S|\boldsymbol{\theta})=\frac{\rho(\phi;S)D(\phi;\boldsymbol{\theta})}{Z(\boldsymbol{\theta};S)}$
(3)

where $\rho(\phi;S)$ is the distribution that defines the ensemble $S$ and
$Z(\boldsymbol{\theta};S)$ is the partition function defined by

$Z(\boldsymbol{\theta};S)=\int
d\phi\,\rho(\phi;S)D(\phi;\boldsymbol{\theta}).$ (4)

Suppose we make $n$ measurements of the ensemble energy, obtaining
$\phi_{1},\ldots,\phi_{n}$. We want to know what is the most probable values
of the parameters $\boldsymbol{\theta}$ compatible with those measurements. We
accomplish this by computing the maximum of
$P(\boldsymbol{\theta}|\phi_{1},\ldots,\phi_{n};S)$. In practice, it is better
to compute the maximum of the logarithm of
$P(\boldsymbol{\theta}|\phi_{1},\ldots,\phi_{n};S)$, that is,

$\partialderivative{\theta_{k}}\ln
P(\boldsymbol{\theta}|\phi_{1},\ldots\phi_{n};S)=0\qquad\qquad k=1,\ldots,m.$
(5)

Because the measurements are considered to be statistically independent, we
have

$P(\phi_{1},\ldots,\phi_{n}|\boldsymbol{\theta},S)=\prod_{i=1}^{n}\frac{\rho(\phi_{i};S)D(\phi_{i};\boldsymbol{\theta})}{Z(\boldsymbol{\theta};S)}$
(6)

and, by using Bayes theorem [23] we have

$\begin{split}P(\boldsymbol{\theta}|\phi_{1},\ldots\phi_{n};S)&=\frac{P(\boldsymbol{\theta}|I_{0})P(\phi_{1},\ldots,\phi_{n}|\boldsymbol{\theta},S)}{P(\phi_{1},\ldots,\phi_{n};S)}\\\
&=\frac{P(\boldsymbol{\theta}|I_{0})}{P(\phi_{1},\ldots,\phi_{n};S)}\prod_{i=1}^{n}\frac{\rho(\phi_{i};S)D(\phi_{i};\boldsymbol{\theta})}{Z(\boldsymbol{\theta};S)}\\\
&=\frac{1}{\eta}P(\boldsymbol{\theta}|I_{0})\frac{\prod_{i=1}^{n}D(\phi_{i};\boldsymbol{\theta})}{Z(\boldsymbol{\theta};S)^{n}}\end{split}$
(7)

where

$\eta\mathrel{\mathop{:}}=\frac{\prod_{i=1}^{n}\rho(\phi_{i};S)}{P(\phi_{1},\ldots,\phi_{n}|I_{0})}$
(8)

is an irrelevant normalization constant that does not depend on
$\boldsymbol{\theta}$.

In this way, Eq. 5 can be written as

$\partialderivative{\theta_{k}}\ln
P(\boldsymbol{\theta}|\phi_{1},\ldots\phi_{n};S)=\partialderivative{\theta_{k}}\ln
P(\boldsymbol{\theta}|I_{0})+\sum_{i=1}^{n}\partialderivative{\theta_{k}}\ln
D(\phi_{i};\boldsymbol{\theta})-n\partialderivative{\theta_{k}}\ln
Z(\boldsymbol{\theta};S)=0.$ (9)

In the limit $n\rightarrow\infty$ the first term on the right-hand side of Eq.
9 is negligible, and Eq. 5 becomes

$\partialderivative{\theta_{k}}\ln
Z(\boldsymbol{\theta};S)=\frac{1}{n}\sum_{i=1}^{n}\partialderivative{\theta_{k}}\ln
D(\phi_{i};\boldsymbol{\theta})$ (10)

for $k=1,\ldots,m$. Notice that the right-hand side is the average of a known
function of the measurements, while the left-hand side is given by

$\begin{split}\partialderivative{\theta_{k}}\ln
Z(\boldsymbol{\theta};S)&=\frac{1}{Z(\boldsymbol{\theta};S)}\int
d\phi\,\rho(\phi;S)\partialderivative{\theta_{k}}D(\phi;\boldsymbol{\theta})\\\
&=\frac{1}{Z(\boldsymbol{\theta})}\int
d\phi\,\rho(\phi;S)D(\phi;\boldsymbol{\theta})\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})\\\
&=\expectationvalue{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}_{S},\end{split}$ (11)

that is, the expectation value over the ensemble of the same function that is
averaged on the right-hand side. In other words, we have that the condition
that fixes $\boldsymbol{\theta}$ is simply

$\expectationvalue{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}_{S}=\overline{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}$ (12)

for $k=1,\ldots,m$. This condition can be understood as the validity of the
law of large numbers, namely the identification of the expected value with the
data average, for the logarithmic derivatives of the density of states.

A straightforward generalization is to consider $r$ groups of
$n_{1},\ldots,n_{r}$ measurements obtained from different ensembles
$S_{1},\ldots,S_{r}$ respectively. Then Eq. 7 generalizes to

$P(\boldsymbol{\theta}|\phi_{1},\ldots\phi_{n};S_{1},\ldots,S_{r})=\frac{1}{\eta}P(\boldsymbol{\theta})\frac{\prod_{i=1}^{N}D(\phi_{i};\boldsymbol{\theta})}{\prod_{j=1}^{r}[Z(\boldsymbol{\theta};S_{j})]^{n_{j}}}$
(13)

where $N=\sum_{j=1}^{r}n_{j}$ is the total number of measurements, and Eq. 12
generalizes to

$\sum_{j=1}^{r}\frac{n_{j}}{N}\expectationvalue{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}_{S_{j}}=\overline{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}$ (14)

for $k=1,\ldots,m$ where the average in the right hand side is performed over
all $N$ measurements. This algorithm is based on the solution of Eq. 14 by
numerical methods. An equivalent problem is to search for the global minimum
of

$\epsilon^{2}(\boldsymbol{\theta})=\sum_{k=1}^{m}\left[\sum_{j=1}^{r}\frac{n_{j}}{N}\expectationvalue{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}_{S_{j}}-\overline{\partialderivative{\theta_{k}}\ln
D(\phi;\boldsymbol{\theta})}\right]^{2}$ (15)

by using numerical minimization methods such as particle swarm optimization
(PSO) [24].

Now we can formulate a practical method of computing the DOS, using Eq. 15.
Firstly, we select a set of canonical ensembles $S_{1},\ldots,S_{r}$ such that

$\rho(\phi;S_{j})\propto\exp(-\beta_{j}\phi)$

and we take $n_{j}$ samples of energy respectively. The ensemble information
($\beta_{j}$) and the measurements ($\phi_{i}$) are stored. Secondly, we
choose a suitable function $\ln D(\phi;\boldsymbol{\theta})$ and then we run a
minimization algorithm on the sampled data in order to find the global minimum
of $\epsilon^{2}(\boldsymbol{\theta})$. Once the global minimum has been
found, the parameters $\boldsymbol{\theta}$ that minimize
$\epsilon^{2}(\boldsymbol{\theta})$ are the ones that maximize the probability
that $D(\phi;\boldsymbol{\theta})$ is the true DOS, within the chosen
parametrization. If the function chosen is suitable enough, we expect this
function to be very close to the true DOS.

Alternatively, we could have chosen a different set of ensembles, for
instance, the microcanonical configurational ensemble

$\rho(\phi;S_{j})\propto(E_{j}-\phi)^{\frac{d\mathcal{N}}{2}-1}$

for different total energies $E_{1},\ldots,E_{r}$ where $d$ is the number of
spatial dimensions.

##  3 Implementation

The algorithm consists of two parts. The first one is to generate measurements
of the ensembles energy. We do this by using a canonical Monte Carlo
Metropolis algorithm. We set an initial state
$\mathbf{x}=(x_{1},\ldots,x_{\mathcal{N}})$ and we propose a transition to a
new state $\mathbf{x}^{\prime}$ which is accepted with probability

$P(\mathbf{x}\rightarrow\mathbf{x}^{\prime}|\beta)=\min\left(1,e^{-\beta(\Phi(\mathbf{x}^{\prime})-\Phi(\mathbf{x}))}\right)$
(16)

where $\Phi(\mathbf{x})$ is the configurational energy of state $\mathbf{x}$.
We vary $\beta$ in order to produce data for different ensembles. We take
samples of the ensemble energy after a certain number of calibration steps
have been performed following each variation of $\beta$. This part of the
algorithm will generate two files, one with the set of temperatures
$\beta_{j}$ that define each ensemble, and one with the sampled energies.

The second part of the algorithm is to run a PSO algorithm on the sampled data
in order to find the global minimum of $\epsilon^{2}(\boldsymbol{\theta})$.
Since the computation of $\epsilon^{2}$ is the computationally expensive part
of the algorithm, it has been coded in C++ and loaded as a library into a
easier to handle Python code.

Furthermore, we can speed up the minimization by considering a very common
special case. Suppose that $\ln D(\phi,\boldsymbol{\theta})$ depends linearly
on the parameters $\boldsymbol{\theta}$, that is

$\ln D(\phi;\boldsymbol{\theta})=\sum_{k=1}^{m}\theta_{k}f_{k}(\phi)$ (17)

where we can interpret the set of functions
$\\{f_{1}(\phi),\ldots,f_{m}(\phi)\\}$ as a basis for the representation of
the density of states. Inserting this into Eq. 15 we obtain

$\epsilon^{2}(\boldsymbol{\theta})=\sum_{k=1}^{m}\left\\{\sum_{j=1}^{r}\frac{n_{j}}{N}\expectationvalue{f_{k}(\phi)}_{S_{j}}-\overline{f_{k}(\phi)}\right\\}^{2},$
(18)

where the terms $\overline{f_{k}(\phi)}$ do not depend on
$\boldsymbol{\theta}$ thus they will only need to be computed once. Since
these are the only terms that depend on the sampled data, we can make the
minimization time independent of the size of this data.

To estimate the statistical error of the DOS reconstruction, we peopose the
use of the Bayesian information criteria (BIC) [25], which is defined by:

$\text{BIC}\mathrel{\mathop{:}}=-2\ln
P(\phi_{1},\ldots,\phi_{N}|\boldsymbol{\theta})+m\ln N$ (19)

and such that we should prefer the model with the lowest BIC to reduce the
error. By using Eq. 6 we have

$\text{BIC}=-2N\Big{(}\overline{\ln\rho(\phi;S)}+\overline{\ln
D(\phi;\boldsymbol{\theta})}-\ln Z(\boldsymbol{\theta};S)\Big{)}+m\ln N.$ (20)

Because the first term inside the parenthesis does not depend on the
parameterization $\boldsymbol{\theta}$ we can instead compare two models
$D_{1}(\phi;\bm{\theta}_{1})$ and $D_{2}(\phi;\boldsymbol{\theta}_{2})$ for
the same set of ensembles and same data through the quantity

$\text{BIC}^{*}=-2\Big{(}N\overline{\ln
D(\phi;\boldsymbol{\theta})}-\sum_{j=1}^{r}n_{j}\ln
Z(\boldsymbol{\theta};S_{j})\Big{)}+m\ln N.$ (21)

##  4 Results

The algorithm and implementation have been tested using a variety of systems
in order to evaluate the accuracy of the suggested methodology.

###  4.1 Harmonic Oscillator

We consider a system of $\mathcal{N}=100$ identical harmonic oscillators,

$\Phi(\mathbf{x})=\frac{1}{2}m\omega^{2}\sum_{i=1}^{\mathcal{N}}x_{i}^{2}.$
(22)

Instead of $\beta$ we have used $k_{B}T=\beta^{-1}$ in the range
$k_{B}T\in[0.05,5]$ at steps $\Delta(k_{B}T)=0.05$. We consider the energy of
the last $1\,000$ steps of simulation—in a total of 10 000—for each $k_{B}T$.
That is, we produce $r=100$ values of $k_{B}T$ and $N=100\,000$ values of
energy $\phi$.

Next, we chose the following test function with one parameter ($m=1$),

$\ln D(\phi;\boldsymbol{\theta})=\theta_{1}\ln\phi.$ (23)

Figure 1: Configurational DOS for a system of $\mathcal{N}$=100 classical
harmonic oscillators.

We run the PSO algorithm on this data using 1760 particles and after 37
iterations we found that $\theta_{1}=49.001$, which closely matches the
analytical result

$D(\phi)\propto\phi^{\mathcal{N}/2-1}.$

The results for this case are presented on figure 1, we observe an excellent
agreement with the WL algorithm result.

###  4.2 Ising model

The Ising model is a well-known model of spins. The model consists in a
lattice filled with spins that can take values $\sigma=\pm 1$. The model
energy is computed according to

$\phi=-J\sum_{\langle i,j\rangle}\delta(\sigma_{i},\sigma_{j}).$ (24)

where the sum is over adjacent pairs and $\delta(\sigma_{i},\sigma_{j})$ is 1
if $\sigma_{i}=\sigma_{j}$ and 0 otherwise. We use an Ising model on a
$10\times 10$ lattice and we approximate the DOS using two 1-parameter test
functions. Since the DOS is defined up to a constant, we require that the test
functions vanish at $\phi=\pm 200J$ which are the energies less populated. The
test functions will be

$\displaystyle\log D_{1}(\phi;\boldsymbol{\theta})$
$\displaystyle=\theta_{1}\left[\left(\frac{\phi}{200}\right)^{2}-1\right]$
(25) $\displaystyle\log D_{2}(\phi;\boldsymbol{\theta})$
$\displaystyle=\theta_{1}\cos\left(\frac{\phi}{2\cdot 200}\right)$ (26)

We found the parameters $\theta_{1}=-69.2799$ for the first function and
$\theta_{1}=67.1554$ for the second one. The BIC is smaller for the cosine
function $3.40277\times 10^{7}<3.41584\times 10^{7}$. The results for this
case are presented on figure 2, we observe a better agreement with WL
algorithm on the function with the lower BIC value.

Figure 2: Configurational DOS for an Ising system of 100 spins on a $10\times
10$ lattice using a quadratic and a cosine approximation for $\ln D(\phi)$.

###  4.3 Lennard-Jones

Finally, we study a system that needs 2 parameters to get a good
approximation. A Lennard–Jones gas is a gas of particles whose interaction
energy is

$\phi=\sum_{i,j}4\varepsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{6}-\left(\frac{\sigma}{r_{ij}}\right)\right],$
(27)

where the sum is over every pair. For this system we consider the addition of
a logarithm and a linear function.

$\log
D(\phi;\boldsymbol{\theta})=\theta_{1}\log(\phi-\phi_{0})+\theta_{2}\phi$ (28)

where $\phi_{0}$ is the minimum energy.

We found $\theta_{1}=1.52994\times 10^{2}$ and $\theta_{2}=1.09103$. The
drawback here is that this model is unable to produce the phase transition. In
spite of this, as we can observe on figure 3, there is a clear insight in the
figure whit a results that behave different on the two main branches of the WL
algorithm. For energies lower than -530 , the Bayesian estimation is above the
WL approximation. On the other side, with energies above the -530 the Bayesian
estimation is below the WL results. This suggest that our model it could be
splitted on different branches when we know beforehand about the presence of a
phase transition.

Figure 3: Configurational DOS for a Lennard-Jones system of 120 particles
assuming $\log D(\phi)$ is of the form shown in Eq. 28.

##  5 Conclusions

We have presented a new technique to determine the Density of States.
Different simulation examples have been used: Ising model, Lenard-Jones, and
classical harmonic oscillator. Our results shows a good agreement with the
known theoretical results. Moreover this methodology shows an improvement on
computational efforts to be required to determine the DOS. The requirement of
an input function, its an strength when some insights are known beforehand,
otherwise a large functional base can always be used. If the functional form
is improved the BIC error procedure will return the most effective evaluation
of the DOS.

##  6 Acknowledgments

FM thanks to Universidad Andrés Bello for providing the doctorate scholarship
and acknowledges funding from ANID BECAS/DOCTORADO NACIONAL 21212267. SD
acknowledges partial funding from ANID PIA ACT172101 grant. Computational work
was supported by the supercomputing infrastructures of the NLHPC (ECM-02), and
FENIX (UNAB).

## References

  * [1] W.A. Harrison.  Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond.  Dover Books on Physics. Dover Publications, 1989.
  * [2] C. J. Gommes, Y. Jiao, and S. Torquato.  Density of states for a specified correlation function and the energy landscape.  Phys. Rev. Lett., 108:080601, Feb 2012.
  * [3] P. Sarker, T. Harrington, C. Toher, C. Oses, M. Samiee, J.-P. Maria, D. W. Brenner, K. S. Vecchio, and S. Curtarolo.  High-entropy high-hardness metal carbides discovered by entropy descriptors.  Nature Communications, 9:1–10, 2018.
  * [4] E. J. Garboczi and M. F. Thorpe.  Density of states for random-central-force elastic networks.  Phys. Rev. B, 32:4513–4518, Oct 1985.
  * [5] John C. Mauro, Roger J. Loucks, Jitendra Balakrishnan, and Srikanth Raghavan.  Monte carlo method for computing density of states and quench probability of potential energy and enthalpy landscapes.  The Journal of Chemical Physics, 126(19):194103, 2007.
  * [6] S.-Q. Yu, J.-B. Wang, D. Ding, S.R. Johnson, D. Vasileska, and Y.-H. Zhang.  Impact of electronic density of states on electroluminescence refrigeration.  Solid-State Electronics, 51(10):1387–1390, 2007.  Special Issue: Papers Selected from the NGC2007 Conference.
  * [7] N. Nakai, P. Miranović, M. Ichioka, and K. Machida.  Theoretical study on the field dependence of the zero energy density of states in an anisotropic gap superconductors.  Journal of Physics and Chemistry of Solids, 66(8):1362–1364, 2005\.  Proceedings of the ISSP International Symposium (ISSP-9)on Quantum Condensed System.
  * [8] Helge Rosé, Werner Ebeling, and Torsten Asselmeyer.  The density of states — a measure of the difficulty of optimisation problems.  In Hans-Michael Voigt, Werner Ebeling, Ingo Rechenberg, and Hans-Paul Schwefel, editors, Parallel Problem Solving from Nature — PPSN IV, pages 208–217, Berlin, Heidelberg, 1996. Springer Berlin Heidelberg.
  * [9] N. Rathore, T. A. Knotts IV, and J. J. de Pablo.  Configurational temperature density of states simulations of proteins.  Biophys. J., 85:3963–3968, 2003.
  * [10] G. J. Rodgers and A. J. Bray.  Density of states of a sparse random matrix.  Phys. Rev. B, 37:3557–3562, Mar 1988.
  * [11] M S Methfessel, M H Boon, and F M Mueller.  Analytic-quadratic method of calculating the density of states.  Journal of Physics C: Solid State Physics, 16:L949–L954, 09 1983\.
  * [12] Hainam Do and Richard J. Wheatley.  Density of states partitioning method for calculating the free energy of solids.  Journal of Chemical Theory and Computation, 9(1):165–171, 2013\.  PMID: 26589019.
  * [13] Thomas Aquino Ariasoca, Sholihun, and Iman Santoso.  Trotter-suzuki-time propagation method for calculating the density of states of disordered graphene.  Computational Materials Science, 156:434–440, 2019.
  * [14] F. Wang and D. P. Landau.  Efficient, multiple-range random walk algorithm to calculate the density of states.  Phys. Rev. Lett., 86:2050–2053, 2001.
  * [15] P. Dayal, S. Trebst, S. Wessel, D. Würtz, M. Troyer, S. Sabhapandit, and S. N. Coppersmith.  Performance limitations of flat-histogram methods.  Phys. Rev. Lett., 92:097201, Mar 2004.
  * [16] Thomas Prellberg and Jarosław Krawczyk.  Flat histogram version of the pruned and enriched rosenbluth method.  Phys. Rev. Lett., 92:120602, Mar 2004.
  * [17] Faming Liang.  A theory on flat histogram monte carlo algorithms.  Journal of Statistical Physics, 122:511–529, 01 2006.
  * [18] D. J. Earl and M. W. Deem.  Parallel tempering: Theory, applications, and new perspectives.  Phys. Chem. Chem. Phys., 7:3910–3916, 2005.
  * [19] C. Micheletti, A. Laio, and M. Parrinello.  Reconstructing the density of states by history-dependent metadynamics.  Phys. Rev. Lett., 92:170601, 2004.
  * [20] A. Laio and F. L. Gervasio.  Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and materials science.  Rep. Prog. Phys., 71:126601, 2008.
  * [21] A. Barducci, G. Bussi, and M. Parrinello.  Well-tempered metadynamics: A smoothly converging and tunable free-energy method.  Phys. Rev. Lett., 100:020603, 2008.
  * [22] N. Rathore, T. A. Knotts, and J. J. de Pablo.  Density of states simulations of proteins.  J. Chem. Phys., 118(9):4285–4290, 2003.
  * [23] D. S. Sivia and J. Skilling.  Data Analysis: A Bayesian Tutorial.  Oxford University Press, 2006.
  * [24] A. E. Olsson.  Particle Swarm Optimization: Theory, techniques and applications.  Nova Science Publishers, 2011.
  * [25] G. E. Schwarz.  Estimating the dimension of a model.  Ann. Statist., 6:461–464, 1978.
  * [26] M. Habeck.  Bayesian reconstruction of the density of states.  Phys. Rev. Lett., 98:200601, 2007.
  * [27] F. Moreno, S. Davis, C. Loyola, and J. Peralta.  Ordered metastable states in the Potts model and their connection with the superheated solid state.  Phys. A, 509:361–368, 2018.
  * [28] C. Farías and S. Davis.  Multiple metastable states in an off-lattice Potts model.  Phys. A, 581:126215, 2021.
  * [29] S. Davis, J. Jain, and B. Bora.  Computational statistical mechanics of a confined, three-dimensional Coulomb gas.  Phys. Rev. E, 102:042137, 2020.


FIRST-ORDER PHASE TRANSITION IN THE FRAMEWORK OF THE CELL FLUID MODEL: REGIONS
OF CHEMICAL POTENTIAL VARIATION AND THE CORRESPONDING DENSITIES

I.V. Pylyuk111e-mail<EMAIL_ADDRESS>M.P. Kozlovskii

Institute for Condensed Matter Physics
of the National Academy of Sciences of Ukraine,
1 Svientsitskii Str., 79011 Lviv, Ukraine

A microscopic description is given for the behavior of the fluid system in an
immediate vicinity of its critical point, where theoretical and experimental
researches are difficult to carry out. For the temperatures $T<T_{c}$, the
regions of chemical potential and density variations are singled out and
analyzed. The equation of state of the cell fluid model in terms of
temperature-chemical potential is written using the Heaviside functions. This
equation is also given in terms of the temperature and density variables. As a
result of the study of the relationship between the density and the chemical
potential, an equation for the binodal curve is obtained in a narrow
neighborhood of the critical point.

PACS numbers: 05.70.Ce, 64.60.F-, 64.70.F-

Keywords: cell fluid model, chemical potential, density, equation of state,
binodal

##  1 Introduction

The theoretical and experimental study of the behavior of liquids and their
mixtures in a vicinity of their critical point (see, e.g., works [1, 2, 3, 4,
5, 6, 7, 8]) is an important and challenging task. In our previous works [9,
10], the behavior of fluid was studied in an immediate vicinity of the
critical point, and in works [11, 12, 13] beyond this vicinity. As a result, a
wide region near the critical point has been covered. The mathematical
description was carried out in the framework of the cell fluid model using the
grand canonical ensemble. The whole volume $V$ of a system consisting of $N$
interacting particles was conventionally divided into $N_{v}$ cells, each of
the volume $v=V/N_{v}=c^{3}$, where $c$ is the linear cell size. Note that
unlike the lattice gas model where the cell is assumed to contain no more than
one particle, in this approach the cell can include more than one particle.

In works [9, 10], the analysis was performed in the framework of the
collective-variables approach using the renormalization group transformation
[14]. An analytical procedure for calculating the grand partition function and
the thermodynamic potential of the cell fluid model was developed in works [9,
10] in the approximation of a non-Gaussian (quartic) distribution for order
parameter fluctuations without involving the hard-spheres reference system.
The formation of the reference system as a part of the repulsive component of
the interaction potential made it possible to take into account all kinds of
interaction (both short- and long-range) from the same position of the
collective-variables approach.

The role of interaction potential in this work is played by the Morse
potential. The interaction potential parameters, which are given in works [9,
10] and are necessary for quantitative estimates, correspond to the data for
sodium. They were taken from work [15], which was devoted to the study of
vapor-liquid equilibrium curves for metals using Monte Carlo simulation and
Morse potential. Taking this potential into account, the vapor-liquid
coexistence curves for metals were also analyzed in the framework of the
integral-equations approach [16].

The Morse potential is widely used when studying the melting and laser
ablation processes using computer simulation [17, 18, 19]. There are works in
the literature devoted to the study of the structural properties of the Morse
and Lennard-Jones clusters, as well as to a comparison between them [19, 20,
21, 22]. In some cases, modifications of the Morse function were made [23, 24,
25, 26] in order to improve the numerical results. Although the Morse
potential was traditionally used to model covalently bound diatomic molecules
[27, 28, 29], it is also applied to estimate non-bounding interactions [30,
31]. This potential is qualitatively similar to the Lennard-Jones potential,
but they are quite different from the quantitative point of view. The Lennard-
Jones and Morse potentials can be compared directly using a mathematical
relationship enabling the point of energy minimum to be located at the same
position [32, 33]. In addition, it was shown that either of the potentials can
be derived from the other [34].

This work complements the study of the critical behavior of the Morse fluid
that was performed in works [9, 10]. In particular, for temperatures lower
than the critical one, solutions of a certain cubic equation are obtained.
They govern the quantities entering the equation of state of the cell fluid
model. The behavior of the equation solutions depending on the chemical
potential value is analyzed in the immediate vicinity of the critical point,
and the regions where the chemical potential changes are singled out. Each of
those regions is considered separately. Expressions for the boundary densities
(the densities at the regions’ boundaries) are obtained. For this purpose, a
nonlinear equation, which relates the density to the chemical potential, is
used. The equation of state of the cell fluid model and the binodal equation
are presented.

##  2 Chemical potential and density changes at temperatures below the
critical one

The equation of state of the cell fluid model for temperatures $T<T_{c}$ [10]
contains quantities dependent on the solution $\sigma^{\prime}_{0}$ of the
equation

$(\sigma^{\prime}_{0})^{3}+p^{\prime}\sigma^{\prime}_{0}+q^{\prime}=0.$ (1)

(the solution is also given in work [10]). The coefficients

$p^{\prime}=6\frac{r_{n^{\prime}_{p}+2}}{u_{n^{\prime}_{p}+2}},\quad
q^{\prime}=-6\frac{s^{5/2}}{u_{n^{\prime}_{p}+2}}\frac{\tilde{h}}{(\tilde{h}^{2}+h^{2}_{cm})^{1/2}}$

include the quantities $r_{n^{\prime}_{p}+2}$ and $u_{n^{\prime}_{p}+2}$
determining the long-wavelength part of the grand partition function of the
model. The quantity $\tilde{h}$ is proportional to the chemical potential $M$,
and the quantity $h_{cm}$ is characterized by the renormalized relative
temperature $\tau=(T-T_{c})/T_{c}$ (see work [10]). The renormalization group
parameter $s$ determines the separation of the phase space of collective
variables into layers.

The form of the solutions of Eq. (1) depends on the sign of the discriminant

$Q=(p^{\prime}/3)^{3}+(q^{\prime}/2)^{2}.$ (2)

If $Q>0$, the single real solution $\sigma^{\prime}_{0}$ of Eq. (1), according
to Cardano’s formula, looks like

$\displaystyle\sigma^{\prime}_{0b}=A+B,$ $\displaystyle
A=(-q^{\prime}/2+Q^{1/2})^{1/3},$ $\displaystyle
B=(-q^{\prime}/2-Q^{1/2})^{1/3}.$ (3)

If $Q<0$, there are three real solutions (the quantity $\sigma^{\prime}_{0}$
acquires three possible real values)

$\displaystyle\sigma^{\prime}_{01}=2(-p^{\prime}/3)^{1/2}\cos(\alpha_{r}/3),$
$\displaystyle\sigma^{\prime}_{02,03}=-2(-p^{\prime}/3)^{1/2}\cos(\alpha_{r}/3\pm\pi/3),$
(4)

where $\alpha_{r}$ is determined from the equation

$\cos\alpha_{r}=-\frac{q^{\prime}}{2(-p^{\prime}/3)^{3/2}}.$ (5)

If the discriminant is negative, solutions (4) can be rewritten as follows:

$\displaystyle\sigma^{\prime}_{01}=2\sigma_{0r}\cos\frac{\alpha_{r}}{3},\quad\sigma^{\prime}_{02}=-2\sigma_{0r}\cos\left(\frac{\alpha_{r}}{3}+\frac{\pi}{3}\right),$
$\displaystyle\sigma^{\prime}_{03}=-2\sigma_{0r}\cos\left(\frac{\alpha_{r}}{3}-\frac{\pi}{3}\right).$
(6)

Here,

$\displaystyle\sigma_{0r}=\left(-\frac{2r_{n^{\prime}_{p}+2}}{u_{n^{\prime}_{p}+2}}\right)^{1/2},$
$\displaystyle\alpha_{r}=\arccos\left(\frac{M\left(\tilde{h}^{2}_{q}+h^{2}_{cm}\right)^{1/2}}{M_{q}\left(\tilde{h}^{2}+h^{2}_{cm}\right)^{1/2}}\right).$
(7)

The chemical potential $M_{q}$ is determined from the condition $Q=0$ and
satisfies the equality

$M_{q}=\Biggl{[}-\frac{8r^{3}_{n^{\prime}_{p}+2}(1+\alpha^{2}_{mq})}{9u_{n^{\prime}_{p}+2}s^{5}\beta
W(0)}\Biggr{]}^{1/2}h_{cm},$ (8)

where

$\alpha_{mq}=\tilde{h}_{q}/h_{cm},\quad\tilde{h}_{q}=M_{q}\left(\beta
W(0)\right)^{1/2}.$ (9)

Here, $\beta=1/(kT)$ is the inverse temperature, and $W(0)$ is the Fourier
transform of the effective interaction potential [9] at the zero wave vector.
For all $|M|<M_{q}$, the discriminant $Q<0$ and, therefore, there are three
real roots of Eq. (1) in this interval of $M$-values.

The dependence of the solutions of the cubic equation (1) on the chemical
potential $M$ at $T<T_{c}$ is shown in Figs. 1 and 2.

Figure 1: Solutions of cubic equation (1) as functions of the chemical
potential $M$ at $\tau=-0.005$. Curves 1, 2, 3, and 4 correspond to
$\sigma^{\prime}_{0}=\sigma^{\prime}_{01}$,
$\sigma^{\prime}_{0}=\sigma^{\prime}_{02}$,
$\sigma^{\prime}_{0}=\sigma^{\prime}_{03}$, and
$\sigma^{\prime}_{0}=\sigma^{\prime}_{0b}$, respectively. Figure 2: Dependence
of the solutions of cubic equation (1) on the chemical potential $M$ at
various values of the relative temperature ($\tau=-0.005$, $\tau=-0.007$,
$\tau=-0.009$).

Figure 2 makes it possible to trace the shift of $M_{q}$ (the joint points of
the thin and thick solid curves) as $\tau$ changes. As one can see, the
absolute value of $M_{q}$ decreases with the reduction of $|\tau|$.

In the case $|M|>M_{q}$, as was at $T>T_{c}$, Eq. (1) has a single real
solution (3). For the latter, at $M=-M_{q}$ and taking into account the
equality $Q=0$, we obtain

$\sigma_{0bq}^{{}^{\prime}(-)}\\!=\\!2\left[-3s^{5/2}\frac{M_{q}(\beta
W(0))^{1/2}}{u_{n^{\prime}_{p}+2}h_{cm}(1+\alpha^{2}_{mq})^{1/2}}\right]^{1/3}\\!\\!\\!=\\!\\!-2\sigma_{0r},$
(10)

whereas at $M=M_{q}$ we have

$\sigma_{0bq}^{{}^{\prime}(+)}=2\left[3s^{5/2}\frac{M_{q}(\beta
W(0))^{1/2}}{u_{n^{\prime}_{p}+2}h_{cm}(1+\alpha^{2}_{mq})^{1/2}}\right]^{1/3}=2\sigma_{0r}.$
(11)

Let us analyze the asymptotics of solutions (6) at $|M|=M_{q}$. If $M=-M_{q}$,
we obtain $\cos\alpha_{rq}^{(-)}=-1$, $\alpha_{rq}^{(-)}=\pi$, and

$\displaystyle\sigma_{01}^{{}^{\prime}(-)}=2\sigma_{0r}\cos\frac{\pi}{3}=\sigma_{0r},$
$\displaystyle\sigma_{02}^{{}^{\prime}(-)}=-2\sigma_{0r}\cos\left(\frac{2\pi}{3}\right)=\sigma_{0r},$
$\displaystyle\sigma_{03}^{{}^{\prime}(-)}=-2\sigma_{0r}\cos 0=-2\sigma_{0r}.$
(12)

The case $M=M_{q}$ brings us to the formulas $\cos\alpha_{rq}^{(+)}=1$,
$\alpha_{rq}^{(+)}=0$, and

$\displaystyle\sigma_{01}^{{}^{\prime}(+)}=2\sigma_{0r}\cos 0=2\sigma_{0r},$
$\displaystyle\sigma_{02}^{{}^{\prime}(+)}=-2\sigma_{0r}\cos\frac{\pi}{3}=-\sigma_{0r},$
$\displaystyle\sigma_{03}^{{}^{\prime}(+)}=-2\sigma_{0r}\cos\left(-\frac{\pi}{3}\right)=-\sigma_{0r}.$
(13)

Thus, if $M=-M_{q}$, the solution is
$\sigma_{0bq}^{{}^{\prime}(-)}=-2\sigma_{0r}$ (see Eq. (10)), which coincides
with $\sigma_{03}^{{}^{\prime}(-)}$ from Eq. (12). If $M=M_{q}$, the solution
is $\sigma_{0bq}^{{}^{\prime}(+)}=2\sigma_{0r}$ (see Eq. (11)) and it
coincides with $\sigma_{01}^{{}^{\prime}(+)}$ from Eq. (13).

The conclusion drawn from the above calculations is as follows. As the
chemical potential increases to $-M_{q}$ from the side of negative values, Eq.
(1) has a single solution given by Eq. (3) (region I, gas phase, in Fig. 3).
At $-M_{q}<M<0$, it transforms into the solution

Figure 3: Regions of chemical potential variation and the corresponding
densities for temperatures below the critical one.

$\sigma^{\prime}_{03}$ from Eq. (6), which is valid up to $M=-0$ (region II,
transient gas phase). For $M=-M_{q}$, we obtain
$\sigma^{\prime}_{03}=\sigma_{03}^{{}^{\prime}(-)}=-2\sigma_{0r}$ (see Eq.
(12)), and for $M=-0$, we arrive at the expressions $\cos\alpha_{r0}^{(-)}=0$,
$\alpha_{r0}^{(-)}=\pi/2$, and

$\lim_{M\rightarrow-0}\sigma^{\prime}_{03}=\sigma_{030}^{{}^{\prime}(-)}=-2\sigma_{0r}\cos\left(-\frac{\pi}{6}\right)=-\sqrt{3}\sigma_{0r}.$
(14)

On the other hand, as $M$ decreases to $M_{q}$ from the side of positive
values, there is the single solution (3) (region IV, fluid phase). At
$0<M<M_{q}$, this solution transforms into the solution $\sigma^{\prime}_{01}$
from Eq. (6), which is valid up to $M=+0$ (region III, transient fluid phase).
For $M=M_{q}$, we have
$\sigma^{\prime}_{01}=\sigma_{01}^{{}^{\prime}(+)}=2\sigma_{0r}$ (see Eq.
(13)), and for $M=+0$ we obtain

$\lim_{M\rightarrow+0}\sigma^{\prime}_{01}=\sigma_{010}^{{}^{\prime}(+)}=2\sigma_{0r}\cos\left(\frac{\pi}{6}\right)=\sqrt{3}\sigma_{0r}.$
(15)

According to the chemical potential value $M$, the equation of state of the
cell fluid model at $T<T_{c}$ (see work [10]) can be written in the form

$\displaystyle\frac{Pv}{kT}=P_{a}^{(-)}(T)+E_{\mu}+D_{13}(\sigma^{\prime}_{0b})\left[\Theta(-M-M_{q})+\right.$
$\displaystyle\left.+\Theta(M-M_{q})\right]+D_{13}(\sigma^{\prime}_{03})\Theta(-M)\Theta(M+M_{q})+$
$\displaystyle+D_{13}(\sigma^{\prime}_{01})\Theta(M)\Theta(M_{q}-M).$ (16)

Here, the quantity

$\displaystyle
D_{13}(\sigma^{\prime}_{0})=\left(\gamma_{s}^{(-)}-e_{2}^{(-)}\right)\left(\tilde{h}^{2}+h^{2}_{cm}\right)^{\frac{d}{d+2}}+$
$\displaystyle+e_{0}^{(-)}\tilde{h}\left(\tilde{h}^{2}+h^{2}_{cm}\right)^{\frac{d-2}{2(d+2)}}$
(17)

depends on the solution $\sigma^{\prime}_{0}$ of Eq. (1); $d=3$ is the space
dimension; $v$ is the cell volume; $\Theta(M)$ is the Heaviside function,
which is equal to unity if $M>0$, to zero if $M<0$, and to $1/2$ if $M=0$; the
quantity $P_{a}^{(-)}(T)$ depends analytically on the temperature; the
coefficient $\gamma_{s}^{(-)}$ characterizes the non-analytical contribution
to the thermodynamic potential; the quantities $e_{0}^{(-)}$ and $e_{2}^{(-)}$
depend on the roots of the cubic equation (1). Expressions for all those
quantities, as well as for $E_{\mu}$, are given in work [10].

The equation of state (16) makes it possible to study the dependence of the
pressure $P$ on the chemical potential $M$ and the relative temperature
$\tau$. This equation can be rewritten in terms of the temperature and density
variables. For this purpose, the chemical potential expressed via the
temperature $\tau$ and the average density $\bar{n}$ should be substituted
into Eq. (16), and the intervals of chemical potential values in the Heaviside
functions have to be replaced by the corresponding density values from the
regions shown in Fig. 3. Let us consider each of the regions separately.

Region I ($M\leq-M_{q}$). Here the solution $\sigma^{\prime}_{0}$ of Eq. (1)
looks like $\sigma^{\prime}_{0b}$ from Eq. (3). If $M=-M_{q}$, then expression
(10), where $\sigma_{0r}$ is given by the relationship from Eq. (7), is valid
for $\sigma^{\prime}_{0b}$. Solution (10) for $\sigma_{0bq}^{{}^{\prime}(-)}$
coincides with $\sigma_{03}^{{}^{\prime}(-)}$ from Eq. (12). On the other
hand, the equality [10]

$b_{3}^{(-)}M^{1/5}=\bar{n}-n_{g}+M,$ (18)

holds, which couples the average density $\bar{n}$ with the chemical potential
(in this case, $M=-M_{q}$) and the quantity
$\sigma_{00}^{(-)}=f(\sigma^{\prime}_{0b})$, which is included into the
coefficient $b_{3}^{(-)}$. Note that $n_{g}$ is determined via the
coefficients in the initial expression for the grand partition function, and
$\sigma_{00}^{(-)}$ is a function of the quantity
$\alpha_{m}=\tilde{h}/h_{cm}$, which includes the initial chemical potential
$\mu$ (included into $M$) and the relative temperature $\tau$. From Eq. (18),
we can determine the density $\bar{n}_{12}$ (the boundary density between
regions I and II) corresponding to the value $M=-M_{q}$. Neglecting the last
term in Eq. (18), we obtain

$\displaystyle\bar{n}_{12}=n_{g}+b_{3}^{(-)}M^{1/5}\big{|}_{\begin{subarray}{c}M=-M_{q}\end{subarray}}=n_{g}+\left[(1+\alpha^{2}_{mq})^{1/2}h_{cm}\right]^{1/5}\sigma_{00}^{(-)}(\sigma_{03}^{{}^{\prime}(-)}),$
$\displaystyle\sigma_{03}^{{}^{\prime}(-)}=-2\sigma_{0r}.$ (19)

Region II ($-M_{q}<M\leq-0$). At $M=-0$, the equality (14) holds true and the
boundary density $\bar{n}_{20}=\lim_{M\rightarrow-0}\bar{n}$ takes the form

$\displaystyle\bar{n}_{20}=n_{g}+\lim_{M\rightarrow-0}\left[(1+\alpha^{2}_{m})^{1/2}h_{cm}\right]^{1/5}\sigma_{00}^{(-)}=n_{g}+h_{cm}^{1/5}\sigma_{00}^{(-)}(\sigma_{030}^{{}^{\prime}(-)}),$
$\displaystyle\sigma_{030}^{{}^{\prime}(-)}=-\sqrt{3}\sigma_{0r}.$ (20)

Region III ($+0\leq M<M_{q}$). This region starts from the value $M=+0$, where
$\sigma_{010}^{{}^{\prime}(+)}=\sqrt{3}\sigma_{0r}$ and, accordingly,

$\bar{n}_{03}=n_{g}+h_{cm}^{1/5}\sigma_{00}^{(-)}(\sigma_{010}^{{}^{\prime}(+)}).$
(21)

The chemical potential $M$ in region III acquires values less than $M_{q}$.

Region IV ($M\geq M_{q}$). This region starts from the value $M=M_{q}$, which
corresponds to $\sigma_{01}^{{}^{\prime}(+)}=2\sigma_{0r}$, so that

$\bar{n}_{34}=n_{g}+\left[(1+\alpha^{2}_{mq})^{1/2}h_{cm}\right]^{1/5}\sigma_{00}^{(-)}(\sigma_{01}^{{}^{\prime}(+)}).$
(22)

The initial boundary density $\bar{n}_{34}$ in region IV increases to a
certain value $\bar{n}_{max}$, which corresponds to $M_{max}$. At
$\bar{n}>\bar{n}_{max}$, the chemical potential $M$ decreases with the
increasing density $\bar{n}$, which does not reflect the physical nature of
the phenomenon (a similar picture is observed at $\bar{n}<\bar{n}_{min}$).

The determination of the boundary densities $\bar{n}_{12}$, $\bar{n}_{20}$,
$\bar{n}_{03}$, and $\bar{n}_{34}$ makes it possible to write the equation of
state (16) in the form

$\displaystyle\frac{Pv}{kT}=P_{a}^{(-)}(T)+E_{\mu}+D_{13}(\sigma^{\prime}_{0b})\left[\Theta(\bar{n}_{12}-\bar{n})+\right.$
$\displaystyle\left.+\Theta(\bar{n}-\bar{n}_{34})\right]+D_{13}(\sigma^{\prime}_{03})\Theta(\bar{n}-\bar{n}_{12})\Theta(-\bar{n}+\bar{n}_{20})+$
$\displaystyle+D_{13}(\sigma^{\prime}_{01})\Theta(\bar{n}-\bar{n}_{03})\Theta(\bar{n}_{34}-\bar{n}),$
(23)

where

$D_{13}(\sigma^{\prime}_{0})\\!=\\!\left(\frac{\bar{n}\\!-\\!n_{g}}{\sigma_{00}^{(-)}}\right)^{6}\\!\\!\left[\\!e_{0}^{(-)}\\!\frac{\alpha_{m}}{(1+\alpha_{m}^{2})^{1/2}}\\!+\\!\gamma_{s}^{(-)}\\!-\\!e_{2}^{(-)}\\!\right].$
(24)

The dependence of the pressure $P$ (see Eq. (23)) on $\bar{n}$ for various
$\tau$ is shown in Fig. 4.

Figure 4: Pressure as a function of average density for various values of the
relative temperature.

##  3 Relationship between the density and the chemical potential of fluid.
Limiting cases

Nonlinear equation (18), which describes the relationship between the density
$\bar{n}$ and the chemical potential $M$, can be rewritten in the form [10]

$\bar{n}=n_{g}-M+\sigma_{00}^{(-)}\left(\tilde{h}^{2}+h^{2}_{cm}\right)^{\frac{d-2}{2(d+2)}}.$
(25)

The general form of equation (25) or (18) makes it possible to change in a
natural way to the cases when either of the variables (the temperature or the
chemical potential) is crucial for the description of the order parameter
behavior.

Let us describe the behavior of $\bar{n}$ for some limiting cases. One of them
is the absence of chemical potential $M$ (i.e., $M=0$ and hence $\tilde{h}=0$)
and $T\neq T_{c}$. Then, we obtain

$\sigma_{00}^{(-)}(M=0)=\frac{e_{0}^{(-)}}{(\beta W(0))^{1/2}}+e_{020}^{(-)},$
(26)

where

$e_{020}^{(-)}=e_{02}^{(-)}(M=0)=\frac{1}{(\beta
W(0))^{1/2}}\frac{f_{Iv}}{s^{3}}.$

The expression for $f_{Iv}$ is given in work [10]. From Eq. (25), we obtain
the dependence

$\bar{n}=n_{g}+\sigma_{00}^{(-)}(M=0)\tilde{\tau}_{1}^{\beta},$ (27)

where the critical exponent $\beta=\nu/2$.

Another limiting case is $M\neq 0$ and $T=T_{c}$. The density $\bar{n}$ in Eq.
(25) at $T=T_{c}$ satisfies the equality

$\bar{n}=n_{g}-M+\sigma_{00}^{(-)}(T_{c})\tilde{h}^{1/\delta},$ (28)

where

$\sigma_{00}^{(-)}(T_{c})=\frac{6}{5}\frac{1}{(\beta_{c}W(0))^{1/2}}\left[e_{0}^{(-)}+\gamma_{s}^{(-)}-e_{2}^{(-)}\right]$
(29)

and the critical exponent $\delta=5$.

In the general case, i.e., at $M\neq 0$ and $T\neq T_{c}$, Eq. (25) can be
presented as follows:

$\bar{n}=n_{g}-M+\sigma_{00}^{(-)}\left(\tilde{h}^{2}+\tilde{\tau}_{1}^{2\beta\delta}\right)^{1/(2\delta)}.$
(30)

Note that $M\ll 1$, and $\tilde{h}\sim M$. Therefore, the second summand, $M$,
in the right-hand sides of Eqs. (25), (28), and (30) is much smaller than the
third term and can be neglected.

##  4 Binodal equation

The binodal equation can be obtained from Eq. (25) by putting $M=0$. Then we
arrive at Eq. (27). Now, substituting the expression
$\tilde{\tau}_{1}=-\tau\frac{c_{11}}{q}E_{2}^{n_{0}}$, we obtain

$\bar{n}=n_{g}+\sigma_{00}^{(-)}(M=0)\left(-\tau\frac{c_{11}}{q}E_{2}^{n_{0}}\right)^{\beta}.$
(31)

Here, $E_{2}$ is one of the eigenvalues of the matrix for the linear
transformation of the renormalization group, the quantity $c_{11}$
characterizes one of the coefficients in the solutions of recurrence relations
for the $\rho^{4}$-model [10, 9], $n_{0}$ is the difference between the exit
points from the critical fluctuation regime at $T>T_{c}$ and $T<T_{c}$, and
$q$ is associated with the averaging of the wave vector square. For more
information on those parameters, see work [10].

Let us solve Eq. (31) with respect to the temperature. Taking the equalities
$\tau=T/T_{c}-1$ and $\beta=\nu/2$ into account, we obtain

$\left[\frac{\left(\frac{\bar{n}}{n_{g}}-1\right)n_{g}}{\sigma_{00}^{(-)}(M=0)}\right]^{2/\nu}\frac{q}{c_{11}E_{2}^{n_{0}}}=-\frac{T}{T_{c}}+1$
(32)

or

$\frac{T}{T_{c}}=1-\left\\{\left[\frac{\left(\frac{\bar{n}}{n_{g}}-1\right)n_{g}}{\sigma_{00}^{(-)}(M=0)}\right]^{2}\right\\}^{1/\nu}\frac{q}{c_{11}E_{2}^{n_{0}}}.$
(33)

On the basis of Eq. (33), we can plot the binodal curve in the temperature
versus density plane (see Fig. 5).

Figure 5: Coexistence curve (binodal curve) obtained in an immediate vicinity
of the critical point taking into account the interaction potential parameters
that are characteristic of sodium. The solid curve (the dome) was plotted
according to the obtained binodal equation, and the dotted curve is a result
of the zero-mode approximation [12].

This curve agrees with the data predicted for sodium by extrapolating the
results of computer simulation [15] to $T/T_{c}\approx 1$ (see work [10]).

The spinodal equation, which describes the limiting states of the system that
determine the boundaries of the instability region, can be found from the
extremum condition

$\frac{\partial(Pv/kT)}{\partial\bar{n}}\bigg{|}_{T}=0$

for the equation of state (23), where we should substitute the chemical
potential $M$ expressed from Eq. (18) in terms of the average density
$\bar{n}$.

##  5 Conclusions

In this paper, the cell model was used to study the behavior of fluid in a
close vicinity of its critical point. This is an interesting (because of
fundamental and applied aspects) and difficult (because of a substantial role
of fluctuation effects) issue for analysis. The study of the relationship
between the density and the chemical potential at temperatures $T<T_{c}$ made
it possible to determine the corresponding densities in the regions where the
chemical potential changes and obtain both the equation of state and the
binodal equation. Solutions of a certain cubic equation were presented, which
the equation of state of the cell fluid model depends on. Their analysis made
it possible to describe the transition from one solution to another when the
chemical potential tends to zero. On the basis of the obtained equation of
state, the pressure variation with the density growth at various temperatures
has been illustrated graphically. Using the binodal equation, a binodal curve
in a narrow temperature interval was constructed for the microscopic
parameters of the Morse potential that are inherent to sodium. As compared
with the case of the zero-mode approximation, the obtained dome of the
coexistence curve is wider and agrees better with the results of computer
simulation [15].

## References

  * [1] C.-L. Lee, G. Stell, J.S. Høye. A simple SCOZA for simple fluids. J. Mol. Liq. 112, 13 (2004).
  * [2] C.E. Bertrand, J.F. Nicoll, M.A. Anisimov. Comparison of complete scaling and a field-theoretic treatment of asymmetric fluid criticality. Phys. Rev. E 85, 031131 (2012).
  * [3] A. Parola, L. Reatto. Recent developments of the hierarchical reference theory of fluids and its relation to the renormalization group. Mol. Phys. 110, 2859 (2012).
  * [4] I.R. Yukhnovskii. The phase transition of the first order in the critical region of the gas-liquid system. Condens. Matter Phys. 17, 43001 (2014).
  * [5] T.J. Yoon, Y.-W. Lee. Current theoretical opinions and perspectives on the fundamental description of supercritical fluids. J. Supercrit. Fluids 134, 21 (2018).
  * [6] L.F. Vega. Perspectives on molecular modeling of supercritical fluids: From equations of state to molecular simulations. Recent advances, remaining challenges and opportunities. J. Supercrit. Fluids 134, 41 (2018).
  * [7] A. Oleinikova, L. Bulavin, V. Pipich. Critical anomaly of shear viscosity in a mixture with an ionic impurity. Chem. Phys. Lett. 278, 121 (1997).
  * [8] S. Pittois, B. Van Roie, C. Glorieux, J. Thoen. Thermal conductivity, thermal effusivity, and specific heat capacity near the lower critical point of the binary liquid mixture n-butoxyethanol–water. J.Chem. Phys. 121, 1866 (2004).
  * [9] M.P. Kozlovskii, I.V. Pylyuk, O.A. Dobush. The equation of state of a cell fluid model in the supercritical region. Condens. Matter Phys. 21, 43502 (2018).
  * [10] I.V. Pylyuk. Fluid critical behavior at liquid–gas phase transition: Analytic method for microscopic description. J. Mol. Liq. 310, 112933 (2020).
  * [11] M. Kozlovskii, O. Dobush. Representation of the grand partition function of the cell model: The state equation in the mean-field approximation. J. Mol. Liq. 215, 58 (2016).
  * [12] M.P. Kozlovskii, O.A. Dobush, I.V. Pylyuk. Using a cell fluid model for the description of a phase transition in simple liquid alkali metals. Ukr. J. Phys. 62, 865 (2017).
  * [13] I.V. Pylyuk, O.A. Dobush. Equation of state of a cell fluid model with allowance for Gaussian fluctuations of the order parameter. Ukr. J. Phys. 65, 1080 (2020).
  * [14] I.R. Yukhnovskii. Phase Transitions of the Second Order. Collective Variables Method (World Scientific, 1987).
  * [15] J.K. Singh, J. Adhikari, S.K. Kwak. Vapor–liquid phase coexistence curves for Morse fluids. Fluid Phase Equilib. 248, 1 (2006).
  * [16] E.M. Apfelbaum. The calculation of vapor-liquid coexistence curve of Morse fluid: Application to iron. J. Chem. Phys. 134, 194506 (2011).
  * [17] Y. Zhou, M. Karplus, K.D. Ball, R.S. Berry. The distance fluctuation criterion for melting: Comparison of square-well and Morse potential models for clusters and homopolymers. J. Chem. Phys. 116, 2323 (2002).
  * [18] X. Xu, C. Cheng, I.H. Chowdhury. Molecular dynamics study of phase change mechanisms during femtosecond laser ablation. J. Heat Transfer 126, 727 (2004).
  * [19] I. Last, Y. Levy, J. Jortner. Beyond the Rayleigh instability limit for multicharged finite systems: From fission to Coulomb explosion. Proc. Natl. Acad. Sci. USA 99, 9107 (2002).
  * [20] J.P.K. Doye, D.J. Wales. The structure and stability of atomic liquids: From clusters to bulk. Science 271, 484 (1996).
  * [21] J.P.K. Doye, R.H. Leary, M. Locatelli, F. Schoen. Global optimization of Morse clusters by potential energy transformations. INFORMS J. Comput. 16, 371 (2004).
  * [22] A. Tekin, M. Yurtsever. Molecular dynamics simulation of phase transitions in binary LJ clusters. Turk. J. Chem. 26, 627 (2002).
  * [23] C.-I. Chou, C.-L. Ho, B. Hu, H. Lee. Morse-type Frenkel-Kontorova model. Phys. Rev. E 57, 2747 (1998).
  * [24] A. Strachan, T. Cagin, W.A. Goddard, III. Phase diagram of MgO from density-functional theory and molecular-dynamics simulations. Phys. Rev. B 60, 15084 (1999).
  * [25] T.-C. Lim. Approximate relationships between the Generalized Morse and the Extended-Rydberg potential energy functions. Acta Chim. Slov. 52, 149 (2005).
  * [26] A. Del Sol Mesa, C. Quesne, Yu.F. Smirnov. Generalized Morse potential: Symmetry and satellite potentials. J. Phys. A: Math. Gen. 31, 321 (1998).
  * [27] P.M. Morse, E.C.G. Stueckelberg. Diatomic molecules according to the wave mechanics I: Electronic levels of the hydrogen molecular ion. Phys. Rev. 33, 932 (1929).
  * [28] A.S. Leal, C. Gouvea dos Santos, C.M. Quintella, H.H.R. Schor. A theoretical model for the scattering of $I_{2}$ molecule from a perfluoropolyeter liquid surface. J. Braz. Chem. Soc. 10, 359 (1999).
  * [29] V. Constantoudis, C.A. Nicolaides. Stabilization and relative phase effects in a dichromatically driven diatomic Morse molecule: Interpretation based on nonlinear classical dynamics. J. Chem. Phys. 122, 084118 (2005).
  * [30] A.I. Milchev, A.A. Milchev. Wetting behavior of nanodroplets: The limits of Young’s rule validity. Europhys. Lett. 56, 695 (2001).
  * [31] D. Osorio-Gonzalez, M. Mayorga, J. Orozco, L. Romero-Salazar. Entropy and thermalization of particles in liquids. J. Chem. Phys. 118, 6989 (2003).
  * [32] P. Shah, C. Chakravarty. Instantaneous normal mode analysis of Morse liquids. J. Chem. Phys. 116, 10825 (2002).
  * [33] H. Okumura, F. Yonezawa. Liquid–vapor coexistence curves of several interatomic model potentials. J. Chem. Phys. 113, 9162 (2000).
  * [34] T.-C. Lim. The relationship between Lennard-Jones (12-6) and Morse potential functions. Z. Naturforsch. 58a, 615 (2003).