11institutetext: RWTH Aachen University, Aachen, Germany
11email<EMAIL_ADDRESS>
# Exact Probabilistic Inference Using Generating Functions††thanks: Extended
abstract accepted by LAFI 2023 – the Languages for Inference workshop co-
located with POPL 2023.

Lutz Klinkenberg and Tobias Winkler and Mingshuai Chen and Joost-Pieter Katoen

Probabilistic programs are typically normal-looking programs describing
posterior probability distributions. They intrinsically code up randomized
algorithms and have long been at the heart of modern machine learning and
approximate computing. We explore the theory of _generating functions_ [19]
and investigate its usage in the exact quantitative reasoning of probabilistic
programs. Important topics include the exact representation of program
semantics [13], proving exact program equivalence [5], and – as our main focus
in this extended abstract – exact probabilistic inference.

In probabilistic programming, _inference_ aims to derive a program’s posterior
distribution. In contrast to approximate inference, inferring _exact_
distributions comes with several benefits [8], e.g., no loss of precision,
natural support for symbolic parameters, and efficiency on models with certain
structures. Exact probabilistic inference, however, is a notoriously hard task
[6, 12, 17, 18]. The challenges mainly arise from three program constructs:
(1) unbounded while-loops and/or recursion,  (2) infinite-support
distributions, and  (3) conditioning (via posterior observations).  We present
our ongoing research in addressing these challenges (with a focus on
conditioning) leveraging generating functions and show their potential in
facilitating exact probabilistic inference for discrete probabilistic
programs.

##  1 Inference in Probabilistic Programs

### State-of-the-Art.

Most existing probabilistic programming languages implement _sampling_ -based
inference algorithms rooted in the principles of Monte Carlo [15], thereby
yielding numerical approximations of the exact results, see, e.g., [9]. In
terms of semantics, many probabilistic systems employ _probability density
function_ (PDF) representations of distributions, e.g., ($\lambda$)PSI [7, 8],
AQUA [11], Hakaru [16], and the density compiler in [3, 4]. These systems are
dedicated to inference (with conditioning) for programs encoding (joint
discrete-)continuous distributions. Reasoning about the underlying PDF
representations, however, amounts to resolving complex integral expressions in
order to answer inference queries, thus confining these techniques either to
(semi-)numerical methods [3, 4, 11, 16] or exact methods yet limited to
bounded looping behaviors [7, 8]. Dice [10] employs weighted model counting to
enable exact inference for discrete probabilistic programs, yet is also
confined to statically bounded loops. The tool Mora [1, 2] supports exact
inference for various types of Bayesian networks, but relies on a restricted
form of intermediate representation known as prob-solvable loops.

### The PGF Approach.

Klinkenberg et al. [13] provide a program semantics that allows for exact
quantitative reasoning about probabilistic programs without conditioning. They
exploit a denotational approach à la Kozen [14] and treat a probabilistic
program as a _distribution transformer_ , i.e., mapping a distribution over
the inputs (the prior) into a distribution after execution of the program (the
posterior). In [13], the domain of discrete distributions is represented in
terms of _probability generating functions_ (PGFs), which are a special kind
of generating functions [19]. This representation comes with several benefits:
(a) it naturally encodes common, infinite-support distributions (and
variations thereof) like the geometric or Poisson distribution in compact,
_closed-form_ representations;  (b) it allows for compositional reasoning and,
in particular, in contrast to representations in terms of density or mass
functions, the effective computation of (high-order) moments;  (c) tail
bounds, concentration bounds, and other properties of interest can be
extracted with relative ease from a PGF; and  (d) expressions containing
parameters, both for probabilities and for assigning new values to program
variables, are naturally supported.  Some successfully implemented ideas based
on PGFs, e.g., for deciding probabilistic equivalence and for proving non-
almost-sure termination, are presented in [5, 13], which address especially
the aforementioned challenges (1) and (2) for exact probabilistic inference
_without_ conditioning.

##  2 Taming Conditioning Using PGFs

The creation of generative models is a challenging task, as these models
oftentimes need expert domain knowledge. Therefore, the concept of
_conditioning as a first-class language element_ is crucial as it allows for a
natural and intuitive approach to the creation of models. Our current research
aims to _extend the PGF approach towards exact inference for probabilistic
programs with conditioning – thus addressing challenges(1), (2), and (3) – and
to push the limits of automation as far as possible_. To this end, we are in
the process of developing an exact, symbolic inference engine based on the
open-source, PGF-based tool Prodigy [5]. We illustrate below its current
capability to cater for conditioning via two examples.

$\displaystyle\\{{w\coloneqq{0}}\\}[{\nicefrac{{5}}{{7}}}]\\{{w\coloneqq{1}}\\}\fatsemi$
$\displaystyle\textup{{{if}}}({w=0})\\{\ c\coloneqq{\mathtt{poisson}({6})}\
\\}$ $\displaystyle\textup{{{else}}}\ \\{\ c\coloneqq{\mathtt{poisson}({2})}\
\\}\fatsemi$ $\displaystyle\textup{{{observe}}}({c=5})$
Prog. 1: The telephone operator. $\displaystyle x\coloneqq{1}\fatsemi$
$\displaystyle\textup{{{while}}}({x=1})\\{$
$\displaystyle\quad\\{{c\coloneqq{c+1}}\\}[{\nicefrac{{1}}{{2}}}]\\{{x\coloneqq{0}}\\}\fatsemi\
\\}$ $\displaystyle\textup{{{observe}}}({c\equiv 1\\!\\!\\!\\!\pmod{2}})$
Prog. 2: The odd geometric distribution.

### Conditioning in Loop-Free Programs.

Prog. 1 is a loop-free probabilistic program encoding an infinite-support
distribution. It describes a telephone operator who is unaware of whether
today is a weekday or weekend. The operator’s initial belief is that with
probability $\nicefrac{{5}}{{7}}$ it is a weekday ($w=0$) and thus with
probability $\nicefrac{{2}}{{7}}$ weekend ($w=1$). Usually, on weekdays there
are 6 incoming calls per hour on average; on weekends this rate decreases to 2
calls – both rates are subject to a Poisson distribution. The operator
observes 5 calls in the last hour. The inference task is to compute the
distribution in which the initial belief is updated based on the posterior
observation. Prodigy can automatically infer that
$\textup{Pr}(w=0)=\frac{1215}{1215+2\cdot e^{4}}\approx 0.9178$.

### Conditioning Outside of Loops.

Prog. 2 describes an iterative algorithm that repeatedly flips a fair coin –
while counting the number of trials – until seeing heads, and observes that
this number is odd. Whereas Prog. 1 can be handled by ($\lambda$)PSI [7, 8],
Prog. 2 cannot, as ($\lambda$)PSI do not support programs with unbounded
looping behaviors. However, given a suitable invariant as described in [5],
Prodigy is able to reason about the posterior distribution of Prog. 2 in an
automated fashion using the _second-order PGF_ (SOP) technique [5]: the
resulting posterior distribution for any input with $c=0$ is $\frac{3\cdot
c}{(4-c^{2})}$ which encodes precisely a closed-form solution for the
generating function $\sum_{n=0}^{\infty}-3\cdot
c^{n}\cdot\left(2^{-2-n}\cdot\left(-1+(-1)^{n}\right)\right)$.

##  3 Future Directions

A natural question is whether we can tackle exact inference when conditioning
occurs _inside_ of a loop. As argued in [17], more advanced inference
techniques are required to answer this question. In fact, to the best of our
knowledge, there is no (semi-)automated exact inference technique that allows
for the presence of observe statements inside a (possibly unbounded) loop (an
exception could be the potentially automatable conditional weakest
preexpectation calculus [17]). This is precisely our current research focus.
One promising idea is to develop a non-trivial syntactic restriction of the
programming language, where the more advanced SOP technique [5] can be
generalized to address conditioning inside loops.

The possibility to incorporate symbolic parameters in PGF representations can
enable the application of well-established optimization methods, e.g.,
maximum-likelihood estimations and parameter fitting, to the inference for
probabilistic programs. Other interesting future directions include deciding
equivalence of probabilistic programs with conditioning, amending our method
to continuous distributions using characteristic functions, and exploring the
potential of PGFs in differentiable programming.

## References

  * [1] Ezio Bartocci, Laura Kovács, and Miroslav Stankovic.  Analysis of Bayesian networks via prob-solvable loops.  In ICTAC, volume 12545 of LNCS, pages 221–241. Springer, 2020.
  * [2] Ezio Bartocci, Laura Kovács, and Miroslav Stankovic.  Mora \- Automatic generation of moment-based invariants.  In TACAS (1), volume 12078 of LNCS, pages 492–498. Springer, 2020.
  * [3] Sooraj Bhat, Ashish Agarwal, Richard W. Vuduc, and Alexander G. Gray.  A type theory for probability density functions.  In POPL, pages 545–556. ACM, 2012.
  * [4] Sooraj Bhat, Johannes Borgström, Andrew D. Gordon, and Claudio V. Russo.  Deriving probability density functions from probabilistic functional programs.  Log. Methods Comput. Sci., 13(2), 2017.
  * [5] Mingshuai Chen, Joost-Pieter Katoen, Lutz Klinkenberg, and Tobias Winkler.  Does a program yield the right distribution? Verifying probabilistic programs via generating functions.  In CAV (1), volume 13371 of LNCS, pages 79–101. Springer, 2022.
  * [6] Gregory F. Cooper.  The computational complexity of probabilistic inference using Bayesian belief networks.  Artif. Intell., 42(2-3):393–405, 1990.
  * [7] Timon Gehr, Sasa Misailovic, and Martin T. Vechev.  PSI: Exact symbolic inference for probabilistic programs.  In CAV (1), volume 9779 of LNCS, pages 62–83. Springer, 2016.
  * [8] Timon Gehr, Samuel Steffen, and Martin T. Vechev.  $\lambda$PSI: Exact inference for higher-order probabilistic programs.  In PLDI, pages 883–897. ACM, 2020.
  * [9] Bradley Gram-Hansen.  Extending probabilistic programming systems and applying them to real-world simulators.  PhD thesis, University of Oxford, 2021.
  * [10] Steven Holtzen, Guy Van den Broeck, and Todd D. Millstein.  Scaling exact inference for discrete probabilistic programs.  Proc. ACM Program. Lang., 4(OOPSLA):140:1–140:31, 2020.
  * [11] Zixin Huang, Saikat Dutta, and Sasa Misailovic.  AQUA: Automated quantized inference for probabilistic programs.  In ATVA, volume 12971 of LNCS, pages 229–246. Springer, 2021.
  * [12] Benjamin Lucien Kaminski, Joost-Pieter Katoen, and Christoph Matheja.  On the hardness of analyzing probabilistic programs.  Acta Inform., 56(3):255–285, 2019.
  * [13] Lutz Klinkenberg, Kevin Batz, Benjamin Lucien Kaminski, Joost-Pieter Katoen, Joshua Moerman, and Tobias Winkler.  Generating functions for probabilistic programs.  In LOPSTR, volume 12561 of LNCS, pages 231–248. Springer, 2020.
  * [14] Dexter Kozen.  Semantics of probabilistic programs.  J. Comput. Syst. Sci., 22(3):328–350, 1981.
  * [15] Nicholas Metropolis and Stanisław Ulam.  The Monte Carlo method.  J. Am. Stat. Assoc., 44(247):335–341, 1949.
  * [16] Praveen Narayanan, Jacques Carette, Wren Romano, Chung-chieh Shan, and Robert Zinkov.  Probabilistic inference by program transformation in Hakaru (system description).  In FLOPS, volume 9613 of LNCS, pages 62–79. Springer, 2016\.
  * [17] Federico Olmedo, Friedrich Gretz, Nils Jansen, Benjamin Lucien Kaminski, Joost-Pieter Katoen, and Annabelle McIver.  Conditioning in probabilistic programming.  ACM Trans. Program. Lang. Syst., 40(1):4:1–4:50, 2018.
  * [18] Dan Roth.  On the hardness of approximate reasoning.  Artif. Intell., 82(1):273–302, 1996.
  * [19] Herbert S Wilf.  Generatingfunctionology.  CRC press, 2005.


# Elastic properties of three-dimensional Yukawa or dust crystals from
molecular dynamic simulations

Sandeep Kumar  sandeepshukla1112<EMAIL_ADDRESS>Department of
Physics, University of South Florida, Tampa, Florida 33620, USA

###### Abstract

I calculate elastic properties of three-dimensional Yukawa or dust crystals
using molecular dynamic simulations. The elastic properties are computed by
deforming (compressing/expanding) a cubic dust crystal along different
directions. The bulk modulus, shear modulus, and Poisson’s ratio are
determined as a function of shielding parameter $\kappa$ and strong coupling
parameter $\Gamma$. Furthermore, I calculate sound wave velocity and shear
wave velocity using bulk and shear modulus values and compare them with the
velocities obtained from the slope of longitudinal and transverse wave
dispersion relations. It is found that the sound and shear velocities obtained
using bulk and shear modulus agree with the previously observed experimental
values.

##  I Introduction

Dusty plasmas are composed of electrons, ions, neutrals, and immersed dust
particles. It is observed in both natural and laboratory settings. In nature,
dusty plasmas are observed in Saturn’s ring, interstellar clouds, and cometary
tails Shukla and Mamun (2015). In laboratory, it is present in the plasmas of
fusion devices, rocket exhaust, and created in experimental devices under
controlled conditions Chu and Lin (1994); Chaubey and Goree (2022). The dusty
plasma offers a model system to study generic phenomena such as self-
organization and transport at particle level Morfill and Ivlev (2009). Both
theoretical and experimental studies have been carried out to study generic
phenomena in dusty plasmas, such as crystallization Thomas  _et al._ (1994);
Hamaguchi  _et al._ (1997); Maity and Das (2019), single particle dynamics
Maity  _et al._ (2018); Deshwal  _et al._ (2022), solitons Kumar  _et al._
(2017); Bandyopadhyay  _et al._ (2008); Donkó  _et al._ (2020); Kumar and
Sharma (2023), shocks Sharma  _et al._ (2016); Arora and Maity (2023), spiral
waves Kumar  _et al._ (2018); Kumar and Das (2018), vortices Bailung  _et al._
(2020); Dharodi (2020), and Mach cones Samsonov  _et al._ (1999).

A typical micron-sized dust particle carries an electronic charge ranging from
approximately -10,000e to -20,000e and possesses a mass approximately
$10^{13}$ to $10^{14}$ times that of ions. Dusty plasma can be modeled by a
system of point particles interacting through a Yukawa potential, also known
as a shielded Coulomb potential, described by the following expression Konopka
_et al._ (2000):

$\displaystyle{U(r)=\frac{Q^{2}}{4\pi\epsilon_{0}r}\exp(-\frac{r}{\lambda_{D}})}\,,$
(1)

where $Q=-Z_{d}e$ is the charge on a typical dust particle, $r$ is the
separation between two dust particles, and $\lambda_{D}$ is the Debye length
of background plasma. The Yukawa system can be characterized in terms of two
dimensionless parameters $\Gamma={Q^{2}}/{4\pi\epsilon_{0}ak_{B}T_{d}}$ (known
as the strong coupling parameter) and $\kappa={a}/{\lambda_{D}}$ (known as the
shielding parameter). Here, $T_{d}$ and $a$ are the dust temperature and the
Wigner-Seitz (WS) radius, respectively. The Yukawa inter-particle interaction
is also used to model other systems, including charged colloids Palberg  _et
al._ (1995); Terao and Nakayama (1999), electrolytes Levin (2002); Lee  _et
al._ (2017), and strongly coupled electron-ion plasmas Lyon  _et al._ (2013);
Vorberger  _et al._ (2012).

Elastic properties are crucial material characteristics that indicate the
stiffness of a material. They are linked to the propagation of waves and their
velocities along specific crystallographic directions within a material. Wave
propagation and viscoelasticity in dusty plasmas have been extensively
investigated through experiments and simulations over the past few decades
Morfill and Ivlev (2009). Earlier studies have been carried out to calculate
the elastic properties of Yukawa systems. Robbins  _et al._ (1988) and Khrapak
and Klumov (2020) have conducted simulations where sound velocities and radial
distribution functions have been utilized in the calculations of elastic
properties of Yukawa systems, respectively. Wang  _et al._ (2019) have
determined the shear modulus of two-dimensional dusty plasmas by considering
the viscoelasticity of the liquid state. More recently, Kozhberov (2022) have
studied the elastic properties of Yukawa crystals using analytical
expressions, incorporating the electrostatic energies of deformed crystals in
the calculations. The investigation into the elastic properties of Yukawa or
dust crystals continues to be an area of sustained interest.

In this paper, I calculate the elastic properties of dusty plasma by deforming
(compressing/expanding) a cubic dust crystal along different directions. The
stress and strain of the deformed crystal are used in the calculation of the
elastic stiffness constants of dusty plasma. Furthermore, these elastic
constants are utilized to calculate the bulk modulus, shear modulus, and
Poisson’s ratio of the Yukawa crystal. I determine the bulk modulus, shear
modulus, and Poisson’s ratio as functions of the shielding parameter $\kappa$
and the strong coupling parameter $\Gamma$. These elastic properties are
calculated using experimental dusty plasma parameters. I assess the accuracy
of my elastic property calculations by comparing the sound velocity and shear
velocity obtained using it with those obtained from the slope of longitudinal
and transverse wave dispersion relations.

The paper is organized as follows. Section II provides details of molecular
dynamic (MD) simulations. In Section III, I report the results for elastic
properties as a function of strong coupling parameter $\Gamma$ and shielding
parameter $\kappa$. Section IV contains a summary of the work.

##  II Simulation Method

A three-dimensional cubic box containing point particles is created to study
the elastic properties of dust crystals. Yukawa interaction potential is taken
among the dust particles, which mimics the screening from background electrons
and ions. The state-of-the-art open-source code LAMMPS Thompson  _et al._
(2022) is used for the molecular dynamic simulations. The simulation box
contains 10235 dust particles within dimensions of
35$a$$\times$35$a$$\times$35$a$ (ranging from 0 to 35$a$) along the X, Y, and
Z directions, respectively, where $a$ is average inter-particle distance.
Experimental parameters Bailung  _et al._ (2020) are employed in the
simulations: $Q$ = -20000e (where e denotes the charge of an electron), dust
mass m = 1.7$\times$$10^{-13}$ Kg, and average inter-particle distance $a$ =
6$\times$$10^{-4}$ m. The average inter-particle distance sets the number
density of dust particles $n$ = 1.1$\times$$10^{9}$ per $\text{m}^{3}$,
corresponding to a mass density $\rho$ of 1.87$\times$$10^{-4}$
$\text{Kg}/\text{m}^{3}$. In the simulations, the value of shielding parameter
$\kappa$ is increased by decreasing the screening length $\lambda_{D}$ of the
Yukawa interaction, and the value of strong coupling parameter $\Gamma$ is
increased by decreasing the temperature of dust particles.

The particles are equilibrated at a given temperature using Langevin dynamics,
which reads as

$m\ddot{\mathbf{r}_{\mathrm{i}}}=-\sum_{j}\mathbf{\nabla}U_{ij}+\mathbf{F}_{f}+\mathbf{F}_{r}\,,$
(2)

where $U_{ij}$ is interaction potential, $\mathbf{F}_{f}$ is frictional force
on the particles, and $\mathbf{F}_{r}$ is random force (kicks) on the
particles. I choose simulation time step of 100 $\mu$s, which ensures a fine
discretization along the temporal domain and good resolution of the underlying
dust kinetics. The dust particles are evolved for the 40000 time steps for the
equilibration of the system. Fluctuations in temperature, pressure, and total
energy of the system over time are monitored to verify system equilibration.

Following equilibration, the system is prepared for the calculation of elastic
properties under isothermal conditions. The cubic box is transformed into a
triclinic box and equilibrated again for 450000 time steps. Subsequently, the
triclinic box undergoes deformation (Voigt deformation) along various
directions. After each expansion and contraction, the triclinic box is
equilibrated for 50000 time steps. The change in the stress of crystal
$P_{x}$, $P_{y}$, $P_{z}$, $P_{xy}$, $P_{xz}$, $P_{yz}$ with deformation
(strain) along X, Y, Z, XY, XZ, and YZ directions (Voigt deformation
component), respectively, are used in the calculations of elastic stiffness
constants.

The convergence of elastic properties with magnitude of deformation
(expansion/compression) is checked carefully. A deformation magnitude of 0.075
(in strain units) is employed in my calculations. The average of stresses are
taken at 1000 time steps. The finite size effects on elastic properties are
examined, which reveals that 10235 atoms are adequate to mitigate such
effects.

##  III Results

The purpose of this study is to calculate elastic properties of three-
dimensional Yukawa or dust crystal as a function of $\kappa$ and $\Gamma$.
Specifically, I compute the bulk modulus, shear modulus, and Poisson’s ratio
to characterize the elastic properties.

###  III.1 Bulk Modulus

The bulk modulus $K$ is a measure of the resistance of a material to an
applied bulk compression. I calculate bulk modulus using the following
expression Kittel and McEuen (2018); Robbins  _et al._ (1988):

$K=\frac{C_{11}+2C_{12}}{3}\,.$ (3)

In the above expression, $C_{11}$ and $C_{12}$ are the elastic stiffness
constants. The bulk modulus at $\Gamma$ = 2000 as a function of screening
parameter $\kappa$ is displayed in Fig. 1. It decreases with an increase in
the $\kappa$. This reduction occurs because as $\kappa$ increases, the
interaction potential length decreases, making compression of the material
easier. For large $\kappa$, the Yukawa potential becomes extremely short-
range, leading to interactions among particles resembling those in a hard
sphere system. The characteristics of bulk modulus at $\kappa$ = 2 with
varying strong coupling parameter $\Gamma$ is shown in Fig. 2. The bulk
modulus decreases with an increase in $\Gamma$. In this study, $\Gamma$ is
increased by reducing the temperature of the dust particles, making it easier
to compress the particles due to reduced thermal pressure. However, at higher
$\Gamma$ values, the rigidity driven by the interaction potential overtakes
the rigidity driven by temperature, resulting in no further change in bulk
modulus with decreasing temperature (see Fig. 2).

Figure 1: Bulk modulus $K$ of dust crystal at the strong coupling parameter
$\Gamma$ = $2000$ as a function of shielding parameter $\kappa$. Figure 2:
Bulk modulus $K$ of dust crystal at shielding parameter $\kappa$ = $2$ as a
function of strong coupling parameter $\Gamma$.

To validate the accuracy of my bulk modulus calculations, I compute the
longitudinal sound velocity $C_{S}$ using it. I utilize the relationship
$C_{S}$ = $\sqrt{K/\rho}$ for the calculation of sound velocity Kittel and
McEuen (2018). At $\kappa$ = 2.0 and $\Gamma$ = 3000, the sound velocity
obtained using bulk modulus is 2.01 cm/s, which closely matches the sound
velocity value of 2.15 cm/s obtained from the slope of the longitudinal wave
dispersion relation, see Fig. S1 in supplementary material. The details of the
sound velocity calculation from the dispersion relation are provided in the
supplementary material sup  and have also been recently reported by Kumar  _et
al._ (2023). The experimentally measured and theoretically calculated values
of sound velocity by Bailung  _et al._ (2020) at the same parameters are 4.5
cm/s and 3.06 cm/s, respectively. From these comparisons, it is evident that
my bulk modulus values yield correct sound velocity in the dusty plasma
medium. The strongly coupled dusty plasmas are extremely soft, so that the
bulk modulus of solid dusty plasma is much smaller than those of typical
solids like metals. For example, bulk modulus of aluminum at ambient condition
is 79 GPa Simmons (1965).

###  III.2 Shear Modulus

The shear modulus $G$ is a measure of a material’s resistance to an applied
shearing deformation. I calculate the shear modulus using the following
formulas Kittel and McEuen (2018); Robbins  _et al._ (1988):

$G_{1}=C_{44}$ (4)

and

$G_{2}=\frac{C_{11}-C_{12}}{2}\,.$ (5)

Here, $G_{1}$ and $G_{2}$ are the shear modulus along 100 and 110
crystallographic directions, respectively. The shear modulus $G_{1}$ as a
function of $\Gamma$ at different $\kappa$ values is displayed in Fig. 3.
Additionally, the shear modulus as a function of $\kappa$ at $\Gamma$ = 2000
along the 100 and 110 crystallographic directions are shown in Fig. 4. The
values of $G_{1}$ and $G_{2}$ differ at small $\kappa$ values because the
particles are in a crystalline phase, resulting in different particle
arrangements along the 100 and 110 directions, leading to different shear
strains along these directions. However, at high $\kappa$ values, the
particles are in a liquid phase and exhibit a uniform distribution along all
directions resulting in $G_{1}$ and $G_{2}$ are equal. In the crystalline
phase, the shear modulus decreases with an increase in the shielding
parameter, which is a direct consequence of the reduction in inter-particle
interactions. The shear modulus initially decreases with an increase in the
strong coupling parameter $\Gamma$, and later it increases with further
increases in the strong coupling parameter due to an increase in inter-
particle interactions. This increase in inter-particle interaction results in
an increase in resistance to an applied shear deformation. At very high
$\Gamma$ values, temperature-driven shear deformation becomes very low and is
suppressed by inter-particle driven shear resistance, consequently causing the
shear modulus to saturate (see Fig. 3).

Figure 3: Shear modulus $G_{1}$ of dust crystal at different shielding
parameters $\kappa$ as a function of strong coupling parameter $\Gamma$.
Figure 4: Shear modulus of dust crystal along crystallographic directions of
100 and 110 as a function of shielding parameter $\kappa$. The strong coupling
parameter $\Gamma$ of the Yukawa crystal is 2000.

At $\kappa$ = 2.0 and $\Gamma$ = 3000, the shear velocity obtained from the
shear modulus relation $C_{T}$ = $\sqrt{G/\rho}$ Kittel and McEuen (2018) is
3.61 mm/s and 3.59 mm/s along the 100 and 110 directions, respectively. These
values are close to the experimentally measured value of 4.2 mm/s reported by
Pramanik  _et al._ (2002). Nunomura  _et al._ (2000) have also reported
similar values of shear velocity in their experiments. From the slope of the
transverse wave dispersion relation, the obtained shear velocity is 5.36 mm/s,
see Fig. S2 in the supplementary material sup . Nevertheless, my shear modulus
values are providing correct shear velocity in the dusty plasma medium.
Similar to the bulk modulus, the shear modulus is much smaller than that of
typical solids. For instance, the shear modulus for aluminum at ambient
conditions is 25.5 GPa Archer  _et al._ (1972).

###  III.3 Poisson’s Ratio

The Poisson’s ratio $\nu$ quantifies the deformation (compression/expansion)
of a material in directions perpendicular to a specific direction of loading.
It is defined as the negative ratio of transverse strain to axial strain. I
calculate Poisson’s ratio using the elastic constants, which is given by

$\nu=\frac{1}{1+\frac{C_{11}}{C_{12}}}\,.$ (6)

The Poisson’s ratio at $\Gamma$ = 2000 as a function of $\kappa$ is displayed
in Fig. 5. Initially, Poisson’s ratio decreases with an increase in $\kappa$
because transverse strain decreases with the reduction in interaction length
of potential. However, above $\kappa$ > 2.5, it increases because the Yukawa
interaction becomes extremely short-range, resulting in interactions among
particles resembling those in a hard sphere system that leads to much less
axial strain. The characteristics of Poisson’s ratio at $\kappa$ = 2 as a
function of $\Gamma$ is shown in Fig. 6. With an increase in the strong
coupling parameter $\Gamma$, Poisson’s ratio decreases, which is a direct
consequence of less strain in the transverse direction due to the reduced
thermal motion (thermal pressure) of particles at low temperatures.

Figure 5: Poisson’s ratio of dust crystal at coupling parameter $\Gamma$ =
$2000$ as a function of shielding parameter $\kappa$. Figure 6: Poisson’s
ratio of dust crystal at shielding parameter $\kappa$ = 2 as a function of
coupling parameter $\Gamma$.

For dusty plasma, the bulk modulus is much larger compared to the shear
modulus, so they can be regarded as effectively incompressible, as it is
easier to change shape than to compress. Soft materials with similar
characteristics are also considered incompressible in nature Jastrzebski
(1976). A perfectly incompressible isotropic material has a Poisson’s ratio of
exactly 0.5; for dusty plasma, it is approximately 0.5. However, the
incompressible nature of dusty plasma changes with variations in $\kappa$ and
$\Gamma$, see Fig. 5 and Fig. 6. For comparison, it is worth mentioning that
the Poisson’s ratio of natural rubber and aluminum is 0.4999 and 0.33,
respectively Mott and Roland (2009).

##  IV Summary

I have calculated bulk modulus, shear modulus, and Poisson’s ratio of a cubic
Yukawa or dusty crystal in the molecular dynamic simulations. These properties
are computed by deforming (compressing/expanding) the dust crystal along
different directions. The characteristics of elastic properties as a function
of shielding parameter and strong coupling parameter have been studied. The
shear modulus has also been calculated along the crystallographic directions
of 100 and 110, and it is found that the values are different at small
screening parameters. Furthermore, I have calculated sound and shear wave
velocities using bulk and shear modulus and compared them with the velocities
obtained from the slope of longitudinal and transverse wave dispersion
relations. It is found that the sound and shear velocities obtained using bulk
and shear modulus agree with the previously observed experimental values.

Strongly coupled dusty plasmas exhibit solid-like traits, but they are
extremely soft. As a result, the bulk modulus and shear modulus of solid dusty
plasmas are much smaller than those of typical solids, such as metals. In
fact, the bulk modulus is much larger compared to the shear modulus, which
effectively renders dusty plasma as incompressible. This is because it is
easier to change shape than to compress a dusty plasma system.

###### Acknowledgements.

The author thanks Dr. Attila Cangi (CASUS-HZDR, Görlitz, Germany) and Dr.
Svetoslav Nikolov (SNL, Albuquerque, USA) for helpful discussions. The
computations were performed on a Bull Cluster at the Center for Information
Services and High Performance Computing (ZIH) at Technische Universität
Dresden, on the cluster Hemera at Helmholtz-Zentrum Dresden-Rossendorf (HZDR).

## References

  * Shukla and Mamun (2015) Padma K Shukla and AA Mamun,  _Introduction to dusty plasma physics_ (CRC press, 2015).
  * Chu and Lin (1994) JH Chu and I Lin, “Direct observation of coulomb crystals and liquids in strongly coupled rf dusty plasmas,” Phys. Rev. Lett. 72, 4009 (1994).
  * Chaubey and Goree (2022) Neeraj Chaubey and J Goree, “Coulomb expansion of a thin dust cloud observed experimentally under afterglow plasma conditions,” Phys. Plasmas 29 (2022).
  * Morfill and Ivlev (2009) Gregor E Morfill and Alexei V Ivlev, “Complex plasmas: An interdisciplinary research field,” Rev. Mod. Phys. 81, 1353 (2009).
  * Thomas  _et al._ (1994) H. Thomas, G. E. Morfill, V. Demmel, J. Goree, B. Feuerbacher, and D. Möhlmann, “Plasma crystal: Coulomb crystallization in a dusty plasma,” Phys. Rev. Lett. 73, 652–655 (1994).
  * Hamaguchi  _et al._ (1997) S Hamaguchi, RT Farouki, and DHE Dubin, “Triple point of yukawa systems,” Phys. Rev. E 56, 4671 (1997).
  * Maity and Das (2019) Srimanta Maity and Amita Das, “Molecular dynamics study of crystal formation and structural phase transition in yukawa system for dusty plasma medium,” Phys. Plasmas 26 (2019).
  * Maity  _et al._ (2018) Srimanta Maity, Amita Das, Sandeep Kumar, and Sanat Kumar Tiwari, “Interplay of single particle and collective response in molecular dynamics simulation of dusty plasma system,” Phys. Plasmas 25 (2018).
  * Deshwal  _et al._ (2022) Priya Deshwal, Mamta Yadav, Chaitanya Prasad, Shantam Sridev, Yash Ahuja, Srimanta Maity, and Amita Das, “Chaotic dynamics of small-sized charged yukawa dust clusters,” Chaos 32 (2022).
  * Kumar  _et al._ (2017) Sandeep Kumar, Sanat Kumar Tiwari, and Amita Das, “Observation of the korteweg-de vries soliton in molecular dynamics simulations of a dusty plasma medium,” Phys. Plasmas 24, 033711 (2017).
  * Bandyopadhyay  _et al._ (2008) P Bandyopadhyay, G Prasad, A Sen, and PK Kaw, “Experimental study of nonlinear dust acoustic solitary waves in a dusty plasma,” Phys. Rev. Lett. 101, 065006 (2008).
  * Donkó  _et al._ (2020) Zoltán Donkó, Peter Hartmann, Ranna U Masheyeva, and Karlygash N Dzhumagulova, “Molecular dynamics investigation of soliton propagation in a two-dimensional yukawa liquid,” Contrib. Plasma Phys. 60, e201900197 (2020).
  * Kumar and Sharma (2023) Prince Kumar and Devendra Sharma, “Quasi-localized charge approximation approach for the nonlinear structures in strongly coupled yukawa systems,” Phys. Plasmas 30 (2023).
  * Sharma  _et al._ (2016) Sumita K Sharma, A Boruah, Y Nakamura, and H Bailung, “Observation of dust acoustic shock wave in a strongly coupled dusty plasma,” Phys. Plasmas 23 (2016).
  * Arora and Maity (2023) Garima Arora and Srimanta Maity, “Self-excited converging shock structure in a complex plasma medium,” Phys. Rev. E 108, 055209 (2023).
  * Kumar  _et al._ (2018) Sandeep Kumar, Bhavesh Patel, and Amita Das, “Spiral waves in driven dusty plasma medium: Generalized hydrodynamic fluid description,” Phys. Plasmas 25, 043701 (2018).
  * Kumar and Das (2018) Sandeep Kumar and Amita Das, “Spiral waves in driven strongly coupled yukawa systems,” Phys. Rev. E 97, 063202 (2018).
  * Bailung  _et al._ (2020) Yoshiko Bailung, Bidyut Chutia, T Deka, A Boruah, Sumita K Sharma, Sandeep Kumar, Joyanti Chutia, Y Nakamura, and H Bailung, “Vortex formation in a strongly coupled dusty plasma flow past an obstacle,” Phys. Plasmas 27 (2020).
  * Dharodi (2020) Vikram S Dharodi, “Rotating vortices in two-dimensional inhomogeneous strongly coupled dusty plasmas: Shear and spiral density waves,” Phys. Rev. E 102, 043216 (2020).
  * Samsonov  _et al._ (1999) D Samsonov, J Goree, ZW Ma, A Bhattacharjee, HM Thomas, and GE Morfill, “Mach cones in a coulomb lattice and a dusty plasma,” Phys. Rev. Lett. 83, 3649 (1999).
  * Konopka  _et al._ (2000) U. Konopka, G. E. Morfill, and L. Ratke, “Measurement of the interaction potential of microspheres in the sheath of a rf discharge,” Phys. Rev. Lett. 84, 891–894 (2000).
  * Palberg  _et al._ (1995) T. Palberg, W. Mönch, F. Bitzer, R. Piazza, and T. Bellini, “Freezing transition for colloids with adjustable charge: A test of charge renormalization,” Phys. Rev. Lett. 74, 4555–4558 (1995).
  * Terao and Nakayama (1999) Takamichi Terao and Tsuneyoshi Nakayama, “Crystallization in quasi-two-dimensional colloidal systems at an air-water interface,” Phys. Rev. E 60, 7157–7162 (1999).
  * Levin (2002) Yan Levin, “Electrostatic correlations: from plasma to biology,” Rep. Prog. Phys. 65, 1577 (2002).
  * Lee  _et al._ (2017) Alpha A. Lee, Carla S. Perez-Martinez, Alexander M. Smith, and Susan Perkin, “Scaling analysis of the screening length in concentrated electrolytes,” Phys. Rev. Lett. 119, 026002 (2017).
  * Lyon  _et al._ (2013) M. Lyon, S. D. Bergeson, and M. S. Murillo, “Limit of strong ion coupling due to electron shielding,” Phys. Rev. E 87, 033101 (2013).
  * Vorberger  _et al._ (2012) Jan Vorberger, Z Donko, IM Tkachenko, and Dirk O Gericke, “Dynamic ion structure factor of warm dense matter,” Phys. Rev. Lett. 109, 225001 (2012).
  * Robbins  _et al._ (1988) Mark O Robbins, Kurt Kremer, and Gary S Grest, “Phase diagram and dynamics of yukawa systems,” J. Chem. Phys. 88, 3286–3312 (1988).
  * Khrapak and Klumov (2020) Sergey A Khrapak and Boris A Klumov, “Instantaneous shear modulus of yukawa fluids across coupling regimes,” Phys. Plasmas 27 (2020).
  * Wang  _et al._ (2019) Kang Wang, Dong Huang, and Yan Feng, “Shear modulus of two-dimensional yukawa or dusty-plasma solids obtained from the viscoelasticity in the liquid state,” Phys. Rev. E 99, 063206 (2019).
  * Kozhberov (2022) AA Kozhberov, “Elastic properties of yukawa crystals,” Phys. Plasmas 29 (2022).
  * Thompson  _et al._ (2022) Aidan P Thompson, H Metin Aktulga, Richard Berger, Dan S Bolintineanu, W Michael Brown, Paul S Crozier, Pieter J In’t Veld, Axel Kohlmeyer, Stan G Moore, Trung Dac Nguyen, _et al._ , “Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comput. Phys. Commun. 271, 108171 (2022).
  * Kittel and McEuen (2018) Charles Kittel and Paul McEuen,  _Introduction to solid state physics_ (John Wiley & Sons, 2018).
  * (34) See Supplemental Material at [URL will be inserted by publisher].
  * Kumar  _et al._ (2023) Sandeep Kumar, Hossein Tahmasbi, Kushal Ramakrishna, Mani Lokamani, Svetoslav Nikolov, Julien Tranchida, Mitchell A. Wood, and Attila Cangi, “Transferable interatomic potential for aluminum from ambient conditions to warm dense matter,” Phys. Rev. Res. 5, 033162 (2023).
  * Simmons (1965) Gene Simmons, “Single crystal elastic constants and calculated aggregate properties,” J. Grad. Res. 34, 1 (1965).
  * Pramanik  _et al._ (2002) J Pramanik, G Prasad, A Sen, and PK Kaw, “Experimental observations of transverse shear waves in strongly coupled dusty plasmas,” Phys. Rev. Lett. 88, 175001 (2002).
  * Nunomura  _et al._ (2000) S Nunomura, D Samsonov, and J Goree, “Transverse waves in a two-dimensional screened-coulomb crystal (dusty plasma),” Phys. Rev. Lett. 84, 5141 (2000).
  * Archer  _et al._ (1972) Robert R. Archer, Stephen H. Crandall, Norman C. Dahl, and Thomas J. Lardner,  _An introduction to the mechanics of solids_ , 2nd ed. (McGraw-Hill Kogakusha, 1972).
  * Jastrzebski (1976) Zbigniew D Jastrzebski,  _Nature and properties of engineering materials_ (John Wiley and Sons, Inc., New York, 1976).
  * Mott and Roland (2009) PH Mott and CM Roland, “Limits to poisson’s ratio in isotropic materials,” Phys. Rev. B 80, 132104 (2009).