#  Calculation of Thermodynamic Characteristics and
Sound Velocity for Two-Dimensional Yukawa Fluids Based on a
Two-Step Approximation for the Radial Distribution Function

Ilnaz I. Fairushin  Department of Computational Physics, Institute of Physics,
Kazan Federal University, 420008 Kazan, Russia  Anatolii V. Mokshin
Department of Computational Physics, Institute of Physics, Kazan Federal
University, 420008 Kazan, Russia  Department of Computational Physics,
Institute of Physics, Kazan Federal University, 420008 Kazan, Russia
Department of Computational Physics, Institute of Physics, Kazan Federal
University, 420008 Kazan, Russia

###### Abstract

We propose a simple two-step approximation for the radial distribution
function of a one-component two-dimensional Yukawa fluid. This approximation
is specified by the key parameters of the system: coupling parameter and
screening parameter. On the basis of this approximation, analytical
expressions are obtained for the same thermodynamic quantities as internal
energy, internal pressure, excess entropy in the two-particle approximation,
and also longitudinal sound velocity. The theoretical results show an
agreement with the results obtained in the case of a true radial distribution
function.

The model system of charged particles with a screened Coulomb (Yukawa)
interaction is widely used for modeling and predicting the physical properties
of a wide variety of real fluids such as simple neutral fluids, liquid metals,
colloidal solutions, microemulsions, etc. [1, 2]. This system represents
similarly charged particles of the same mass (ions, dusty, or colloidal
particles) surrounded by a background of the opposite sign, which in practice
usually consists of particles of a smaller mass (electrons or ions in the case
of dusty plasma or colloidal solutions) [2]. The interaction potential energy
$u$ of a pair of charged particles on the distance $r$ in this case is written
as

$u(r)=\frac{(Ze)^{2}\exp(-r/\lambda_{s})}{4\pi\varepsilon_{0}r}\,,$ (1)

where $Z$ is the charge number, $e$ is the electron charge, $\varepsilon_{0}$
is the vacuum permittivity, and $\lambda_{s}$ is the Debye screening length
associated with the presence of a neutralizing medium, which depends on the
concentration and temperature of the background particles and determines how
the interaction of the main particles will differ from the simple Coulomb
interaction.

In the last decade, researchers have increased their attention on two-
dimensional Yukawa systems [1, 2, 3, 4, 5, 6, 7, 8, 9]. Firstly, this is due
to the fact that such systems can be easily implemented in experiments. These
include experiments with a monolayer of charged macroparticles in a low-
temperature plasma [5, 6] and experiments with colloidal solutions with
microsized particles [9]. Within the framework of these experiments, it is
possible to directly monitor the dynamics of individual particles, as well as
to study the features of phase transitions in two-dimensional systems [5].
Second, the Yukawa interaction potential is expressed by a simple analytical
expression. Thus, a system of particles interacting through a potential of the
form (1) can be considered as a convenient model of a one-component substance,
on the example of which it is convenient to test one or another microscopic
theory describing the structure or dynamics of simple substances. This is
especially important for liquid-like disordered systems, which, unlike gases
and crystalline bodies, are characterized by the absence of a suitable small
parameter for the development of an appropriate theoretical description [10,
11, 12, 13]. The key characteristic of the liquid structure and its short-
range order is the radial distribution function (RDF) $g(r)$. Knowing the RDF
for systems with a known interparticle interaction potential $u(r)$, one can
directly calculate such thermodynamic parameters of many-particle systems as
internal energy, internal pressure, and excess entropy in the pairwise
approximation [14, 15, 16, 17, 18, 19, 20]. Note that the RDF is a special
case of the radial basis functions, which are widely used in various problems
of fluid mechanics [21, 22].

The specificity of the interparticle interaction in the case of the Yukawa
system is determined by two key dimensionless parameters: coupling parameter
$\Gamma$ and screening parameter $\kappa$ [1, 2, 23]. The coupling parameter

$\Gamma=\frac{(Ze)^{2}}{4\pi\varepsilon_{0}ak_{B}T}\ $ (2)

is the ratio of the average potential energy of interaction (without
screening) to the average energy of the particle thermal motion. In the
expression (2), the quantity $a=(\pi\rho)^{-1/2}$ is half the average
interparticle distance or the so-called radius of the Wigner–Seitz cell,
$\rho$ is the number particles per unit area of a two-dimensional system,
$k_{B}$ is the Boltzmann constant, and $T$ is the absolute temperature of the
system. The screening parameter $\kappa$ is defined as the ratio of $a$ to the
Debye length $\lambda_{s}$:

$\kappa=\frac{a}{\lambda_{s}}\,.$ (3)

The time scale of the charge density fluctuations in the system is determined
by the plasma frequency

$\omega_{p}=\sqrt{\frac{Z^{2}e^{2}\rho}{2a\varepsilon_{0}m}}\,,$ (4)

where $m$ is the particle mass.

In [24], on the bases of the quasi-localized charge approximation (QLCA) [25,
26, 27], simple analytical expressions are obtained that describe the
dispersions of longitudinal and transverse collective excitations in a three-
dimensional Yukawa fluid. This goal was achieved by the authors by using the
so-called one-step (1st) approximation for the RDF:

$g^{(1\rm st)}(x)=\theta(x-x_{\rm eff})\\\ ,$ (5)

where $\theta(x)$ is the Heaviside function, $x_{\rm eff}$ is the effective
radius of the particle, which is an input parameter in this approximation, and
which is determined from the known data on the internal pressure and internal
energy systems. It should be noted that the one-step approximation (5)
actually corresponds to the function $g(x)$ of a highly rarefied gas of
absolutely hard spheres of radius $x_{\rm eff}$. This approximation ignores
the presence of local short-range order in liquids, which manifests itself in
a characteristic maximum in the function $g(x)$.

In this paper, we propose a two-step approximation for the function $g(x)$ of
a two-dimensional Yukawa system. It should be noted that, initially, this
approximation was proposed for the three-dimensional system [28]. Here, two
key parameters of the Yukawa system are used as input parameters: $\Gamma$ and
$\kappa$. Within the framework of this approximation, the internal energy,
internal pressure, excess entropy in the two-particle approximation, and the
longitudinal sound velocity of a two-dimensional Yukawa fluid are calculated.
The ($\Gamma,\;\kappa$)-states of the Yukawa fluid will be considered, where
$\Gamma$ = 20; 50; and 100 and $\kappa$ = 1; 1.5; and 2. The theoretical
results are compared with the calculation results based on the true RDF
obtained by us using molecular dynamics (MD) simulations. The equilibrium
molecular dynamics simulations of the Yukawa fluid for $\Gamma=20,50,$ and
$100$ and $\kappa=1,1.5,$ and $2$ were carried out using the computational
package LAMMPS [29]. The simulation was performed for a system consisting of
2500 particles interacting through a potential (1) in a square cell, on which
periodic boundary conditions were imposed. The calculations were performed in
the $NVT$ ensemble. The particle motion equations were integrated in
accordance with the Verlet algorithm with a time integration step $t_{\rm
step}=0.01/\omega_{p}$. Averaging over 10,000 time steps was used to calculate
the RDF.

The two-step (2st) approximation for the function $g(x)$ is provided as:

$g^{(2\rm st)}(x)=g_{m}\theta(x-x_{1})\theta(x_{2}-x)+\theta(x-x_{2}).$ (6)

Here, the distances $x_{1}$ and $x_{2}$ determine the position and width of
the first maximum of the function $g(x)$ on the $x$ axis, and $g_{m}$ is the
height of this maximum. The presence of the maxima corresponding to the
second, third, and other co-ordinations are ignored within the framework of
the two-step approximation (6). Following [30], we define the distance $x_{1}$
as the size of the region of absence of interparticle correlations, which is
provided by the condition $g(x_{1})=0.5$. Further, the distance $x_{1}$ can be
related to the coupling parameter $\Gamma$ and the screening parameter
$\kappa$ [30]:

$x^{2}_{1}=\frac{1}{b_{1}}\ln\frac{\Gamma-b_{2}(\kappa)}{b_{3}(\kappa)},$ (7)

where

$b_{1}=2.434,$ (8) $b_{2}(\kappa)=-5.21+6.866\kappa-2.492\kappa^{2}\,,$ (9)
$b_{3}(\kappa)=0.712-0.572\kappa+0.437\kappa^{2}\,.$ (10)

The values of the constant coefficients in these polynomials can be determined
by numerically solving a system of seven non-linear Equations (7) written for
seven different states with $\Gamma=20$, $50$, and $100$; $\kappa=1$, and $2$,
and also $\Gamma=50$, $\kappa=1.5$. In this case, the distances $x_{1}$ for
these states were determined from the true $g(r)$ obtained from the results of
our MD simulation. The solution of the system of equations was carried out
using the modified Newton method, the accuracy of the numerical solution of
the system of equations was at least $99.7$ % (this is quite sufficient for
the purposes of this work). On the other hand, the distance $x_{2}$, which
determines the size of the first coordination shell within the framework of
the (6) approximation, can be found from the condition of the charge
neutrality of the system under consideration, which in the two-dimensional
case is written as

$\int_{0}^{\infty}[1-g(x)]xdx=\frac{1}{2}.$ (11)

From the expression (11) and taking into account the relation (6), we obtain:

$x^{2}_{2}=\frac{x^{2}_{1}g_{m}-1}{g_{m}-1}.$ (12)

The value of $g_{m}$ can be determined from the relation found in [30],

$\Gamma=a_{1}(\kappa)+a_{2}(\kappa)g_{m}+a_{3}(\kappa)g_{m}^{2},$ (13)

where the $\kappa$-dependence of the parameters $a_{1}$, $a_{2}$, and $a_{3}$
are providing by using a second degree polynomial

$a_{\xi}(\kappa)=c_{1}^{(\xi)}+c_{2}^{(\xi)}\kappa+c_{3}^{(\xi)}\kappa^{2},\,\,\xi=1,\,2,\,3\,.$
(14)

For the case of a two-dimensional Yukawa fluid, the values of the
dimensionless parameters $c_{1}^{(\xi)}$, $c_{2}^{(\xi)}$, and $c_{3}^{(\xi)}$
were are found by solving a system of nine cases of Equation (13) written for
states with $\Gamma=20$, $50$, and $100$, and $\kappa=1$, $1.5$, and $2$. The
solution was also carried out using the modified Newton method; the error of
the numerical solution of this system of equations in this case did not exceed
$0.03$ %. The values for $g_{m}$ were collected from MD simulation data. As a
result of solving this system, it was obtained that: $c_{1}^{(1)}=-248.56$,
$c_{2}^{(1)}=369.596$, and $c_{3}^{(1)}=-126.792$; $c_{1}^{(2)}=272.541$,
$c_{2}^{(2)}=-427.101$, and $c_{3}^{(2)}=137.019$; and $c_{1}^{(3)}=-58.808$,
$c_{2}^{(3)}=112.002$, and $c_{3}^{(3)}=-28.326$. The relation (13) is
actually the equation of state for the equilibrium liquid phase of the two-
dimensional Yukawa system, from which we find

$g_{m}=\frac{-a_{2}(\kappa)+\sqrt{a_{2}^{2}(\kappa)-4a_{1}(\kappa)a_{3}(\kappa)+4a_{3}(\kappa)\Gamma}}{2a_{3}(\kappa)}.$
(15)

Thus, the value of $x_{2}$ for a particular $(\Gamma,\;\kappa)$-state is
determined by solving the system of Equations (7), (12), and (15).

The RDF is included in microscopic expressions for many physical
characteristics. Thus, the excess internal energy $U_{\rm ex}$ of a two-
dimensional Yukawa fluid (in units of $k_{B}T$) is [14, 15, 16, 17, 18]

$U_{\rm ex}=\Gamma\int_{0}^{\infty}\ \exp(-\kappa x)g(x)\,dx.$ (16)

Taking into account the approximation (6), from the expression (16) we obtain

$U_{\rm ex}^{(2\rm st)}=\frac{\Gamma g_{m}}{\kappa}\biggl{[}\exp(-\kappa
x_{1})-\frac{g_{m}-1}{g_{m}}\exp(-\kappa x_{2})\biggr{]}.$ (17)

Further, for the excess internal pressure $P_{\rm ex}$ of a two-dimensional
Yukawa fluid (in units of $\rho k_{B}T$), we have: [14, 15, 16, 17, 18]

$P_{\rm ex}=\frac{\Gamma}{2}\int_{0}^{\infty}\ \bigl{(}\kappa
x+1\bigr{)}\exp(-\kappa x)g(x)\,dx\,.$ (18)

Hence, taking into account the approximation (6), we find

$P_{\rm ex}^{(2\rm st)}=\frac{\Gamma g_{m}}{2\kappa}\biggl{[}\frac{\kappa
x_{1}+2}{\exp(-\kappa x_{1})}-\frac{g_{m}-1}{g_{m}}\frac{\kappa
x_{2}+2}{\exp(\kappa x_{2})}\biggr{]}.$ (19)

The microscopic expression for the excess entropy $S_{ex2}$ in units of $\rho
k_{B}$ in the two-particle approximation does not explicitly contain the
interaction potential $u(r)$, and for a two-dimensional isotropic system has
the form [19]:

$S_{\rm ex2}=-\int_{0}^{\infty}[g(x)\ln g(x)+1-g(x)]x\,dx.$ (20)

Then, within the approximation (6), from the expression (20), we obtain

$S_{\rm ex2}^{(2\rm st)}=-\frac{1}{2(g_{m}-1)}\left[x_{1}^{2}g_{m}\ln
g_{m}-\left(g_{m}\ln g_{m}+1-g_{m}\right)\right].$ (21)

Knowing the RDF using the QLCA model, for a two-dimensional Yukawa fluid, one
can calculate the longitudinal sound velocity $c_{L}$ (in units of thermal
velocity$v_{\rm{th}}=\sqrt{k_{B}T/m}$) [17]:

$c_{L}^{2}=\frac{\Gamma}{8}\int_{0}^{\infty}\bigl{(}3(\kappa x)^{2}+5\kappa
x+5\bigr{)}\exp(-\kappa x)g(x)\,dx.$ (22)

Within the framework of the approximation (6), from the expression (22), we
obtain

$\begin{gathered}c_{L}^{2(2\rm st)}=\frac{\Gamma
g_{m}}{8\kappa}\biggl{[}\frac{3(\kappa x_{1})^{2}+11\kappa
x_{1}+16}{\exp(\kappa x_{1})}\\\
-\left(\frac{g_{m}-1}{g_{m}}\right)\frac{3(\kappa x_{2})^{2}+11\kappa
x_{2}+16}{\exp(\kappa x_{2})}\biggr{]}.\end{gathered}$ (23)

As can be seen from the expressions (17), (19), (21), and (23), the values
$U_{ex}$, $P_{ex}$, $S_{ex2}$, and $c_{L}$ can be directly calculated for a
given ($\Gamma,\;\kappa$)-state. On the other hand, if the true function
$g(x)$ is known, then the quantities can be estimated using microscopic
expressions (16), (18), (20), and (22).

The results of the numerical calculations of the reduced excess internal
energy $U_{ex}$, the reduced excess internal pressure $P_{ex}$, and the
reduced excess entropy $S_{ex2}$ performed within the (6) approximation for
$g(x)$, as well as using the true $g(x)$ from the MD simulation, are presented
in Table 1. This table also shows the relative correspondences between the
theoretical results and simulation results. For most
($\Gamma,\;\kappa$)-states, the differences between the theoretical results
and simulation data for $U_{ex}$ and $P_{ex}$ do not exceed $2$ %. The largest
discrepancies corresponding to $3.818$ % for $U_{ex}$ and $2.912$ % for
$P_{ex}$ are observed for states with maximum $\kappa=2$, which can be
explained by violation of the charge neutrality condition (11) in the case of
states with this value of the screening parameter $\kappa$. Further, the
entropy $S_{ex2}$ is very structure sensitive. This may explain the weak
agreement between the results of theoretical calculations and the data of the
MD simulation.

Table 1: Reduced excess internal energy $U_{\rm ex}$, reduced excess internal pressure $P_{\rm ex}$, and reduced excess entropy $S_{\rm ex2}$ of a 2D Yukawa fluid found using $g(x)$ from MD simulations. The same quantities ($U_{\rm ex}^{(2\rm st)}$, $P_{\rm ex}^{(2\rm st)}$, and $S_{\rm ex2}^{(2\rm st)}$) are calculated based on the expressions (6), (19), and (21). The relative correspondences $\delta_{U_{\rm ex}}$ , $\delta_{P_{\rm ex}}$, and $\delta_{S_{\rm ex2}}$ of these quantities in % between theoretical results and simulation results. $\kappa$ | $\Gamma$ |  $U_{\rm ex}$ | $U_{\rm ex}^{(2\rm st)}$ | $\delta_{U_{\rm ex}}$ | $P_{\rm ex}$ | $P_{\rm ex}^{(2\rm st)}$ | $\delta_{P_{\rm ex}}$ | $S_{\rm ex2}$ | $S_{\rm ex2}^{(2\rm st)}$ | $\delta_{S_{\rm ex2}}$
---|---|---|---|---|---|---|---|---|---|---
1 | 20 |  6.762 | 6.767 | 0.072 | 10.554 | 10.536 | 0.165 |  $-$0.643 |  $-$0.782 | 21.576
1 | 50 |  15.901 | 15.951 | 0.319 | 25.518 | 25.461 | 0.225 |  $-$1.111 |  $-$1.068 | 3.883
1 | 100 |  30.943 | 30.641 | 0.977 | 50.294 | 49.719 | 1.142 |  $-$1.803 |  $-$1.348 | 25.226
1.5 | 20 |  2.732 | 2.708 | 0.871 | 4.958 | 4.919 | 0.797 |  $-$0.542 |  $-$0.679 | 25.377
1.5 | 50 |  5.963 | 6.068 | 1.751 | 11.403 | 11.470 | 0.591 |  $-$0.897 |  $-$0.941 | 4.927
1.5 | 100 |  11.165 | 11.266 | 0.910 | 21.912 | 21.853 | 0.267 |  $-$1.402 |  $-$1.190 | 15.117
2 | 20 |  1.336 | 1.285 | 3.818 | 2.697 | 2.619 | 2.912 |  $-$0.454 |  $-$0.581 | 27.931
2 | 50 |  2.647 | 2.699 | 1.982 | 5.766 | 5.823 | 0.986 |  $-$0.716 |  $-$0.808 | 12.893
2 | 100 |  4.665 | 4.815 | 3.225 | 10.601 | 10.756 | 1.468 |  $-$1.053 |  $-$1.026 | 2.562

Figure 1 shows the dependence of the reduced excess internal energy $U_{ex}$
and reduced excess internal pressure $P_{ex}$ on the coupling parameter
$\Gamma$ for the values of the screening parameter $\kappa$ outside the range
that was used to construct a two-step approximation for RDF (6). Here, we also
show the data from the work [20], in which the energy and pressure values of
the two-dimensional Yukawa fluid were calculated in a wide range of changes in
the value of the $\Gamma$ and $\kappa$. It can be seen that, even at these
$\kappa$ values, the analytical formulas obtained in this work for the direct
calculation of reduced excess internal energy and pressure of two-dimensional
Yukawa fluids are in good agreement with the simulation data.

Figure 1: Dependence of the reduced excess internal energy $U_{ex}$ and the
reduced excess internal pressure $P_{ex}$ on the coupled parameter $\Gamma$
for the values of the screening parameter $\kappa=0.5$ (a) and $\kappa=2.4$
(b), which are constructed using expressions (17) and (19). Symbols show data
from work [20].

Table 2 shows the results of the numerical calculations of the reduced
longitudinal sound velocity $c_{L}$, performed within the approximation (6)
for $g(x)$, as well as using the true $g(x)$ from MD modeling. This table also
shows the relative correspondences between theoretical results and simulation
results. The greatest discrepancy is observed for the state of the two-
dimensional Yukawa liquid with $\Gamma=20$ and $\kappa=2$, which is the
closest of all those considered to the gas phase. This feature is related to
the fact that the approximation of a quasi-localized charge better describes
the collective properties of a Yukawa liquid with states close to crystalline,
i.e., with large $\Gamma$ and small $\kappa$ [25].

Table 2: Longitudinal sound velocity in units of thermal velocity $c_{L}$ of a 2D Yukawa fluid, found using $g(x)$ from MD simulations. The same value $c_{L}^{(2\rm st)}$ calculated based on the expression (23). The relative correspondences $\delta_{c_{L}}$ of this value in % between theoretical and simulation results are also provided. $\kappa$ | $\Gamma$ |  $c_{L}$ | $c_{L}^{(2\rm st)}$ | $\delta_{c_{L}}$
---|---|---|---|---
1 | 20 |  5.196 | 5.189 | 0.139
1 | 50 |  8.145 | 8.13 | 0.184
1 | 100 |  11.475 | 11.422 | 0.465
1.5 | 20 |  3.777 | 3.761 | 0.434
1.5 | 50 |  5.831 | 5.827 | 0.077
1.5 | 100 |  8.156 | 8.122 | 0.422
2 | 20 |  2.917 | 2.881 | 1.266
2 | 50 |  4.391 | 4.392 | 0.030
2 | 100 |  6.051 | 6.055 | 0.063

In order to check the validity of the relation (23) in the case of values of
the coupling and screening parameters outside the range used in this work when
constructing the two-step approximation for the $g(x)$, we calculated the
longitudinal sound velocity for the gamma–kappa states collected from [17].
The results are presented in Table 3. It can be seen that, even at
$\Gamma=1033$ and $\kappa=3$, formula (23) provides a deviation of less than
4% from the value of s calculated using the true $g(x)$.

Table 3: The same as in Table 2, except for the third column; here are the data from work [17]. $\kappa$ | $\Gamma$ |  $c_{L}$ | $c_{L}^{(2\rm st)}$ | $\delta_{c_{L}}$
---|---|---|---|---
1 | 163 |  14.62 | 14.51 | 0.75
2 | 362 |  11.25 | 11.02 | 2.04
3 | 1033 |  10.23 | 9.87 | 3.52

The results of this work indicate the following. The thermodynamic
characteristics and the sound velocity for a two-dimensional Yukawa fluid can
be theoretically calculated from microscopic expressions, where the
characteristic parameters of the Yukawa system (coupling parameter and
screening parameter) are used as input parameters. The two-step approximation
proposed in this work for the RDF provides a good agreement with the
simulation results for such quantities as the internal energy, internal
pressure, and longitudinal sound velocity. If we compare the results for the
two-dimensional Yukawa fluid with the three-dimensional case [28], we can see
that the obtained analytical expressions (17) and (19) are simpler and that,
at the same time, they provide the same accuracy when reproducing the
simulations results. In addition, as in the case of the three-dimensional
Yukawa fluid in the two-dimensional case, to correctly calculate the excess
entropy, it is necessary to use a more accurate model for $g(r)$.

###### Acknowledgements.

This work was supported by the Russian Science Foundation (Project No.
19-12-00022). The authors are grateful to S.A. Khrapak for helpful
discussions.

## References

  * [1] Fortov, V.E.; Ivlev, A.V.; Khrapak, S.A.; Khrapak, A.G.; Morfill, G.E. Phys. Rep. 2005, 421, 1–103.
  * [2] Dubin, D.H.E.; O’Neill, T.M. Rev. Mod. Phys. 1999, 71, 87–172.
  * [3] Klumov, B.A.; Khrapak, S.A. Results Phys. 2020, 17, 103020.
  * [4] Klumov, B.A. JETP Lett. 2022, 115, 108–113.
  * [5] Vasilieva, E.V.; Petrov, O.F.; Vasiliev, M.M. Sci. Rep. 2021, 11, 523.
  * [6] Kononov, E.A.; Vasiliev, M.M.; Petrov, O.F.; Vasilieva, E.V. Nanomaterials 2021, 11, 2931.
  * [7] Fairushin, I.I.; Vasiliev, M.M.; Petrov, O.F. Molecules 2021, 26, 6974.
  * [8] Fairushin, I.I.; Petrov, O.F.; Vasiliev, M.M. J. Exp. Theor. Phys. 2020, 130, 477–481.
  * [9] Lin, B.J.; Chen, L.J. J. Chem. Phys. 2007, 126, 34706.
  * [10] Hansen, J.P.; McDonald, I.R. Theory of Simple Liquids; Academic Press: London, UK, 2006.
  * [11] Mokshin, A.V. Theor. Math. Phys. 2015, 183, 449–477.
  * [12] Khrapak, S. Appl. Sci. 2022, 12, 7939.
  * [13] Tareyeva, E.E.; Ryzhov, V.N. Theor. Math. Phys. 2016, 189, 1806–1817.
  * [14] Farouki, R.T.; Hamaguchi, S. J. Chem. Phys. 1994, 101, 9885–9893.
  * [15] Hartmann, P.; Kalman, G.J.; Donkó, Z.; Kutasi, K. Phys. Rev. E 2005, 72, 026409.
  * [16] Khrapak, S.A.; Kryuchkov, N.P.; Yurchenko, S.O.; Thomas, H.M. J. Chem. Phys. 2015 ,142, 194903.
  * [17] Khrapak, S.A. Phys. Plasm. 2016, 23, 024504.
  * [18] Filippov, A.V.; Reshetnyak, V.V.; Starostin, A.N.; Tkachenko, I.M.; Fortov, V.E. JETP Lett. 2019, 110, 659–666.
  * [19] Baranyai, A.; Evans, D.J. Phys. Rev. A 1989, 40, 3817.
  * [20] Kryuchkov, N.P.; Khrapak, S.A.; Yurchenko, S.O. J. Chem. Phys. 2017, 146, 134702.
  * [21] Geronzi, L.; Gasparotti, E.; Capellini, K.; Cella, U.; Groth, C.; Porziani, S.; Chiappa, A.; Celi, S.; Biancolini, M.E. J. Comput. Sci. 2021, 51, 101327.
  * [22] Groth, C.; Porziani, S.; Biancolini, M.E. Fluids 2021, 6, 314.
  * [23] Mokshin, A.V.; Fairushin, I.I.; Tkachenko, I.M. Phys. Rev. E 2022, 105, 025204.
  * [24] Khrapak, S.; Klumov, B.; Couedel, L.; Thomas, H. Phys. Plasm. 2016, 23, 023702.
  * [25] Kalman, G.J.; Golden, K.I. Phys. Rev. A 1990, 41, 5516–5527.
  * [26] Golden, K.I.; Kalman, G.J. Phys. Plasm. 2000, 7, 14–32.
  * [27] Donkó, Z.; Kalman, G.J.; Hartmann, P. J. Phys. Condens. Matter 2008, 20, 413101.
  * [28] Fairushin, I.I.; Khrapak, S.A.; Mokshin, A.V. Results Phys. 2020, 19, 103359.
  * [29] Plimpton, S. J. Comput. Phys. 1995, 117, 1.
  * [30] Ott, T.; Bonitz, M.; Stanton, L.; Murillo, M.S. Phys. Plasm. 2014, 21, 113704.