# A new discrete calculus of variations and its applications in statistical
physics

Q. H. Liu<EMAIL_ADDRESS>School for Theoretical Physics, School of
Physics and Electronics, Hunan University, Changsha 410082, China

###### Abstract

For a discrete function $f\left(x\right)$ on a discrete set, the finite
difference can be either forward and backward. However, we observe that if
$f\left(x\right)$ is a sum of two functions
$f\left(x\right)=f_{1}\left(x\right)+f_{2}\left(x\right)$ defined on the
discrete set, the first order difference of $\Delta f\left(x\right)$ is
equivocal for we may have
$\Delta^{f}f_{1}\left(x\right)+\Delta^{b}f_{2}\left(x\right)$ where
$\Delta^{f}$ and $\Delta^{b}$ denotes the forward and backward difference
respectively. Thus, the first order variation equation for this function
$f\left(x\right)$ gives many solutions which include both true and false one.
A proper formalism of the discrete calculus of variations is proposed to
single out the true one by examination of the second order variations, and is
capable of yielding the exact form of the distributions for Boltzmann, Bose
and Fermi system without requiring the numbers of particle to be infinitely
large. The advantage and peculiarity of our formalism are explicitly
illustrated by the derivation of the Bose distribution.

statistical distribution, discrete calculus of variations, most probable
distribution

_Introduction_ The method of _most probable distribution_(MPD)__ is the most
common way used to derive various statistical distributions in mathematics,
physics, chemistry, materials science and computational science, etc. The
original data are usually discrete rather than continuous, so we must be able
to deal with the differences, difference quotients and sums of discrete
functions, instead of the differentials, derivatives and integrations of the
continuous functions. However, we immediately run into a problem: For a
discrete function that has extremals, the difference can be forward, backward,
central, and even more complicated combinations of these differences, so that
the first order variation equation leads to many solutions which can not be
all true. By the true solution, we mean that it is realistic or physical.
There has been no satisfactory procedure to single out the true one in its own
right. Thus the derivation of the exact form of the distribution functions for
Boltzmann, Bose and Fermi system obtained from the MPD has been an open
problem for nearly one century. chem ; PLA  In present work, I report a
discrete calculus of variations to resolve this problem.

_A new discrete calculus of variations_ Let us first define _an_ _individual
function_ (IF) $f\left(x\right)$ of variable $x$ on an interval $C$, which
amounts to that the exponential of the function $f\left(x\right)$ as $\exp
f\left(x\right)$ can not be further factorized. We can therefore define the
sum of two IFs $f\left(x\right)=f_{1}\left(x\right)+f_{2}\left(x\right)$, and
the sum of many IFs. Secondly we define the _unidirectional differences_ of
the IF __ $f\left(x\right)$, which means that when taking the finite
differences of different orders for the function __ $f\left(x\right)$, we take
_forward or backward_ differences in all first, second, and higher order
variation. For instance, the _forward directional differences_ of the IF __
$f\left(x\right)$ are
$\left(\Delta^{f}\right)f\left(x\right)=f\left(x+h\right)-f\left(x\right)$,
($h\succ 0$),
$\left(\Delta^{f}\right)^{2}f\left(x\right)=f\left(x+2h\right)-f\left(x+h\right)-\left(f\left(x+h\right)-f\left(x\right)\right)=f\left(x+2h\right)-2f\left(x+h\right)+f\left(x\right)$,
and so forth; and the _backward unidirectional differences_ of the IF __
$f\left(x\right)$ are
$\left(\Delta^{b}\right)f\left(x\right)=f\left(x\right)-f\left(x-k\right)$,
($k\succ 0$),
$\left(\Delta^{b}\right)^{2}f\left(x\right)=f\left(x\right)-f\left(x-k\right)-\left(f\left(x-k\right)-f\left(x-2k\right)\right)=f\left(x\right)-2f\left(x-k\right)+f\left(x-2k\right)$,
and so forth. Thus, the _forward directional difference quotients_ are
$\left(\Delta^{f}\right)f\left(x\right)/h$,
$\left(\Delta^{f}\right)^{2}f\left(x\right)/h^{2}$, …; and similarly, the
_backward directional difference quotients_ are
$\left(\Delta^{b}\right)f\left(x\right)/k$,
$\left(\Delta^{b}\right)^{2}f\left(x\right)/h^{2}$, …. Specially for $h=1$, we
have
$\left(\Delta^{f}\right)^{i}f\left(x\right)/h=\left(\Delta^{f}\right)^{i}f\left(x\right)$,
($i=1,2,3,...$); and similarly for $k=1$, we have
$\left(\Delta^{b}\right)^{i}f\left(x\right)/k=\left(\Delta^{b}\right)^{i}f\left(x\right)$.
The _unidirectional differences_ are similar to the treatment of the parameter
$\lambda\in\mathbb{R}$ in the usual calculus of variation in which for the
given variational parameter $\lambda$, we must keep using this parameter
$\lambda$ in all orders of variations. arnold  However, the different point of
the new procedure is that, once
$f\left(x\right)=f_{1}\left(x\right)+f_{2}\left(x\right)$, $\Delta f$ does
mean four combinations:
$\left(\Delta^{f}\right)f_{1}\left(x\right)+\left(\Delta^{f}\right)f_{2}\left(x\right)$,
$\left(\Delta^{f}\right)f_{1}\left(x\right)+\left(\Delta^{b}\right)f_{2}\left(x\right)$,
$\left(\Delta^{b}\right)f_{1}\left(x\right)+\left(\Delta^{f}\right)f_{2}\left(x\right)$,
and
$\left(\Delta^{b}\right)f_{1}\left(x\right)+\left(\Delta^{b}\right)f_{2}\left(x\right)$,
which can be simply named by 1f2f, 1f2b, 1b2f and 1b2b, respectively.
Moreover, the second order variations are must be taken _unidirectionally_ ,
so $\Delta^{2}f$ for the combination 1f2f, 1f2b, 1b2f and 1b2b gives,
$\left(\Delta^{f}\right)^{2}f_{1}\left(x\right)+\left(\Delta^{f}\right)^{2}f_{2}\left(x\right)$,
$\left(\Delta^{f}\right)^{2}f_{1}\left(x\right)+\left(\Delta^{b}\right)^{2}f_{2}\left(x\right)$,
$\left(\Delta^{b}\right)^{2}f_{1}\left(x\right)+\left(\Delta^{f}\right)^{2}f_{2}\left(x\right)$,
and
$\left(\Delta^{b}\right)^{2}f_{1}\left(x\right)+\left(\Delta^{b}\right)^{2}f_{2}\left(x\right)$,
respectively.

The problem is formulated in the following. We at first deal with a discrete
IF $\Psi\left(n\right)$ defined, for convenience, on the interval of semi-
positive integers $n\in\mathbb{Z}^{+}$, which has the local maxima and minima,
subject to some equality constraints
$\mathbf{\psi=}\left(\psi_{1},\psi_{2},\psi_{3},...\right)=\mathbf{0}$. For
finding the local maxima and minima of the function $\Psi\left(n\right)$, we
construct a functional $\Phi$,

$\Phi=\Psi\left\\{n\right\\}+\mathbf{\alpha\cdot\psi,}$ (1)

where $\mathbf{\alpha=}\left(\alpha_{1},\alpha_{2},\alpha_{3},...\right)$ are
Lagrange multipliers each of which $\alpha_{i}$ companies a constraint
condition $\psi_{i}=0$. The local maxima and minima satisfy,

$\delta\Phi=0.$ (2)

Since the discreteness of the function $\Psi\left(n\right)$, the smallest
finite change of $n$ is $\Delta n=1$, and we have accordingly two differences
(or difference quotients) of the function $\Psi\left\\{n\right\\}$ in the
following: $\left(\delta^{f}\right)\Psi\left\\{n\right\\}/\delta
n=\Psi\left\\{n+1\right\\}-\Psi\left\\{n\right\\}$
$=\left(\Delta^{f}\right)\Psi\left\\{n\right\\}$ for the forward difference
and $\left(\delta^{b}\right)\Psi\left\\{n\right\\}/\delta
n=\Psi\left\\{n\right\\}-\Psi\left\\{n-1\right\\}=\left(\Delta^{b}\right)\Psi\left\\{n\right\\}$
for the backward difference. In consequence, each of two relations
$n^{\mu}=n^{\mu}\left(\mathbf{\alpha}\right)$ ($\mu=1,2$) solves the equation
(2), constituting the solution pool which contains both spurious and true one.
Then how to single out the true one?

Assuming that all relations $n^{\mu}=n^{\mu}\left(\mathbf{\alpha}\right)$ in
the solution pool have different values of the second order differences in the
sense of the _unidirectional differences_ of $\Psi\left\\{n\right\\}$, and
further assuming that $n^{1}=n^{1}\left(\mathbf{\alpha}\right)$ is obtained
from the forward difference variation as
$\left(\Delta^{f}\right)\Psi\left\\{n^{1}\right\\}+\mathbf{\alpha\cdot}\Delta^{f}\mathbf{\psi}\left\\{n^{1}\right\\}=0$,
the another solution $n^{2}=n^{2}\left(\mathbf{\alpha}\right)$ must then be
obtained from the backward difference variation as
$\left(\Delta^{b}\right)\Psi\left\\{n^{2}\right\\}+\mathbf{\alpha\cdot}\Delta^{b}\mathbf{\psi}\left\\{n^{2}\right\\}=0$.
To discriminate the difference between $n^{1}\left(\mathbf{\alpha}\right)\ $or
$n^{2}\left(\mathbf{\alpha}\right)$, we examine
$\left(\Delta^{f}\right)^{2}\Psi\left\\{n^{1}\right\\}$ and
$\left(\Delta^{b}\right)^{2}\Psi\left\\{n^{2}\right\\}\left(\neq\left(\Delta^{f}\right)^{2}\Psi\left\\{n^{1}\right\\}\right)$.
Thus, it is reasonable to conjecture that if
$\left|\left(\Delta^{f}\right)^{2}\Psi\left\\{n^{1}\right\\}\right|>\left|\left(\Delta^{b}\right)^{2}\Psi\left\\{n^{2}\right\\}\right|$,
the solution $n^{1}=n^{1}\left(\mathbf{\alpha}\right)$ is true, and vice
versa. In other words, the second order difference of the true solution take
largest value as
$\max\left\\{\left|\Delta^{2}\Psi\left\\{n^{1}\right\\}\right|,\left|\Delta^{2}\Psi\left\\{n^{2}\right\\}\right|\right\\}$,
where $\Delta^{2}$ must be used in the sense of the _unidirectional
differences_. This is another point of the new procedure greatly different
from the conventional one that accepts all solutions once they are stationary.

If the _discrete function_ $\Psi\left(n\right)$ is the sum of two IFs __
$\Psi_{1}\left(n\right)$ and $\Psi_{2}\left(n\right)$,
$\Psi\left(n\right)=\Psi_{1}\left(n\right)+\Psi_{2}\left(n\right)$, there are
four relations $n^{\xi}=n^{\xi}\left(\mathbf{\alpha}\right)$ ($\xi=1,2,3,4$)
from (2) in the solution pool, satisfying, respectively,

$\displaystyle\left(\Delta^{f}\right)\Psi_{1}\left(n^{1}\right)+\left(\Delta^{f}\right)\Psi_{2}\left(n^{1}\right)+\mathbf{\alpha\cdot}\Delta\mathbf{\psi}\left\\{n^{1}\right\\}$
$\displaystyle=$ $\displaystyle 0,\text{ for 1f2f,}$ (3a)
$\displaystyle\left(\Delta^{f}\right)\Psi_{1}\left(n^{2}\right)+\left(\Delta^{b}\right)\Psi_{2}\left(n^{2}\right)+\mathbf{\alpha\cdot}\Delta\mathbf{\psi}\left\\{n^{2}\right\\}$
$\displaystyle=$ $\displaystyle 0,\text{ for 1f2b,}$ (3b)
$\displaystyle\left(\Delta^{b}\right)\Psi_{1}\left(n^{3}\right)+\left(\Delta^{f}\right)\Psi_{2}\left(n^{3}\right)+\mathbf{\alpha\cdot}\Delta\mathbf{\psi}\left\\{n^{3}\right\\}$
$\displaystyle=$ $\displaystyle 0,\text{ for 1b2f,}$ (3c)
$\displaystyle\left(\Delta^{b}\right)\Psi_{1}\left(n^{4}\right)+\left(\Delta^{b}\right)\Psi_{2}\left(n^{4}\right)+\mathbf{\alpha\cdot}\Delta\mathbf{\psi}\left\\{n^{4}\right\\}$
$\displaystyle=$ $\displaystyle 0,\text{ for 1b2b,}$ (3d) where the simple
form of the constraint functions $\mathbf{\psi}$ all proportional to $n$ is
assumed, for simplicity, so that
$\Delta\mathbf{\psi}\left\\{n^{\mu}\right\\}=\left(\Delta^{f}\right)\mathbf{\psi}\left\\{n\right\\}=\left(\Delta^{b}\right)\mathbf{\psi}\left\\{n\right\\}$.
The mixing form of the solution 1f2b (3b) and 1b2f (3c) can not be excluded
but has not been noted before. It is understandable, for a continuous function
$f\left(x\right)=f_{1}\left(x\right)+f_{2}\left(x\right)$, no difference is
possible in differentials ($df=df_{1}+$ $df_{2}$) or derivatives
($df/dx=df_{1}/dx+df_{2}/dx$), when these operations acting on two functions
$f_{1}\left(x\right)$ and $f_{2}\left(x\right)$; and it is then reasonably
infer that when $f\left(x\right)=f_{1}\left(x\right)+f_{2}\left(x\right)$, we
have still $\Delta f\left(x\right)=\Delta f_{1}\left(x\right)+\Delta
f_{2}\left(x\right)$ without mixing of $\Delta^{f}$ and $\Delta^{b}$. The new
procedure conjectures that the true solution
$n^{\xi}=n^{\xi}\left(\mathbf{\alpha}\right)$ is that its second order
difference $\Delta^{2}\Psi\left(n^{\xi}\right)$ takes the largest value in
magnitude among all solutions in the solution pool,
$\left|\Delta^{2}\Psi\left(n^{\xi}\right)\right|=\max\left\\{\left|\Delta^{2}\Psi\left\\{n^{1}\right\\}\right|\text{,}\left|\Delta^{2}\Psi\left\\{n^{2}\right\\}\right|\text{,}\left|\Delta^{2}\Psi\left\\{n^{3}\right\\}\right|\text{,}\left|\Delta^{2}\Psi\left\\{n^{4}\right\\}\right|\right\\}.$
(4)

Five immediate comments follow. 1. The new procedure is essentially a method
of MPD, and importantly during all derivation steps, no one requires that the
variable $n$ be very large. 2. In the new procedure, the correct solution and
the true solution differ; and the former refers to its solving the first order
equation (2) and the latter refers to its maximizing the second order
variations. 3. Once there are degenerate solutions which have no difference up
to second differences, higher order differences must be invoked. 4. The new
procedure can be easily generalized for the _discrete function_
$\Psi\left(n\right)$ that is the sum of more IFs. 5. If the constraint
function $\mathbf{\psi}$ is nonlinear in $n$, which is beyond the scope of
current studies, the problem must be treated on the case-by-case base.

We will simply call the new procedure of the discrete calculus of variations
the _IF method_. In the following, I will illustrate the IF method with a
detailed derivation of the exact form of the Bose distribution, and slightly
discuss the Boltzmann and Fermi distribution.

_IF method for the Bose system_ Considering a system of $N$ noninteracting,
indistinguishable particles confined to a space of volume $V$ and sharing a
given energy $E$. Let $\varepsilon_{i}$ denote the energy of $i$-th level and
$\varepsilon_{1}\prec\varepsilon_{2}\prec\varepsilon_{3}\prec...$, and $g_{i}$
denote the degeneracy of the level. In a particular situation, we may have
$n_{1}$ particles in the first level $\varepsilon_{1}$, $n_{2}$ particles in
the second level $\varepsilon_{2}$, and so on, defining a distribution set
$\left\\{n_{i}\right\\}$ phys . The number of the distinct microstates in set
$\left\\{n_{i}\right\\}$ is then given by,

$\Omega\left\\{n_{i}\right\\}=\prod_{i}\frac{(n_{i}+g_{i}-1)!}{n_{i}!(g_{i}-1)!}.$
(5)

The Bose functional $f$ is,

$f=\sum_{i}\left(\ln(n_{i}+g_{i}-1)!-\ln
n_{i}!-\ln(g_{i}-1)!\right)-\alpha\left(\sum_{i}n_{i}-N\right)-\beta\left(\sum_{i}n_{i}\varepsilon_{i}-E\right).$
(6)

The variational $\delta f$ is,

$\displaystyle\delta f$ $\displaystyle=$
$\displaystyle\sum_{i}\left\\{\delta\ln(n_{i}+g_{i}-1)!-\delta\ln
n_{i}!\right\\}-\delta n_{i}\left(\alpha+\beta\varepsilon_{i}\right)$ (7)
$\displaystyle=$ $\displaystyle\sum_{i}\delta
n_{i}\left\\{\frac{\delta\ln(n_{i}+g_{i}-1)!}{\delta n_{i}}-\frac{\delta\ln
n_{i}!}{\delta n_{i}}-\left(\alpha+\beta\varepsilon_{i}\right)\right\\}.$

Since the independence of variables $n_{i}$, $\delta f=0$ leads to,

$\frac{\delta\ln(n_{i}+g_{i}-1)!}{\delta n_{i}}-\frac{\delta\ln n_{i}!}{\delta
n_{i}}-\left(\alpha+\beta\varepsilon_{i}\right)=0.$ (8)

Thus, there are essentially _two IFs_ $\Psi_{1}=\ln(n_{i}+g_{i}-1)!$ and
$\Psi_{2}=\ln n_{i}!$_._ For $\ln(n_{i}+g_{i}-1)!$ and $\ln n_{i}!$, the
smallest forward and backward differences are given in Table (1). Four
combinations constructing $\delta\ln(n_{i}+g_{i}-1)!/\delta n_{i}-\delta\ln
n_{i}^{\ast}!/\delta n_{i}$ are summarized in Table (2). Accordingly, we have
four solutions presented in Table (3), and all satisfy $\delta f=0$. These
distributions are mutually different when $n_{i}\sim 1$, though all of them
converge to the same one when $n_{i}\gg 1$,

$n_{i}\approx\frac{g_{i}}{e^{\alpha+\beta\varepsilon_{i}}-1}.$ (9)

Table 1: The smallest forward/backward difference of $\Psi_{1}=\ln(n_{i}+g_{i}-1)!$ and $\Psi_{2}=\ln n_{i}!$. | $\Psi_{1}=\ln(n_{i}+g_{i}-1)!$ | $\Psi_{2}=\ln n_{i}!$
---|---|---
forward | $\ln(n_{i}+g_{i})$ | $\ln(n_{i}+1)$
backward | $\ln(n_{i}+g_{i}-1)$ | $\ln(n_{i})$
Table 2: The column and row give the forward/backward differences for $\Psi_{1}=\ln(n_{i}+g_{i}-1)!$ and $\Psi_{2}=\ln n_{i}!$, respectively, and we list the results of four combinations $\Delta\ln(n_{i}+g_{i}-1)!-\Delta\ln n_{i}!$. |  $\Psi_{2}$ forward |  $\Psi_{2}$ backward
---|---|---
$\Psi_{1}$ forward |  $\ln(n_{i}+g_{i})-\ln(n_{i}+1)$ for 1f2f |  $\ln(n_{i}+g_{i})-\ln n_{i}$ for 1f2b
$\Psi_{1}$ backward |  $\ln(n_{i}+g_{i}-1)-\ln(n_{i}+1)$ for 1b2f |  $\ln(n_{i}+g_{i}-1)-\ln n_{i}$ for 1b2b
Table 3: Four solutions |  $\Psi_{2}$ forward |  $\Psi_{2}$ backward
---|---|---
$\Psi_{1}$ forward |  $\frac{n_{i}+g_{i}}{n_{i}+1}=\exp\left(\alpha+\beta\varepsilon_{i}\right)$ for 1f2f |  $\frac{n_{i}+g_{i}}{n_{i}}=$ $\exp\left(\alpha+\beta\varepsilon_{i}\right)$ for 1f2b
$\Psi_{1}$ backward |  $\frac{n_{i}+g_{i}-1}{n_{i}+1}=$ $\exp\left(\alpha+\beta\varepsilon_{i}\right)$ for 1b2f |  $\frac{n_{i}+g_{i}-1}{n_{i}}=\exp\left(\alpha+\beta\varepsilon_{i}\right)$ for 1b2b

The second order variations of $\delta^{2}\ln\Omega\left\\{n_{i}{}\right\\}$
for the four solutions are explicitly shown in the Table (4).

Table 4: The column and row give the forward/backward finite differences for $\ln(n_{i}+g_{i}-1)!$ and $\ln n_{i}!$, respectively, and then forming $\Delta^{2}\ln(n_{i}+g_{i}-1)!-\Delta^{2}\ln n_{i}!$ accordingly. |  $\Psi_{2}$ forward |  $\Psi_{2}$ backward
---|---|---
$\Psi_{1}$ forward |  $-\ln(\frac{n_{i}+2}{n_{i}+1}\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1})$ for 1f2f |  $-\ln(\frac{n_{i}}{n_{i}-1}\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1})$ for 1f2b
$\Psi_{1}$ backward |  $-\ln(\frac{n_{i}+2}{n_{i}+1}\frac{n_{i}+g_{i}-2}{n_{i}+g_{i}-1})$ for 1b2f |  $-\ln(\frac{n_{i}}{n_{i}-1}\frac{n_{i}+g_{i}-2}{n_{i}+g_{i}-1})$ for 1b2b

In order to eliminate the spurious solutions, we examine which one is the
largest in magnitude. It is easily to verify that one combination 1f2b is the
only right one for we have,

$\displaystyle\frac{n_{i}}{n_{i}-1}\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1}-\frac{n_{i}}{n_{i}-1}\frac{n_{i}+g_{i}-2}{n_{i}+g_{i}-1}$
$\displaystyle=$
$\displaystyle\frac{n_{i}}{n_{i}-1}\left(\frac{2}{\left(n_{i}+g_{i}\right)^{2}-1}\right)>0,\text{
}\left(n_{i}>1\right)$ (10a)
$\displaystyle\frac{n_{i}}{n_{i}-1}\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1}-\frac{n_{i}+2}{n_{i}+1}\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1}$
$\displaystyle=$
$\displaystyle\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1}\left(\frac{2}{n_{i}{}^{2}-1}\right)>0,\text{
}\left(n_{i}>1\right)$ (10b)
$\displaystyle\frac{n_{i}}{n_{i}-1}\frac{n_{i}+g_{i}}{n_{i}+g_{i}+1}-\frac{n_{i}+2}{n_{i}+1}\frac{n_{i}+g_{i}-2}{n_{i}+g_{i}-1}$
$\displaystyle=$ $\displaystyle
2\frac{(2n_{i}{}^{2}+2g_{i}n_{i}-2-g_{i}+g_{i}^{2})}{\left(\left(n_{i}+g_{i}\right)^{2}-1\right)\left(n_{i}{}^{2}-1\right)}>0,\text{
}\left(n_{i}>1\right).$ (10c) To note that once $n_{i}=0$, and $n_{i}=1$, we
need to directly invoke the expression for 1f2b in Table (4) and the
$\left|\ln(n_{i}/(n_{i}-1)\right|\rightarrow\infty$, and in final we obtain
the Bose distribution,
$n_{i}=\frac{g_{i}}{e^{\alpha+\beta\varepsilon_{i}}-1}.$ (11)

Two remarks follow. 1. We can follow the similar manner to present detailed
derivations for Boltzmann and Fermi distribution. For the Boltzmann
distribution, there is essentially one IF and the true solution comes from
backward difference; ; and for the Fermi distribution there are essentially
two IFs and the true solution comes from the combination 1b2b. 2, Our
procedure puts no requirement on the total number of particle $N$ that be very
large. Within the grand ensemble theory, we can easily show the Bose
distribution to be given by with $g_{i}=1$,

$n_{i}=\frac{1}{e^{\alpha+\beta\varepsilon_{i}}-1}+\frac{1+N}{e^{\left(\alpha+\beta\varepsilon_{i}\right)\left(1+N\right)}-1}.$
(12)

It differs from ours (11) by an additional term depending on $N$.

_Conclusions and discussions_ In contrast to the continuous calculus of
variations, the discrete one possesses some peculiarities. The function in it
can not be treated in the uniformly forward or backward difference, but must
exhaust all mathematical possibilities. All possible solutions determined by
the first order variation equation include both true and false one. Each of
the solutions has its own second order variation and the true solution takes
the largest value in magnitude. The application of the procedure to the
statistical distributions is successful, giving results identical to those
obtained from the grand ensemble theory. However, the ensemble theory is valid
only when the total number of particle in the system is large. Our approach
does not suffer from such a limitation.

###### Acknowledgements.

This work is financially supported by National Natural Science Foundation of
China under Grant No. 11675051.

## References

  * (1) P. A. H. Wyatt, Elementary statistical mechanics without Stirling’s approximation, J. Chem. Educ. 39, 27(1962).
  * (2) R. K. Niven, Exact Maxwell-Boltzmann, Bose-Einstein and Fermi-Dirac statistics, Phys. Lett. A., 342, 286(2005).
  * (3) V. I. Arnold, Mathematical methods of classical mechanics, (New York: Springer-Verlag, 1978).
  * (4) L. D. Landau and E. M. Lifshitz, Statistical Physics I, 3rd ed. (Oxford: Pergamon Press, 1980).


# Invariant yield and azimuthal anisotropy measurements of strange and multi-
strange hadrons in Au+Au collision at $\sqrt{s_{NN}}$ = 27 and 54.4 GeV at
STAR

Indian Institute of Science Education and Research, Berhampur, India
E-mail A footnote may follow.

###### Abstract:

In this proceedings, we have presented the invariant yield measurements of
strange and multi strange hadrons such as $K^{0}_{S}$, $\Lambda$,
$\bar{\Lambda}$, $\phi$, $\Xi^{-}$, $\bar{\Xi}^{+}$, $\Omega^{-}$,
$\bar{\Omega}^{+}$ in Au+Au collisions at $\sqrt{s_{NN}}$ = 54.4 GeV at
midrapidity ($|y|<0.5$). The second and third-order azimuthal anisotropic flow
of the above mentioned particles at $\sqrt{s_{NN}}$ = 27 and 54.4 GeV are also
discussed. Some of the important observables such as nuclear modification
factor ($R_{CP}$), baryon to meson ratio ($\Omega/\phi$), transverse momentum
dependence of $v_{2}$ and $v_{3}$ are measured. The number of constituent
quark (NCQ) scaling of $v_{2}$ and $v_{3}$ is studied for all these particles.

##  1 Introduction

Quantum chromodynamics predicts the existence of a deconfined state of quarks
and gluons at high temperature and density. Heavy-ion collisions are a unique
tool to create and study such a state of the matter. There are many
observables to probe the QGP phase of the system created in heavy-ion
collisions, some of them are nuclear modification factor, baryon to meson
ratio, anisotropic flow.

The nuclear modification factor is the ratio of the particle yield in central
to peripheral collisions scaled by the number of binary collisions which is
given by

$R_{CP}=\frac{[\frac{dN}{dp_{T}}\frac{1}{N_{coll}}]_{central}}{[\frac{dN}{dp_{T}}\frac{1}{N_{coll}}]_{peripheral}}.$
(1)

Suppression of the ratio at high $p_{T}$ indicates the energy loss by the high
$p_{T}$ partons inside colour field of the QGP medium. The baryon to meson
ratio, $\Omega/\phi$, gives us information about the hadron formation
mechanism during the hadronization. Another important observable to probe the
QGP phase is anisotropic flow. Flow is defined as the collective expansion of
the medium created in heavy-ion collisions. Depending on the initial spatial
geometry and event-by-event fluctuations in the distribution of nucleons
inside the colliding nuclei there can be the presence of different orders of
anisotropy in the flow. It has been predicted by model calculations that
anisotropic flow coefficients are sensitive to the equation of the state and
shear viscosity to entropy density ratio ($\eta/s$) [1] . To study the flow,
we use Fourier series expansion of the azimuthal distribution of the invariant
yield given by

$E\frac{d^{3}N}{dp^{3}}=\frac{1}{2\pi}\frac{d^{2}N}{p_{T}dp_{T}dy}\left[1+\sum_{n}2v_{n}\cos
n(\phi-\Psi_{R})\right].$ (2)

Here, $v_{1}$, $v_{2}$ and $v_{3}$ represent directed flow, elliptic flow and
triangular flow respectively. The $n^{th}$-order flow coefficient is given by

$v_{n}=\langle\cos n(\phi-\psi_{n})\rangle,$ (3)

where, $\psi_{n}$ is the $n^{th}$ order event plane. A detailed description of
event plane and its calculation can be found in Ref. [2].

We have used strange and multi-strange hadrons in our study since they are
expected to freeze out earlier and have small hadronic interaction cross
section [3] compared to other charged hadrons. So they get least affected by
the late-stage hadronic interaction phase of the system.

##  2 Analysis details

Since the particles used in our analysis are short-lived and decay before
reaching the detector, we have reconstructed them by using the invariant mass
of their daughter particles. The decay channels and branching ratios of these
particles are given below.

$K_{S}^{0}\rightarrow\pi^{+}+\pi^{-}$ (B.R = 69.20%)

$\Lambda(\bar{\Lambda})\rightarrow p(\bar{p})+\pi^{-}(\pi^{+})$ (B.R = 63.9%)

$\Xi(\bar{\Xi^{+}})\rightarrow\Lambda(\bar{\Lambda})+\pi^{-}(\pi^{+})$ (B.R =
99.887%)

$\Omega(\bar{\Omega^{+}})\rightarrow\Lambda(\bar{\Lambda})+K^{-}(K^{+})$ (B.R
= 67.8%)

$\phi\rightarrow K^{+}+K^{-}$ (B.R = 49.1%)

For particles like $K_{S}^{0}$, $\Lambda$, $\Xi$ and $\Omega$ weak decay
topology cuts are applied to reduce the combinatorial background. For
$\Lambda$ hyperons weak decay feed-down correction is implemented. In the
upper panel of Fig. 1, the invariant mass distribution is shown for $K^{+}$
and $K^{-}$ to reconstruct the $\phi$ mesons as an example. The combinatorial
background is calculated by event-mixing method and subtracted from the
signal+background distribution to get the signal of $\phi$ mesons.

Figure 1: The upper panel shows the signal+background distribution and signal
peak of phi meson in Au+Au collision at $\sqrt{s_{NN}}$ = 54.4 GeV in $p_{T}$
bin 0.4-0.6 GeV/$c$. The lower panel shows $v_{2}$ as a function of invariant
mass of $K^{+}$ and $K^{-}$. The distribution is fitted with a function given
in Eq. 4.

For flow coefficient measurements we have used invariant mass method [4]. As
an example we have explain $v_{2}$ calculation of $\phi$ mesons using this
method. In this method we need to calculate $v_{2}$ as a function of invariant
mass of decay daughters, $K^{+}$ and $K^{-}$ as shown in the lower panel of
Fig. 1. The distribution is fitted with a function given in Eq. 4, where
$v_{2}^{S+B}$ is the combined $v_{2}$ of signal+background. $v_{2}^{S}$ and
$v_{2}^{B}$ represent $v_{2}$ of signal and background respectively.
$v_{2}^{B}$ is parametrized as a first order polynomial of invariant mass.
$v_{2}^{S}$ is a free parameter and can be extracted from the fitting. Same
method was used to calculate $v_{2}$ and $v_{3}$ of other particles as well.

$v_{2}^{S+B}(M_{inv})=v_{2}^{S}\frac{S}{S+B}(M_{inv})+v_{2}^{B}\frac{B}{S+B}(M_{inv}).$
(4)

##  3 Results and Discussion

Figure 2 shows the nuclear modification factor ($R_{CP}$) at $\sqrt{s_{NN}}$ =
54.4 GeV, which is calculated by taking ratio of particle yields in 0-5% and
40-60% for $K_{S}^{0}$, $\Lambda+\bar{\Lambda}$, $\Xi^{-}+\bar{\Xi^{+}}$ and
for $\phi$ mesons. The ratio shows a strong suppression at high $p_{T}$
indicating partonic energy loss inside the medium.

Figure 2: The nuclear modification factor, $R_{CP}$ is plotted as a function
of $p_{T}$ at $\sqrt{s_{NN}}$ = 54.4 GeV for $K_{S}^{0}$,
$\Lambda+\bar{\Lambda}$, $\Xi^{-}+\bar{\Xi^{+}}$ and for $\phi$. The vertical
lines represent the statistical error bars and the rectangular boxes represent
the systematic error bars. Systematic errors are not shown for $\phi$ mesons.
The blue box represents the uncertainty from the determination of the number
of binary collisions.

The ratio, $N[\Omega+\bar{\Omega}^{+}]/2N[\phi]$ is calculated for different
centralities as shown in Fig. 3. The ratio shows an enhancement at
intermediate $p_{T}$ and the enhancement is larger for central collisions than
peripheral collisions which suggests the idea of quark coalescence
hadronization [5]. The ratio at 54.4 GeV is compared with other STAR energies
[6] for most central collisions as shown in Fig. 4.
$N[\Omega+\bar{\Omega}^{+}]/2N[\phi]$ ratio shows the same trend and is
consistent with other energies except for 11.5 GeV data where the larger
statistical error bar restrict us from making any conclusion.

Figure 3: $N[\Omega+\bar{\Omega}^{+}]/2N[\phi]$ is plotted as a function of
$p_{T}$ for different centrality classes at $\sqrt{s_{NN}}$ = 54.4 GeV. Only
statistical error bars are shown. Figure 4:
$N[\Omega+\bar{\Omega}^{+}]/2N[\phi]$ ratio vs. $p_{T}$ for 0-10% centrality
at $\sqrt{s_{NN}}$ = 54.4 GeV is compared with other STAR energies for most
central collisions. For 54.4 GeV only statistical errors are shown. Figure 5:
Left panel shows the $v_{2}$ as a function of $p_{T}$ for 0-80% centrality at
$\sqrt{s_{NN}}$ = 54.4 GeV for $K^{0}_{S}$, $\Lambda$, $\phi$, $\Xi^{-}$ and
$\Omega^{-}$. The right panel shows the same for $v_{3}$. The vertical bars
here are the statistical error bars and the shaded boxes are the systematic
error bars. Figure 6: The left panels shows the NCQ scaled $v_{2}$ as a
function of transverse kinetic energy at $\sqrt{s_{NN}}$ = 27 GeV for 0-80%
centrality. The right panel shows the same for $\sqrt{s_{NN}}$ = 54.4 GeV.
Only statistical error bars are shown.

The elliptic ($v_{2}$) and triangular ($v_{3}$) flow coefficients are measured
for the strange and multi-strange hadrons for minimum bias events at
$\sqrt{s_{NN}}$ = 54.4 GeV as shown in Fig. 5. A mass ordering is observed at
low $p_{T}$ ($<$ 2 GeV/$c$), where the particles with smaller mass show a
higher flow which is due to the radial flow of the system during the QGP
phase. A particle type dependence is observed in both $v_{2}$ and $v_{3}$ at
high $p_{T}$ region. We have studied the NCQ scaled $v_{2}$ as a function of
NCQ scaled transverse kinetic energy at $\sqrt{s_{NN}}$ = 27 and 54.4 GeV for
0-80% centrality as shown in Fig. 6. The NCQ scaled $v_{2}$ falls on a single
curve for all the hadron, which indicates that the collectivity is developed
during the initial QGP phase of the system. Similar behaviour is observed in
case of NCQ scaled $v_{3}$ at $\sqrt{s_{NN}}$ = 54.4 GeV as shown in Fig. 7.

Figure 7: NCQ scaled $v_{3}$ is plotted as a function of transverse kinetic
energy at $\sqrt{s_{NN}}$ = 54.4 GeV for 0-80% centrality. The vertical bars
here are the statistical error bars and the shaded boxes are the systematic
error bars.

##  4 Summary

In summary, we have measured the $p_{T}$ spectra of strange hadrons,
$K^{0}_{S}$, $\Lambda$, $\bar{\Lambda}$, $\phi$, $\Xi^{-}$, $\bar{\Xi}^{+}$,
$\Omega^{-}$, $\bar{\Omega}^{+}$ at midrapidity in Au+Au collisions at
$\sqrt{s_{NN}}$ = 54.4 GeV. The nuclear modification factor, $R_{CP}$ for the
above mentioned particles has been measured which shows a strong suppression
at high $p_{T}$. The baryon to meson ratio,
$N[\Omega+\bar{\Omega}^{+}]/N[\phi]$ shows an enhancement at intermediate
$p_{T}$ which is higher for central collisions than that for peripheral
collisions. The anisotropic flow coefficients, $v_{2}$ and $v_{3}$, have been
measured as a function of $p_{T}$ in 0-80% centrality for all the strange
particles mentioned above. A mass ordering is observed at low $p_{T}$ and a
baryon to meson separation is observed at high $p_{T}$. The NCQ scaling founds
to hold for both $v_{2}$ and $v_{3}$.

## References

  * [1] Chun Shen et al., J. Phys. G: Nucl. Part. Phys. 38, 124045 (2011)
  * [2] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58, 1671 (1998)
  * [3] A. Shor, Phys. Rev. Lett. 54, 1122 (1985)
  * [4] N. Borghini and J.-Y. Ollitrault, Phys. Rev. C 70, 064905 (2004)
  * [5] J. H. Chen et al, Phys. Rev. C 78, 034907 (2008)
  * [6] J. Adam et al. (STAR ), Phys. Rev. C 102, 034909 (2020)