[type=editor, orcid=0000-0003-0488-8513, ]

[1]

[1]

1]organization=National Center for Nuclear Research, Materials Research
Laboratory, addressline=ul. Andrzeja Sołtana 7, city=Otwock-Świerk,
postcode=05-400, country=Poland

[1]Corresponding author

# Dynamics of stress propagation in anharmonic crystals: MD simulations.

Zbigniew Kozioł<EMAIL_ADDRESS>[

###### Abstract

Anharmonic inter-atomic potential $\sim|x|^{n}$, $n>1$, has been used in
molecular dynamics (MD) simulations of stress dynamics of FCC oriented
crystal. The model of chain of masses and springs is found as a convenient and
accurate description of simulation results, with masses representing the
crystallographic planes. The period of oscillations of two planes is found
analytically to be given by Euler’s beta functions, and it’s scaling with non-
linearity parameter and amplitude of oscillations, or applied shear pressure
is discussed on examples of time dependencies of displacement, velocities, and
forces acting on masses (planes). The relationship between physical quantities
(displacement, velocities, pressure) is explained. The dynamics of stress
penetration from the sample surface towards its interior is confirmed to be
given exactly by a series of Bessel functions when n=2 (in harmonic case).
When n $\neq$ 2 the stress dynamics (wave propagation) in bulk material
remains qualitatively of the same nature as in the linear case. In particular,
results suggest that the quasi-linear relationship between frequency and wave
number is preserved. Speed of the transverse sound component, dependent on
sound wave amplitude, is found to be a strongly decreasing function of $n$.
The results are useful in the analysis of any MD simulations under pressure,
as they help to understand the dynamics of pressure retarded effects, as well
help design the proper methodology of performing MD simulations, such as
studies of dynamics of dislocations.

###### keywords:

Molecular dynamics simulations Nonlinear Elastic properties Internal stress
Anharmonic crystals

Table method for creating NN inter-atomic potential $|x|^{n}$ is developed.

Period of oscillations is computed analytically and tested in MD simulations.

Scaling of physical quantities (displacement, velocity, stress) on $n$ is
analyzed.

Profiles of virial stress waves are found and the sound velocity is
determined.

##  1 Introduction

The proper understanding of the dynamics of crystalline solids is a key aspect
of materials science. Tools and techniques for visualization of such dynamics
are in general lacking [1]. As a consequence, materials models have suffered
from a lack of input from experiments. A partial remedy to this is MD
simulations, providing an insight into sub-picosecond time scales and
subatomic spatial resolution.

There is several questions we found when performing MD simulations. For
instance, how to define (strictly, in physical terms) the stress-strain region
when crystal response remains linear? It appears that MD results provide a
description of a dynamic process when a pressure is imposed on a crystal and
there is no simple method developed that would allow to extract linear
response from simulation data. On another hand, how to find out the value of
pressure acting on a dislocation in material, while pressure penetration is a
dynamic process of a comparable time scale as the time of simulation itself?
Additionally, during the pressure front passage through the dislocation, the
interactions between atoms are strongly non-linear, and energy losses occur.
Questions like these inspired us to find a description of the dynamics of
atoms in crystals in a case when a strongly nonlinear, well-defined inter-
atomic potential is used. We choose it in the form $|x|^{n}$, with $n>1$,
since it is given then by a simple analytical function, and in the special
case of $n=2$ it allows a comparison with the most common, classical situation
of harmonic dynamics.

Our first paper on this subject, [2], contains already some ideas about, and
some evidence for, how the stress in MD simulations propagates into the sample
volume. Since then, we concentrated on attempts to understand the physical
mechanism and providing a mathematical description of pressure penetration. At
the same time, we excelled in our technical skills, improved the accuracy of
simulations and data analysis. Seldom do we observe that the retarded nature
of stress propagation is taken into account in simulations, therefore this
aspect, the proper understanding of what is simulated, is of crucial
importance. The whole description of all details that are important in this
kind of MD simulations, and the description of the model we base our analysis
on will be available in a separate paper [3].

The harmonic theory used in the lattice vibrations of solids assumes that the
anharmonic terms in the lattice potential energy expansion are neglected while
the quadratic term is retained. This assumption has several consequences: a)
There is no thermal expansion of solids. b) Interaction between lattice waves
does not exist and a single elastic wave does not decay or change its form
with time. c) Elastic constants do not depend on pressure and temperature.
Hence, understanding the role of anharmonicity is the key factor for
understanding the dynamic and thermodynamic properties of real materials.

Non-linearity in the physics of crystals plays a huge role and it is
omnipresent. Even small, nonlinear contributions to interatomic potentials,
may considerably influence the dynamics of particles, leading to chaotic
motion and may cause seemingly unexpected energy losses [4]. It may cause
localization of vibrations, depending on the strength of the contribution [5].
The first systematic studies of the role of these effects date back to over a
half-century ago with the methods of thermodynamic Green’s functions [6].
Lattice dynamics and conditions of its stability have been addressed as well
[7].

Materials with highly anharmonic potentials are not rare. X-ray diffraction
and terahertz time-domain spectroscopy is used to resolve potential in
molecular crystals [8], or higher harmonics detection in ultrasound
propagation in steel [9] is used, or analysis of phonon oscillation spectra in
graphene and other materials [10]. Higher-order phonon non-linearities are
observed by optical, infrared, and terahertz spectroscopy [11].

In MD simulations most of the potentials available may be treated as
approximately only harmonic, in the best case. For instance, we found in our
studies that Bonny [12] and Artur [13] EAM potentials commonly used in
simulations of steel have a small anharmonic contribution. In the case of
Bonny’s, an anharmonic addition is of the form $|x|^{2.5}$.

Many potentials are explicit anharmonic, like Morse potential [14],
$E(r)\sim(exp(-2\alpha(r-r_{0}))-2exp(-\alpha(r-r_{0})))$, Buckingham
potential: $E(r)\sim(exp(-r/\rho)-C/r^{6})$ [15], or variants of Lennard-Jones
potential [9].

Hence, in this article we study the simplest possible nonlinear inter-atomic
potentials of form $V(x)=\epsilon_{0}x^{n}$, where n>1, mostly it is between
n=1.5 and n=3.60. Since closed form analytical solutions to equations of
motion in that case are not known in general, we perform simulations of
dynamics of particles by using an MD simulation tool, LAMMPS [16]. It is an
overlooked opportunity that MD simulations may be used as a tool for
inspiration and for guiding our understanding of physical processes, largely
playing a role like that of experimental results and it is used often in
studies of lattice dynamics of solids [17].

In section 2 we describe the simulation setup and creation of nonlinear
potentials. Section 3 contains the classical derivation of a period of
oscillations of a two-layer system for two simulation modes considered.
Scaling relations are discussed between physical quantities such as
displacement, velocity, pressure, etc. Section 4 presents MD results for
larger structures with many crystallographic layers. The dynamics of virial
stress profiles (waves of pressure) are discussed and compared to the results
of a theoretical model available for harmonic crystals. Finally, we determine
how the speed of sound is influenced by the non-linearity parameter $n$.

##  2 Potentials, simulations setup, and samples

Figure 1: Distance to NN on the next upper plane. The green atom is located at
the lower plane of these two planes on a). It has 3 NN on the upper plane
(large blue atoms). If to assume that the position of the green atom is
$\vec{a_{0}}$=(0,0,0), we find that (in units of lattice parameter $a$) the
positions of 3 NN (blue) atoms may be written as:
$\vec{a_{1}}=(-1/(2\sqrt{2})$, $1/\sqrt{3}$, $+1/(2\sqrt{2}\sqrt{3}))$,
$\vec{a_{2}}=(0,1/\sqrt{3}$, $-1/(\sqrt{2}\sqrt{3})$ and
$\vec{a_{3}}$=$(+1/(2\sqrt{2})$, $1/\sqrt{3}$, +$1/(2\sqrt{2}\sqrt{3}$)). The
distance between green and blue atoms is $a/\sqrt{2}$. Distance between blue
atoms is also $a/\sqrt{2}$ and it is shown by the red arrow in b). That
distance is the translational period of potential in X-direction.

Since we want to have the simulation conditions similar as much as possible to
the situation we already studied intensively and very accurately, we use a
special configuration suitable when the dynamics of dislocations in FCC steel
is studied.

Hence, we assume the FCC structure of the material oriented in a specific way
[2],[3]: Y direction is along [111] crystallographic axis, X is in
[$1\bar{1}0$] and Z in [$\bar{1}\bar{1}2$]. In dynamical studies either an
initial displacement of a layer of atoms is made in X direction or a pressure
is applied to the upper sample layer in that direction, as shown by the label
"Force" and arrow in Fig. 1 c). Figure 1 a) and b) shows the position of atoms
from two perspectives.

In this orientation of the sample, it follows from symmetries that a possible
in general contribution to forces in the Y direction equals 0 when a shear
force is applied in the X-direction.

In MD simulations we register (with desired time resolution, every few time
steps) quantities like positions, velocities, pressure, stress, etc. Next, we
compute the average values of these quantities for whole separate planes.
Planes are enumerated from 1 at the bottom (lower region in Fig. 1) to N for
the uppermost layer (upper region). The number of atoms in every plane must be
exactly the same, in particular in the case of the regions lower and upper.
Otherwise, simulation results obtained will be noisy and hard for
interpretation. It is important also [3] that the upper region in particular
contains exactly one layer. To achieve these conditions, the sample after
annealing at the desired temperature has the top and bottom layer(s) removed
and the remaining part is checked for the same number of atoms in layers by
using scripts that count atoms. Assigning atoms to layers is done by using
their Y-coordinate. This method works well for samples with up to hundreds of
layers when layer separation distance is properly chosen. Hence, finally, we
must obtain a sample with the number of layers N and the total number of atoms
in the sample must be an exact multiplicity of N. In the case of all samples
studied here, every layer has 4800 atoms (samples have the size of 120
oriented unit cells in the X-direction and 20 in the Z-direction), while N
ranges from 2 to 33.

We use the table pair style in LAMMPS of providing interatomic potential. A
text file is created for that by using Perl scripts. It contains dependence of
force and potential energy as a function of distance between particles. We
aim, in this simple approximation, to account for forces between nearest
neighbors (NN) only. Therefore the potential provided must be limited by an
upper cut-off. We used 2.8 Å for that. Hence, we are also limited to some
maximal values of applied displacements or pressure. The cut-off on the lower
side plays less role in simulations itself, however, it may play a role in
initial sample annealing. We used 0.2 Å as a lower cut-off (such a short cut-
off is not needed, as we found out later). For annealing, it makes sense to
add additional, repulsive forces at short distances. We use for that an extra
contribution of the form $E_{p}(r)=(1/r)^{6}-(1/r_{0})^{6}$, where $r_{0}$=2.2
Å.

A larger outer cutoff distance is not useful in our case. It will cause that
no stable crystallographic structure is created during annealing: introducing
interactions with NNN will smear out the effect of potential minimum. The
reason for that is that large repulsive, largely random forces arise then.
Atoms in these cases tend to granulate in complex, interesting otherwise
structures, but with no long-range order.

We aim to use the crystal structure and sample configuration and properties in
such a way that we could reproduce possibly closely the results for n=2 to
these observed in steel in our previous work. Therefore we choose
$\epsilon_{0}=-4.2$ eV at the minimum of the potential well, as in the case of
the potential of Belland et al. [18] characterized at the NIST111NIST,
Interatomic Potentials Repository repository. To particles, we assign the mass
of one of Fe.

Hence, for the convenience of numerical computation, we use the following
formula for the potential:

$E_{p}(r)=\epsilon_{0}\left(-1+\left|(r-r_{0})\right|^{n}\right),$ (1)

where $r$ and $r_{0}$ are in Å and n could be any positive value larger than
1. Most of the results presented here are for $n$ in the range 1.5-3.5. The
potential is 0 at $r-r_{0}=1$ [Å] and it is $-\epsilon_{0}$ at $r-r_{0}=0$.

Let us concentrate on finding out the potential energy of the green atom in
Fig. 1 a) and b) (all atoms are in the equivalent surroundings). When it is
moved in X-direction for x, the distance between it and 3 blue atoms becomes:

$r_{1}(x)=a\cdot\sqrt{\left(\frac{-1}{2\sqrt{2}}+\frac{x}{a}\right)^{2}+\left(\frac{1}{\sqrt{3}}\right)^{2}+\left(\frac{1}{2\sqrt{2}\sqrt{3}}\right)^{2}}=r_{0}\cdot\sqrt{1-\frac{x}{r_{0}}+\left(\frac{x}{r_{0}}\right)^{2}},$
(2)

$r_{2}(x)=a\cdot\sqrt{\left(\frac{x}{a}\right)^{2}+\left(\frac{1}{\sqrt{3}}\right)^{2}+\left(\frac{-1}{\sqrt{2}\sqrt{3}}\right)^{2}}=r_{0}\cdot\sqrt{1+\left(\frac{x}{r_{0}}\right)^{2}},$
(3)

$r_{3}(x)=a\cdot\sqrt{\left(\frac{1}{2\sqrt{2}}+\frac{x}{a}\right)^{2}+\left(\frac{1}{\sqrt{3}}\right)^{2}+\left(\frac{1}{2\sqrt{2}\sqrt{3}}\right)^{2}}=r_{0}\cdot\sqrt{1+\frac{x}{r_{0}}+\left(\frac{x}{r_{0}}\right)^{2}},$
(4)

where $r_{0}=a/\sqrt{2}$.

Every atom has 12 NN, 6 of them on the same plane and 3 on every of two
neighboring planes. Therefore the total potential energy of an atom ought to
be -12$\epsilon_{0}$. However, we count in this way potential energy twice.
Hence, the potential energy per one atom ought to be 12$\epsilon_{0}$/2=-25.2
eV. Indeed, that is exactly what is reported by LAMMPS.

The potential has been checked in two ways. First, we calculated numerically
E(x) by using Eqs. 1-4 and compared the results with explicit functions of
type $x^{n}$. Next, we did that in LAMMPS, by determining changes in potential
energy when the layers of atoms are displaced relative to each other. In the
range of x between 0 and 0.2 Å, and for n between 1.5 and 3.5, the same
exponent is observed as designed (the least squares fit of potential
determined from MD simulations to function $x^{n}$ gives deviations in $n$
less than 0.1 %), with the overall result described by this kind of function:
$E_{p}(x)=\epsilon_{0}\cdot(x/\xi)^{n}$, where $\epsilon_{0}=8.414$ eV,
$\xi$=2.015 Å, and x is in Å.

As a result of imposed conditions average X for the region lower does not
change with time. Average Y for regions lower and upper remain the same.

##  3 Two-layer system.

###  3.1 Period of oscillations.

The equation of motion for an anharmonic oscillator is nonlinear, making it
more challenging to solve analytically compared to a simple harmonic
oscillator.

One common approach to finding approximate solutions for an anharmonic
oscillator is perturbation theory, where one starts with the solution for a
related harmonic oscillator problem and then adds correction terms to account
for the anharmonicity. The idea is to start with the solutions of the simple
harmonic oscillator and then apply perturbation techniques to account for the
deviation caused by the anharmonicity (nonlinearity) of the potential.

It’s important to note that the higher the value of $n$ and the higher the
order of the perturbation theory considered, the more complex the calculations
become. In some cases, the perturbation series may not converge well, leading
to difficulties in obtaining accurate results.

Another method is to use numerical techniques, such as solving the
differential equation of motion numerically using computer algorithms like the
iterative Runge-Kutta method (see, e.g. [19]).

Analytical solutions for anharmonic oscillators with potentials of that form
are not straightforward and, in general, do not have simple closed-form
expressions like those found in the case of the simple harmonic oscillator.

Finding exact, general solutions of potentials of that form, and approximate
ones, have been discussed by Euler [20], Amore [21], Harko [22]. They use the
idea of performing a non-linear transformation of variables x and t. The
general solution of a class of equations of that type can be expressed as an
infinite sum of components with the hypergeometric function [22]. That however
lacks the simplicity to be analyzed.

Approximate solutions are found for a logarithmic form of anharmonicity, which
is shown to be related to conventional power-law anharmonicity $|x|^{n}$ in
the limit $n\rightarrow 0$, [23].

Variants of the chain of masses and springs model have been used often to
study aspects of non-linearities, chaos, localization/de-localization of
lattice vibrations [9], [24], [25], [26]. The model gain a broad interest
after the Fermi-Pasta-Ulam paradox has been published in 1955 [27].

There is another way of gaining an understanding of solutions to the problem:
molecular dynamics simulations. It is flexible, allows one to approach
problems more and more accurately by introducing or removing some physical
interactions, approximating them gradually to the real.

From the equation of motion, $m\ddot{x}=-dV(x)/dx$, we obtain an integral of
the motion which is the total energy:

$E=\frac{m\dot{x}^{2}}{2}+V(x)$ (5)

The motion of the particle is restricted to the interval $x_{-}<x<x_{+}$,
where the turning points $x_{\pm}$ satisfy $V(x_{\pm})=E$; that is to say,
$\dot{x}=0$ at those points, i.e. the kinetic energy becomes zero and E
converts to the potential energy.

It is well known [21], [28], that the period of the motion may be given by:

$T=\oint dt=\oint\frac{dt}{dx}dx,$ (6)

from which we will obtain the frequency $\Omega=2\pi/T$.

Since (from Eq. 5) $v^{2}=(2/m)(E-V(x))$, for a symmetric potential, we can
write that period is 4 times the integral between minimum and maximum of
reached potential between points $x_{min}$ and $x_{max}$.

$T=4\sqrt{\frac{m}{2E}}\int_{x_{min}}^{x_{max}}\frac{dx}{\sqrt{1-V(x)/E}}$ (7)

We are interested in the following potentials:

$V(x)=\epsilon_{0}\cdot\left|\frac{x}{\xi}\right|^{n},$ (8)

where $\xi$ is a certain length scale (actually, already determined in
previous section). The absolute value of displacement, $\left|x\right|$, is
there only for mathematical rigor. As we will see, in our case we will not
need to integrate in negative x-direction.

We can introduce a variable z (that will change integration from the range
[$x_{min}$:$x_{max}$] to [0:1]):

$z=V(x)/E=\frac{\epsilon_{0}}{E}\cdot\left(\frac{x}{\xi}\right)^{n},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
dx=\left(\frac{E}{\epsilon_{0}}\right)^{1/n}\cdot\frac{\xi dz}{nz^{(n-1)/n}}.$
(9)

Now, we can write Eq. 7 in the form:

$T=4\sqrt{\frac{m}{2E}}\cdot\left(\frac{E}{\epsilon_{0}}\right)^{1/n}\cdot\frac{\xi}{n}\int_{0}^{+1}\frac{dz}{z^{(n-1)/n}\sqrt{1-z}}.$
(10)

By introducing parameters p=1/n and q=1/2, we can write:

$T=4\sqrt{\frac{m}{2E}}\cdot\left(\frac{E}{\epsilon_{0}}\right)^{1/n}\cdot\frac{\xi}{n}\int_{0}^{+1}z^{p-1}(1-z)^{q-1}dz.$
(11)

We recognize that the integral expression in Eq. 11 is the definition of
Euler’s $\beta$ function [29]:

$\beta(p,q)=\int_{0}^{1}x^{p-1}(1-x)^{q-1}=\frac{\Gamma(p)\Gamma(q)}{\Gamma(p+q)},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ p>0,q>0.$ (12)

Hence, we have:

$T=4\xi\sqrt{\frac{m}{2\epsilon_{0}}}\cdot\left(\frac{E}{\epsilon_{0}}\right)^{\frac{2-n}{2n}}\cdot\Psi(n),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Psi(n)=\frac{1}{n}\cdot\frac{\Gamma(1/n)\Gamma(1/2)}{\Gamma(1/n+1/2)}.$ (13)

The introduced $\Psi(n)$ is a slowly decreasing function of $n$, with
$\Psi(1)=2$, $\Psi(2)=\pi/2$.

At the maximum of displacement (at $x=A$, amplitude of oscillations), based on
Eq. 8, potential energy equals the total energy $E$, and we get
$A/\xi=\left(E/\epsilon_{0}\right)^{1/n}$. Therefore, we can relate period $T$
to amplitude $A$:

$T=4\xi\sqrt{\frac{m}{2\epsilon_{0}}}\cdot(A/\xi)^{(2-n)/2}\cdot\Psi(n).$ (14)

Accordingly, the angular frequency can be written as:

$\omega=\frac{\pi}{2\xi}\cdot\sqrt{\frac{2\epsilon_{0}}{m}}\cdot\frac{1}{\Psi(n)}\cdot(A/\xi)^{(n-2)/2}.$
(15)

Hence, we see that while for $n<2$ period $T$ increases with $A$ sub-linearly,
for $n>2$ frequency increases with $A$ faster than linearly.

When $n>0$, there are three cases known to us when a close-form exact analytic
solutions are available: 1) A trivial case when $n=1$, resulting in a force
not dependent on distance between particles, describing motion of mass objects
on an inclined plane. 2) An ideal harmonic oscillator for $n=2$. 3) For $n=4$
a closed-form solution can be expressed with Jacobi elliptic functions [30],
[31].

The incomplete Beta function is when we integrate in Eq. 12 from $t=0$ to time
$t$ less than $T/4$. It is usually denoted as $\beta_{z}(p,q)$:

$\beta_{z}(p,q)=\int_{0}^{z}z^{p-1}(1-z)^{q-1}dz,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ 0\leq z\leq 1.$ (16)

Therefore, Eq. 16 gives us time as a function of $z$, i.e. $t=\beta_{z}(p,q)$.

###  3.2 Displacement method

Two basic methods may be used for studying the nonlinear properties of
anharmonic crystals. First, we would like to know well the dynamics of
interaction between two layers, only. For that we use the displacement method:
a structure as shown in Fig. 1 c) is used with two layers only, the upper and
the lower. The upper layer is moved in X-direction for a certain distance
$\Delta X$ (much less than the periodicity of crystal), away from the
equilibrium position, and at time $t=0$ it is released free while all atoms
are undergoing integration in NVE statistical ensemble. The upper layer is
treated as rigid body, while the lower layer as fixed. This is a realization
of the idea of a classical oscillator where a particle of mass m enters a
periodic in-time motion in a potential well. The quantities determined in
simulations are such as the position of all atoms, their velocities, forces
acting on them, etc. We perform averaging of those quantities for all atoms in
a layer.

Figure 2: Dependence of oscillation frequency $\omega$ on exponent of
potential, for displacements from 0.025 to 0.2 Å, as shown in the legend. The
fitting lines are functions based on Eq. 15,
$\omega(n)=a/\rm{\Psi(n)}\cdot(\Delta X/\xi)^{(n-2)/2}$, where parameters $a$
and $\xi$ are common for all curves.

Scaling of oscillation frequency on exponent of potential, and the initial
displacement $A$ (amplitude of oscillations) is shown in Fig. 2, for a few
values of initial displacements, by using Eq. 15 in the form: $\omega(n,\Delta
X)=a/\rm{\Psi(n)}\cdot(\Delta X/\xi)^{(n-2)/2}$. Hence, here
$a=\pi/2\xi\cdot\sqrt{2\epsilon_{0}/m}$, and $b=\xi$. The least squares
fitting of all the data gives values of $a$=42.35 ps-1 and $\xi$=2.020 Å. The
value of $\xi$ is in agreement with what is found from potential analysis
$E_{p}(x,n)$ in the previous section, where the method of creating potential
tables is described.

We use values $\epsilon_{0}$=8.4 eV and the above value of $\xi$ to compute
$a$: $a=(\pi/2\xi)\cdot\sqrt{2\epsilon_{0}/m}$, and we obtain 41.9 ps-1, very
close to the observed fitting parameter.

Figure 3: Displacement normalized by its initial value (amplitude, 0.2 Å in
this case) is shown as a function of time normalized by T. Additionally,
velocities $v_{x}$ as a function of normalized time (right vertical axis), are
shown. Maxima or minima of $v_{x}$(t) are found exactly at t=T/4 or at t=3T/4,
as indicated by vertical red lines. On the top of $v_{x}$(t) curves denoted in
the legend, we draw thin green lines $v_{x}$(t) computed from $\Delta$X(t)
data by using the formula: $v_{x}(t)=v_{max}(n)\cdot\sqrt{1-\left|\Delta
X(t)/A\right|^{n}}$ (taking care of proper sign, depending on time range).

Since $z$ in Eq. 16 is defined as $(\epsilon_{0}/E)\cdot(x/\xi)^{n}$, we can
find out the dependence $t(x)$. In particular, at $z\ll 1$ we can use an
approximation:
$\beta_{z}(p,q)\approxeq(1/p)z^{p}(1-z)^{q}\approxeq(1/p)z^{p}$. Therefore,
based on Eq. 11, we can write for $t\ll T/4$ (near the potential minimum):

$t(x)\approxeq\sqrt{\frac{m}{2E}}\cdot\left(\frac{E}{\epsilon_{0}}\right)^{1/n}\cdot\frac{\xi}{n}\cdot(n)z^{1/n}=\sqrt{\frac{m}{2E}}\cdot
x.$ (17)

Using that equation, we obtain velocity at potential minimum, which is the
maximal (in value), $v_{max}=\sqrt{2E/m}$. It does not depend on $n$ (as
desired): it depends on total energy $E$. Since
$E={\epsilon_{0}}\cdot(A/\xi)^{n}$, we have a scaling relation:

$v_{max}=\sqrt{2\epsilon_{0}/m}\cdot(A/\xi)^{n/2}.$ (18)

That dependence is valid for data as included in Fig. 3, and it can be
approximated (for amplitude $A$ ranging from 0.005 Å to 0.2 Å) as
$v_{max}(n)=54.4\cdot(A/\xi)^{n/2}$ Å/ps, with $\xi$=2.02 Å. Notice that the
proportionality coefficient, 54.4Å/ps differs from that in $\omega$ scaling in
Fig. 2 by factor 1.286 Å, which is exactly $2\xi/\pi$, from the scaling factor
in Eq. 15.

We obviously would like to have $x(t)$ dependence, that is an inverse
incomplete beta function, usually named as I(z,p,q) when it is regularized. It
is implemented in some available software packages. For our computational
purposes and drawing graphs a method for finding I(z,p,q) is to compute $t(z)$
as a two-column data file and display it by replacing columns, and next to
rise $z$ to power $1/n$. That is what we did in this research. We use Gnuplot
for computing incomplete beta function, defined as ibeta(p,q,z). However,
these dependencies coincide so well with the results of our simulations that
we feel no need for reproducing them. Instead, we will analyze certain other
relationships between various quantities.

Figure 4: Relationship between pressure (force) and displacement. Displacement
raised to power (n-1) reproduces pressure curves (with a certain rescaling
factor). The response has been recorded as an initial displacement of the
upper layer for 0.2 Å for a single period of oscillations, which is normalized
to 1 for all curves. Results are shown for 4 values of n: 1.5, 2.0, 3.0, and
3.6. Three kinds of curves are drawn: displacement of the upper layer $\Delta
X/A$ (left vertical axis), pressure at the surface $P_{xy}$ normalized by its
maximum value (right vertical axis), and displacement raised to power n-1,
$(\Delta X/A)^{n-1}$, normalized in such a way that the curves coincide with
pressure values. $(\Delta X/A)^{n-1}$ curves are drawn on the top of pressure
curves and largely cover them up. The colors in the legend refer to pressure
data.

Figure 3 shows a relationship between velocity and displacement. Displacement
is normalized by its amplitude (left vertical axis) while velocity (right
vertical axis) is not. Both quantities are shown as a function of $t/T$, for 4
values of $n$. Vertical red lines are located at $t=T/4$ and $t=3T/4$,
intersecting $\Delta X(t)$ at $\Delta X$=0 and $v_{x}(t)$ at its minimum and
maximum values. Velocity depends on kinetic energy, $v_{x}=\sqrt{2E_{k}/m}$,
while $E_{k}=E-V(x)$ and $V(x)=\epsilon_{0}(x/\xi)^{n}$. From that, the
relation follows: $v_{x}(t)=v_{max}(n)\cdot\sqrt{1-\left|\Delta
X(t)/A\right|^{n}}$. It is drawn with thin green lines on the top of $v_{x}$
data from MD simulations.

Similarly, we show a relationship between displacement and the surface
pressure $P_{xy}$. In LAMMPS pressure components are computed as a sum of
forces (averaged over the surface area) exerted on the surface from within the
entire sample volume. In our particular case, there are two only layers of
atoms. $P_{xy}$ is registered for the upper layer and it is caused by the
forces acting on it from the lower layer; see Fig. 1 c). This illustration
helps understand the physical significance of pressure in LAMMPS since often
it seems to be not understood properly. Hence, while
$V(x)=\epsilon_{0}(x/\xi)^{n}$, force is given by potential derivative,
$F=dV/dx=n\epsilon_{0}/\xi\cdot(x/\xi)^{n-1}$. Therefore, displacement raised
to a power (n-1) reproduces pressure curves, with a rescaling factor which is
related in this case to $n\epsilon_{0}/\xi$. In Fig. 4, the response has been
recorded to an initial displacement of the upper layer for 0.2 Å for a single
period of oscillations. Results are shown for 4 values of $n$: 1.5, 2.0, 3.0,
and 3.6. Three kinds of curves are drawn: displacement of the upper layer
$\Delta X$ normalized by its amplitude (0.2 Å in this case; left vertical
axis), pressure at the surface $P_{xy}$ normalized by its maximum value (right
vertical axis), and displacement to power $n-1$, $(\Delta X/A)^{n-1}$,
normalized in such a way that the curves coincide with pressure values.
Computed $(\Delta X/A)^{n-1}$ curves are drawn on the top of pressure curves,
by using tiny green lines that largely cover them up.

###  3.3 Heaviside method

The method described in the previous section, with two layers, where one is
displaced initially, is suitable for determining the fundamental oscillation
frequency related to the potential between layers. However, in typical MD
simulations, we rather investigate the response of material to an applied
pressure. The results in that case are, in general, different than those
derived to this point. They would be the same only in the case of a perfect
harmonic potential. If any nonlinearities are present, the effective potential
of interaction between layers after applying pressure becomes non-symmetric
concerning the minimum of the potential well (Fig. 5).

Figure 5: Potential energy curves without the force applied (upper 3 curves)
and with the force of F=1 applied (the bottom 3 ones), for n=1.5, 2,0, and
3.0, when $\epsilon_{0}=1$. The blue arrow indicates maximum displacement
(amplitude $A_{H}$), and red arrows show positions $x_{0}/A_{H}$, i.e. minima
of potential energy when force is applied. In this figure, amplitude $A_{H}$
is equal to 1 for all curves.

Hence, when pressure (force F) is applied in a Heaviside-type dependence on
time, to a sample as shown schematically in Fig. 1 c), the upper layer is
moved from the equilibrium position at $x=0$ and tends to achieve another
equilibrium position at a position $x_{0}$ in Fig. 5. Oscillations occur
around $x_{0}$. That equilibrium would be reached after a (possibly very long)
simulation times when losses were present. In the simplest model we are
studying we find no indication for the presence of any energy losses at
temperature $T=0$ and therefore the equilibrium $x_{0}$ is not supposed to be
reached ever. In simulations, it is easy however to introduce nonphysical
losses during oscillation that allow to determine quickly the positions of
equilibrium, $x_{0}$.

In the presence of a Heaviside type of force $F$, since force is a derivative
of potential, we may modify potential given by Eq. 8 to the following form:

$V(x)=\epsilon_{0}\cdot(x/\xi)^{n}-Fx.$ (19)

Let us assume that the total energy of the particle is zero at $x=0$ and
$x=A_{H}$ (see Fig. 5). We find the amplitude of oscillations $A_{H}$ from the
condition $V(A_{H})=0$ and $x_{0}$ from the condition $dV/dx=0$, resulting in:

$\frac{A_{H}}{\xi}=\left(\frac{\xi
F}{\epsilon_{0}}\right)^{\frac{1}{n-1}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \frac{x_{0}}{\xi}=\left(\frac{\xi
F}{n\epsilon_{0}}\right)^{\frac{1}{n-1}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{x_{0}}{A_{H}}=\left(\frac{1}{n}\right)^{\frac{1}{n-1}}.$ (20)

Notice that $x_{0}/A_{H}$ tends to 0 when $n\searrow 1$, and it tends to 1
when $n\nearrow\infty$. Also, the ratio $x_{0}/A_{H}$ does not depend on the
applied force $F$. It depends on $n$, only.

The depth of the potential well is $E_{0}=V(x_{0})$:

$E_{0}=\epsilon_{0}\cdot(x_{0}/\xi)^{n}-Fx_{0}=(1-n)\epsilon_{0}\cdot\left(\frac{\xi
F}{n\epsilon_{0}}\right)^{\frac{n}{n-1}}=(1-n)\epsilon_{0}\left(\frac{x_{0}}{\xi}\right)^{n}.$
(21)

We may assume now that the total energy of the particle is $|E_{0}|$. At the
bottom of the potential well (at $x=x_{0}$), $|E_{0}|$ converts into fully
kinetic energy, while at $x=0$ and $x=A_{H}$ the maxima of potential energy
are reached during particle oscillation.

This time we are not able to provide exact expressions on the period of
oscillations since we do not know how to perform integration analogous to that
one in Eq. 7 for potential given by Eq. 19. However, we expect that some
scaling properties derived so far will remain valid for the new potential.

The quantity that must not be valid in the case of the Heaviside method is the
period of oscillations and angular frequency since the potential has now a
different analytical form and it is asymmetric. Let us however have an
approximate expression on both. By substituting $|E_{0}|$ instead of $E$ in
eq. 13, we have a new period of oscillations:

$T=4\xi\sqrt{\frac{m}{2\epsilon_{0}}}\cdot(n-1)^{\frac{2-n}{2n}}\cdot\left(\frac{\xi
F}{n\epsilon_{0}}\right)^{\frac{2-n}{2(n-1)}}\cdot\Psi(n),$ (22)

and the angular frequency:

$\omega=\frac{\pi}{2\xi\cdot\Psi(n)}\cdot\sqrt{\frac{2\epsilon_{0}}{m}}\cdot(n-1)^{\frac{n-2}{2n}}\cdot\left(\frac{\xi
F}{n\epsilon_{0}}\right)^{\frac{n-2}{2(n-1)}}.$ (23)

Figure 6: Phase space portraits of the system for different situations. The
dependence of velocity on the displacement is shown. The upper two figures a)
and b) are for the displacement method of simulations, with amplitudes of
displacement as shown in the legend. Both, the displacements and velocities
are normalized by amplitude $A$. The c) and d) parts are when Heaviside
pressure is applied, for pressure values as shown in the legend (in GPa). The
displacements and velocities are normalized by the value of applied pressure
$P^{0}_{xy}$.

A good way of comparing the dynamics in the case of both methods
(displacement, and Heaviside) is by drawing phase space portraits, as in Fig.
6. This is also a convenient way of testing the validity of ongoing
simulations. The diagrams show dependencies of velocity versus displacement.
The upper two figures are for the displacement method, where we see phase
loops symmetric around both X and Y with the center at (0,0). For the
Heaviside method, the loops are symmetric only concerning Y axis.

Atoms movement must be restricted to X-direction only. Without these
restrictions, extremely small asymmetries in the sample structure or
fluctuation effects due to temperature will drive the motion into quasi-
chaotic oscillations (still preserving the attractor points). Such a quasi-
chaotic motion results in sub-harmonic components in any time dependencies,
that is in $\Delta X(t)$, $v_{x}(t)$, $P_{xy}(t)$, etc. That would manifest
itself in phase space diagrams as unclosed loops changing with time and
filling in densely a large portion of phase space.

Asymmetry of potential well manifests strongly in the shape of all time-
dependencies of the Heaviside method. Figure 7 shows how pressure reflects the
force acting on the upper layer that comes from the potential derivative.

Figure 7: The applied pressure is 500 MPa (that value is not significant in
this case). Displacement, the absolute value of velocity, and the surface
pressure are shown, when $n$=3.6. All quantities are normalized by the period
of oscillations $T$, and by their maximal values. Red vertical lines cross
velocity curves at their maximum values. The green horizontal line is drawn at
$\Delta X/A_{H}=x_{0}/A_{H}$. On top of pressure data the computed $(\Delta
X(t)/A_{H})^{n-1}$ is drawn with a tiny green line. .

The red vertical lines intersect velocity curves at their maximum values. The
green horizontal line is drawn at $\Delta X=x_{0}/A_{H}$ computed with Eq. 20.
The point of crossing of red and green lines labeled as $(t_{0},x_{0})$ is the
coordinate of the potential minimum. Velocity there has extreme value. $t_{0}$
is also the inflection point on displacement and pressure curves. We have no
equation on $t_{0}$, it can be determined only approximately by using for
instance this graphical method. When $n>2$, as in this case, $t_{0}/T$ is
larger than 1/4 and $x_{0}/A_{H}$ is larger than 1/2.

The pressure registered at the upper layer, $P_{xy}$, is the result of forces
acting on that layer from the lower layer and it can be written as force
(derivative of potential) divided by the surface area per atom, $S$:

$P_{xy}=-(n\epsilon_{0}/\xi)\cdot(x/\xi)^{n-1}\cdot(1/S).$ (24)

The maximum pressure at the surface is when $x=A_{H}$ is reached, and it is
given by:

$P_{xy,max}=-(n\epsilon_{0}/\xi S)\cdot(A_{H}/\xi)^{n-1}.$ (25)

Therefore, we can write immediately:

$\frac{P_{xy}(t)}{P_{xy,max}}=\left(\frac{\Delta X(t)}{A_{H}}\right)^{n-1}.$
(26)

Hence, normalized displacement raised to a power $(n-1)$ reproduces pressure
curves normalized by the maximum pressure value, as is shown by tiny green
lines drawn over the data points of pressure in Fig. 7. At the point
$(t_{0},x_{0})$ velocity reaches maximum (or minimum) value.

Figure 8: This is similar to Fig. 4 but for simulations under Heaviside type
of constant pressure of 1000 MPa, and similar to Fig. 7. Displacement vertical
axis is on the left. For pressure and $\Delta X^{n-1}$ curves the right
vertical scale is used. Notice the difference in the pressure axis: now it
spans the range from -1 to 0 while in Fig. 4 it is in the range -1 to 1. Green
horizontal lines mark the positions of $x_{0}$ (for n values as shown by green
labels), which is the point of equilibrium around which oscillations occur,
i.e. that is the position of the minimum of potential energy. Vertical red
lines mark the time positions $t_{0}$ when $x_{0}$ is passed through. The red
lines cross with green lines at the same points as the displacement curves
pass the green lines. At these time moments, velocity has maximum (or minimum)
values (not shown, to avoid the complexity of the picture).

Similar results are shown in Fig. 8 for two additional values of $n$.

While we know values of $x_{0}$ from Eq. 20, we do not know values of $t_{0}$,
since we have no exact solutions of equations of motion. The pressure curves
(with points) coincide perfectly with rescaled $(\Delta X/A_{H})^{n-1}$ tiny
curves in Fig. 8.

We can draw some additional conclusions regarding the pressure curves.

At the point $(x_{0},t_{0})$ we have $P_{xy}(x_{0})=(n\epsilon_{0}/\xi
S)\cdot(x_{0}/\xi)^{n-1}$. Therefore
$P_{xy}(x_{0})/P_{xy,max}=\left(x_{0}/A_{H}\right)^{n-1}$. Since
$x_{0}/A_{H}=(1/n)^{1/(n-1)}$ (Eq. 20), we can write:

$\frac{P_{xy}(x_{0})}{P_{xy,max}}=\frac{1}{n}.$ (27)

That relation can be verified by the data in Fig. 8. Note only that green
lines are drawn with the left vertical axis, while pressure is described by
the right axis. Therefore, one should use the vertical red lines indicating
$t_{0}$ to determine the crossing point of $x_{0}$ by pressure curves.

Moreover, $x_{0}$ is the equilibrium position that would have been reached
when losses of energy existed. Hence, pressure at that point, $P_{xy}(x_{0})$
is compensated by the applied pressure $P^{0}_{xy}(x_{0})$. Therefore, the
ratio of the maximal value of pressure oscillations $P_{xy,max}$ to the value
of the surface pressure applied $P^{0}_{xy}$, equals $n$, as we observe indeed
in simulations:

$P_{xy}(x_{0})=P^{0}_{xy},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \frac{P_{xy,max}}{P^{0}_{xy}}=n.$ (28)

Surface area per atom is $S=a^{2}\sqrt{3}/4=5.49\cdot 10^{-20}m^{2}$, where
$a$ is lattice constant, $a$=3.56 Å, The proportionality coefficient between
pressure and displacement in Eq. 24, $\epsilon_{0}/(S\xi)$, is 121.2 GPa when
computed with $\epsilon_{0}$=8.4 eV. This value is consistent with what we
find from the analysis of dependence between $P_{xy,max}$, and $A_{H}$ for
various pressures and $n$’s.

Figure 9 shows the scaling of oscillation frequency and amplitude for
simulations performed under constant pressure, for a broad range of pressure
values.

Based on Eq. 23, we may write:
$\omega(n,P)=a/\Psi(n)\cdot(n-1)^{\frac{n-2}{2n}}\cdot(bP/n)^{\frac{n-2}{2(n-1)}}$,
where P is pressure in GPa. Here the parameter a has the same meaning as in
the case of Fig. 2, a=$\pi/2\xi\cdot\sqrt{\epsilon_{0}/2m}$. By using least
squares fitting to all of the data, we obtain $a$=42.25 ps-1, while in the
case of Fig. 2 that value was 42.35 ps-1. The difference is insignificant and
it is within the range of fitting error. For $b$ we find $b$=0.0195 GPa-1.

$A_{H}(n,P)$ in Fig. 9 is fitted to all the data, also with two only
parameters that are the same for all curves. Based on Eq. 20 we may write:
$A_{H}(n)=\xi\cdot\left(bP\right)^{\frac{1}{n-1}}$, where $b$ has the same
meaning as in expression on $\omega(n,P)$: $b=\xi S/\epsilon_{0}$.

The lines for $A_{H}(n,P)$ are drawn with $\xi$=1.97 Å and $b$=0.0083 GPa-1.
Here, $\xi$ is very close to the value found already for a symmetric
potential, of about 2.020 Å. The value of $b$ is about 2.35 times smaller than
the one from fitting $\omega(n,P)$. However, this value is consistent with
what is obtained from computing $b=\xi S/\epsilon_{0}=0.00824GP^{-1}$, with
$\xi$=2.02 Å and $\epsilon_{0}$ = 8.4 eV. It agrees also with the inverse
parameter, 121.2 GPa, determined already from $P(A_{H})$.

While Eq. 20 on amplitude scaling must be correct, the one on frequency ought
to be treated as an approximation only, as mentioned already when deriving
Eqs. 22 and 23.

Figure 9: Scaling of oscillation frequency $\omega$ with exponent of potential
when simulations are performed under constant pressure, for pressure values as
shown in the legend. The left vertical log scale is for $\omega(n)$ (full
circles) and the right linear scale is for amplitude $A_{H}(n)$ (squares). The
colors of symbols and lines are the same for $\omega(n)$ and for $A_{H}(n)$.

In summary, well-established scaling relations are valid in the case of the
Heaviside method. If we change the force $F$ (pressure $P^{0}_{xy}$) by a
factor $\lambda$, $F\rightarrow\lambda F$, new potential depth can written as
$E_{0,\lambda}=\lambda^{\frac{n}{n-1}}\cdot E_{0}$. When simultaneous change
is made to scales of distance and time, equations of motion remain the same
for new quantities (Eqs. 20, 21, 22):

$\begin{split}x_{0,\lambda}=\lambda^{\frac{1}{n-1}}\cdot
x_{0},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
A_{\lambda}=\lambda^{\frac{1}{n-1}}\cdot A,\\\
T_{\lambda}=\lambda^{\frac{2-n}{2(n-1)}}\cdot T,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\omega_{\lambda}=\lambda^{\frac{n-2}{2(n-1)}}\cdot\omega.\\\ \end{split}$ (29)

For instance, either by analyzing the ratio of $x$ to $t$ or by using the
energy rescaling we find that velocities are transformed in this way:
$v\rightarrow\lambda^{\frac{n}{2(n-1)}}\cdot v$. Acceleration will scale up as
a square of velocity: $a\rightarrow\lambda^{\frac{n}{n-1}}\cdot a$. Since it
is related to force (pressure) acting on the particle, therefore the same
scaling will be valid for pressure as well:
$P_{xy}\rightarrow\lambda^{\frac{n}{n-1}}\cdot P_{xy}$

##  4 Nonlinear crystals

###  4.1 Dependencies on the number of layers N.

The description of the two-layer motion of crystallographic layers is the
introduction only to the study of dynamics of many-layer structures and will
serve us as theoretical background for finding out how stress (sound)
propagates in crystals. The model that we developed is essentially a
1-dimensional model of a chain of masses connected by springs. That model has
been found to describe very well stress penetration as analyzed by MD
simulations in steel and it was very helpful in improving the accuracy of data
analysis from simulations, as well as in developing a methodology for the
simulation itself.

The first ideas about the usefulness of this model were proposed in [2]. There
some evidence has been presented also on the existence of pressure
oscillations within the sample interior. Since that time the model itself has
been greatly improved, as well the techniques of conducting MD simulations and
data analysis have been excelled. Details of this model with supporting data
analysis will be published separately [3].

So far we introduced the period of oscillations (and angular frequency) for
two interacting layers. However, in the case of multi-layer material, we have
yet another period of oscillation and angular frequency. Once the pressure is
applied to the upper surface of the sample, as shown in Fig. 1 c), it causes
oscillations of the upper layer. That one interacts with the second one, etc.
In other words, a pressure wave is induced at the surface and that wave
propagates into the material interior with the speed of sound (in our case it
is the transverse sound mode). Once the wave reaches the lower layer of the
sample, it becomes reflected there (depending on boundary conditions imposed
on that lower layer). The reflected wave interferes with the incoming wave,
and reaches the sample surface, where it is reflected. The process repeats; In
cases when no losses are present, a perfect periodic dependence on time is
found for quantities such as displacement, velocity of layers, pressure at the
surface, etc, at least when we have an ideal harmonic interactions between
layers. That period of oscillations depends however on the size of sample in
direction of wave propagation (Y-direction in our case). As a measure of
sample size we use the number of layers N and the period related to sound wave
traversing the sample will be denoted as $T_{N}$.

Figure 10: Displacement of the upper layer normalized by its maximum value as
a function of time normalized by period of oscillations $T_{N}$, for n=1.5 and
n=2.3, at 1000 MPa, for a sample with N=33 layers. Nearly triangle-like
dependencies of $\Delta X(t)$ are observed for any values of $n$. The inset
shows expanded displacement (not normalized by $\Delta X_{max}$) as a function
of time (not normalized by $T_{N}$) for a large number of curves for different
$n$, for very short times. Red dots there have been computed for $n$=2 by
using Eq. 30.

We will show some results obtained for samples of changing sizes, from N=2 to
N=20. Most simulations are obtained for a sample with N=33.

Typical displacement versus time for one period of oscillations, where both
quantities are normalized, $\Delta X(t)$ by its maximum value $\Delta X_{max}$
and time by period of oscillations $T_{N}$, is shown in Fig. 10. Pressure at
the surface, $P_{xy}(t)$, follows closely $\Delta X(t)$ dependence.

The inset shows expanded displacement as a function of time (not normalized by
$T_{N}$) for very short times. In that short time region, $\Delta X(t)\sim
t^{2}$ is found for any $n$. Moreover, $\Delta X$ scales up linearly with the
value of applied pressure, in the same time range. For $n$=1.1 this
approximation works well below around 0.002 ps. For most curves shown it is
valid below around 0.01 ps or more, and for all curves it is valid for $t/T$
less than around 0.1, where $T$ is the oscillation period for the 2-layer
system. Hence, in that region layer dynamics may be viewed as it was caused by
constant force. This is the same as observed in Figs. 8 and 7 for $t\ll T$.
The red dots in the inset of Fig. 10 have been calculated for $n$=2, by using
Eq. 30, as will be explained.

Figure 11: The velocity of the uppermost layer at 200 MPa pressure applied,
for a sample with N=33 layers, when time is rescaled by the period of
oscillation, for 2 first oscillation periods.  Figure 12: Velocity $v_{x}$
determined from simulations under pressure 1000 MPa is shown for times
$t<T_{N}/2$ for three values of n: 1.5, 2.0, and 2.5, for figures from top to
bottom. There are 9 curves drawn for different numbers of layers, from N=2 to
N=20, as shown with green labels in the top figure. The curves coincide
together for initial times, regardless of the value of N.

The velocity of the uppermost layer is shown in Fig. 11 for two periods of
oscillations. We see that while the amplitude of velocity changes strongly
with $n$, the time dependence $v_{x}(t)$ remains approximately the same for
all curves when drawn with time normalized by $T_{N}$.

The one full period of oscillations takes the time needed for 4 times of
traversing the wave through sample length in Y direction [3]. At $t=T_{N}/4$
and $t=3T_{N}/4$ we usually notice a slight change in curves, when reflection
of wave occurs from the lower layer, especially when $n$ is larger than 2 and
applied pressures are high.

The points when $v_{x}(t)$ passes through value of 0 can be used for an
accurate determination of the half of period of oscillations, $T_{N}/2$ (Fig.
11). We have verified that $T_{N}$ determined in this way depends linearly on
N, $T_{N}=\alpha\cdot(N-1)$, where $\alpha$ depends on the parameter of
potential nonlinearity $n$, i.e. the sound wave propagates similarly as in
case of harmonic potential, when $n$=2. Period of oscillations (and speed of
sound) depend on applied pressure, though.

In Fig. 12 velocity of the uppermost layer is shown for the first half of the
oscillation period for samples with different numbers of layers, as labeled in
part a), ranging from N=2 to N=20. The curves for larger N cover up the data
of smaller samples. We compare the results for 3 values of n. The horizontal
axis range in every case of different $n$ has been chosen in a way that the
vertical red lines cross $v_{x}$=0 axis at $T_{N}/2$ for all values of $n$.
This illustrates that the period of oscillations $T_{N}$ scales up in the same
way with the number N of layers in samples, regardless of nonlinearity
parameter $n$. Let us notice that for N=2 the curves $v_{x}(t)$ look quite the
same, regardless of $n$. This is an inspiration for a more accurate analysis
of $v_{x}(t)$ at the shortest times, and that is the clue towards
understanding of stress propagation into the bulk of material.

###  4.2 Pater’s model of a chain of masses and springs for n=2.0

The analytical model of a chain of masses and springs that describes correctly
our results for $n$=2, is provided in a little-known, old article by Pater
[32].

Pater’s exact analytical results reproduce properly the expected dynamics in
the limit of continuous medium. In that case, for instance, the wave entering
an elastic material exposed to surface Heaviside-type pressure, continues to
penetrate the bulk of the material, preserving its original shape. In the case
however of a discrete medium, as the situation of our interests, the solutions
are given in terms of Bessel functions. Let in our case $m$ enumerate
particles in the chain, starting from number 1 as that one exposed to the
initial Heaviside force. In our situation, crystallographic layers play the
role of these particles. The quantities like displacements $u_{m}$, velocities
$v_{m}$, and force $f_{m}$, acting on a particle $m$, after their proper
amplitude’s normalization, are given by:

Figure 13: Pressure is computed with Pater equations, and it is compared with
MD simulation results of $P_{xy}(t)$ under 1 GPa Heaviside pressure impulse,
applied to a sample with N=16 layers and n=2.00 (the thick green curve in the
background). The period of pressure oscillations in this case is 2.547 ps,
therefore at t=0.63675 ps (which is the scale range of the horizontal axis)
the front of the pressure wave is supposed to reach the opposite side of the
sample, penetrating the whole sample volume. The labels indicate the range of
layers over which the summation of forces in analytical equations is
performed.

$u_{m}=J_{2m}(2\theta)+\sum_{k=2m+2,2m+4,...}^{\infty}(k-2m+1)J_{k}(2\theta),$
(30)

$v_{m}=J_{2m-1}(2\theta)+2\cdot\sum_{k=2m+1,2m+3,...}^{\infty}J_{k}(2\theta),$
(31)

$f_{m}=-\left[1-J_{0}(2\theta)-2\cdot\sum_{k=2,4,...}^{2m}J_{k}(2\theta)+J_{2m}(2\theta)\right],$
(32)

where $J_{m}$ are Bessel functions of the first kind and $\theta=\Omega t$.
The $\Omega$ there, for n=2, is $\sqrt{2\epsilon_{0}/m}$. These equations are
valid for $t<T_{N}/4$, i.e. before the front of the wave reaches the opposite
side of the sample, where reflection occurs. The equations still could be used
for longer times, however, we want to avoid complicating this description.

In MD simulations, pressure at the surface $P_{xy}$ is a volume average of the
internal virial stress, $\langle S_{xy,m}\rangle$, or it can treated as a
summation of forces exerted on the surface from all layers inside the sample
volume. Some authors treat virial stress as a not well-defined quantity, with
no clear physical meaning. We have however checked accurately and verified
that $S_{xy,m}$ is a very useful quantity and it fits well into the physical
description of material properties.

Hence, we identify $f_{m}$ as given by Eq. 32 with virial stress averaged over
all atoms of layer $m$, $\langle S_{xy,n}\rangle$. Therefore we can write:

$P_{xy}(t)=\sum_{k=1,2,...}^{m}f_{k}(t),$ (33)

where the maximum summation range $m$ should not exceed the number of layers N
and the summation is valid for times $t<T_{N}/4$. Computational results of
$P_{xy}(t)$ by using this method are shown in Fig. 13, where comparison with
MD simulations are included (MD results are drawn in the background as a thick
green line). We use successive approximations of $P_{xy}(t)$ based on Eg. 33,
by summing a finite number of contributions from forces $f_{m}(t)$, where $m$
denotes the number of layer (counting from the surface of the sample, where
force is applied). The labels in Fig. 13 indicate the range of layers over
which summation is performed. At very short times, a few summation terms are
sufficient for approximation. Once time becomes longer, more terms must be
included, while the pressure wave penetrates deeper and more layers are
involved in the interaction.

Figure 14: Comparison of analytical results obtained with Bessel functions
(green lines in the background) with MD simulations (dots drawn on top of
these) for velocity (a) and virial stress (b) of layers. Labels denote layers
number, from the uppermost 33 (In MD simulations layers are counted from the
sample bottom). The surface Heaviside impulse pressure applied is 2000 MPa,
the sample is with N=33 layers, and $n$=2.00. Here $\Omega$=24.82/ps, which is
the same as in Fig. 13. The virial pressure (force) is computed with Pater’s
Eq. 32, and the velocity of layers with Eq. 31. The red horizontal line in b)
is drawn at 0.1 of applied surface pressure and points where it is crossed by
$S_{xy}(t)$ curves may be used for determining the speed of sound.

###  4.3 Virial stress propagation.

Figure 15: Comparison of dynamics of stress penetration into the sample
interior for three values of nonlinearity parameter n: 1.5, 2.0, and 3.5, when
the surface applied pressure is 1000 MPa. The upper layer is at around Y=68 Å.
Green lines are drawn for $t=T_{N}/4$. For n=2, this movie illustrates better
the process of stress penetration.

At very short times, less than the time needed for the sound to pass the
distance between two layers, the first only contribution, $f_{1}$ of Eq. 33 is
sufficient for the approximation of surface pressure. In that time range MD
results on $\Delta X(t)$ have a parabolic $t$-dependence, in perfect agreement
with Pater analytical results. That is shown by red dots in the inset of Fig.
10, computed for $n$=2 in the same way (by using Eq. 30 for displacement).

The velocity of layers is computed with Eq. 31 and compared with MD simulation
results in Fig. 14 a), for a sample with N=33, for $n$=2, under a surface
pressure of 2000 MPa. First, by green lines, the results of computation with
Bessel functions are drawn, for layers from the uppermost, labeled as layer
33. On top of these lines, with points, the results of simulations are shown.
There are two only parameters needed for obtaining a coincidence between
analytical and MD simulation curves. These are $\Omega$ in Eq. 31 that
determines the time scale, and the amplitude of velocity, both in agreement
with all the scaling relationships discussed so far.

For the first, uppermost layer velocity can be written as an infinite sum (Eq.
31): $v_{1}=J_{1}(2\theta)+2\cdot J_{3}(2\theta)+2\cdot J_{5}(2\theta)+...$.
However, for a sample with the number of layers equal to N (N $\geq$ 2)
summation must be limited to $J_{k}$, where k=4(N-1)-1. Hence, the first layer
oscillations are shown by green lines for samples with N from 2 to 20. A
proper adjustment of time and amplitude is done so that a good coincidence
with the simulation data is obtained.

While in Fig. 13 we demonstrate that Pater’s equations allow us to reconstruct
surface pressure as a function of time, in Fig. 14 b) we show that virial
stress for layers inside of sample volume can also be computed with these
equations, with a very good agreement between analytical results and MD
simulations data. The labels indicate the layers, as in Figure a). There are
two important details here. 1) In LAMMPS, the entire force is reported as
acted on atoms, while Pater equations provide force acting on one side of the
layer, only. Therefore, to compare the results of LAMMPS simulations and
Pater’s, we need to take a sum of two results from Pater for neighboring
planes. 2) There is an exception from the rule 1). In the case of the
uppermost layer, the force in LAMMPS simulations acts on the layer from one
side only (since there is no layer above it). The same is in the case of the
lowermost layer. Therefore in these cases the Patter result for one layer only
ought to be used. This explains why the virial stress for the layer labeled as
33 (the uppermost) in Fig. 14 b) is about 2 times smaller than that for the
remaining layers.

Calculations like these, with the use of Eqs. 30-32, have however a limited
applicability and can be used for not more than a few tens of layers.
Computing Bessel functions requires the summation of a large number of terms,
which finally leads to computational instabilities when either the number of
layers and/or time becomes too large. In our case, for computing Bessel
functions, we used Cephes Perl libraries222Math::Cephes module in Perl..

While Fig. 14 b) shows the virial stress as a function of time for different
layers, in Fig. 15 we draw virial stress $S_{xy}$ profiles across the sample
length in the Y-direction, for increasing time moments, with the time
separation of 0.001 ps. The curves of $S_{xy}(y,t)$ have been obtained by
applying an approximate cubic splines method to a discrete data set of points
(one point for one layer, with layer’s separation of 2.06 Å). Pressure is
applied at the surface of the sample (with a Y-coordinate of around 68 Å), and
the stress front penetrates the interior, reaches the opposite sample side at
Y=0, where it is reflected, and the reflected wave interferes with the
incoming wave in a constructive way causing an additional increase of virial
stress. The light-green lines show the stress profiles at approximately
$t=T_{N}/4$.

Results shown here suggest that the amplitude of oscillations in quantities
like $v_{x}(t)$, $S_{xy}(t)$ or in $S_{xy}$ as a function of the y-coordinate,
grows up with time or with spacial coordinate. This question can not be
resolved by performing analytical modeling by using Bessel functions as
provided by Pater - the range where accurate approximation of these functions
is available is too narrow. We did preliminary MD simulations performed on
very long samples (up to 1000 layers; 0.2 $\mu$m in size) but still were not
able to find clear conclusions. The nature of oscillations changes at long
times (distances) suggesting rather that a sort of another wavefront builds up
on top of the initial oscillations while the amplitude of the initial ones
decreases. The problem is related to the harmonic case as well when n=2, but
the phenomena are stronger for n increasing well over 2 and for increasing
pressures at the surface. Some of our preliminary studies suggest however that
under high pressures a catastrophic, chaotic scenario of stress propagation is
realized. On another hand, there are theoretical arguments that amplitude of
similar oscillations decreases with time/distance [7]. This subject requires
separate systematic research. Perhaps the new experimental techniques
available [1] could provide also physical evidence for the nature of these
phenomena.

###  4.4 Sound velocity

It is difficult to provide a precise, rigorous criterion for determining the
time when stress waves entering the volume start to increase at a position of
a given layer. That increase in time is initially very slow. The red
horizontal line in Fig. 14 b) is drawn at 0.1 of applied surface pressure and
the points of crossing it by $S_{xy}(t)$ curves have been used for determining
the speed of sound (this is a transverse sound component, with propagation in
direction [111]). We could use as a criterion some another value of stress,
e.g. half of the applied pressure. That would lead to an insignificant
difference in the determined value of sound speed.

Full circles in Fig. 16 a) show sound velocity determined in that way, i.e.
from the criterion of reaching 0.1 of $P^{0}_{xy}$, when pressure applied is
1000 MPa. In the case of linear theory (when $n$=2), it follows from
dispersion relations for propagation of waves in a periodic medium that sound
velocity $c$ is related to angular frequency $\Omega$ (the same which is used
in Pater’s equations and the same as determined by us for oscillations of
2-layers), and to the distance $d$ between layers, $c=d\cdot\Omega$. In our
case $d=a/\sqrt{3}$, where lattice constant $a$ is 3.56 Å. For comparison,
empty circles in Fig. 16 a) show the speed of sound computed with that
equation. It does reproduce perfectly the speed determined for n=2, but for
other values of nonlinearity parameter deviations are significant.

We think that the velocity of stress propagation is influenced largely by the
asymmetry of the potential well, as shown in Fig. 5. The time needed for a
particle to move on the left of the potential minimum is different than the
time for traversing the right part: for $n$>2 the time on the right is shorter
than the time on the left. Hence, we argue that the time on the right
determines the speed of sound. This is well seen in Fig. 16 b), where velocity
normalized by the period of oscillations $T$ for a 2-layer case is shown for a
range of $n$-values.

The point of minimum is marked as $(x_{0},t_{0})$ in Fig. 7. Position $x_{0}$
is known from Eq. 20. Since we have no analytical expression on $x(t)$, we can
not provide an exact formula on $t_{0}$. One of the possible methods of
determining $t_{0}/T$ is to find it from the position of maximum in velocity,
as in Fig. 16 b). The results are shown in Fig. 16 c), together with
$x_{0}/A_{H}$.

When we correct the sound speed computed from equation $c=d\cdot\Omega$ by a
factor $1/(2-4t_{0}/T)$ (this is like assuming that the time $T/2-t_{0}$
should replace $T/4$), the results become very close to these determined from
the measurements of the front wave, as shown by empty squares and the blue
curve (for guiding eyes, only), in Fig. 16 a).

Figure 16:  a) Speed of sound determined as a ratio of the distance between
stress curves for layers 32nd and 1st, like these in Fig. 14 b) and the time
needed for the arrival of pressure wave to layer 1, is shown with full
circles. b) Change with $n$ of the position of velocity maximum in a 2-layer
case has been used to find out the ratio of $t_{0}/T$, as shown in c). The
speed of sound corrected by the factor $1/(2-4t_{0}/T)$ is shown in a) with
empty squares.

##  5 Summary and conclusions

The method has been described of creating a table-style nonlinear inter-atomic
potentials between NN atoms in FCC crystals for use in MD simulations in
LAMMPS. Exact analytical formulas were found on the period of oscillations for
a system of two-layer crystallographic planes, when one of them is initially
displaced, and approximate equations are provided for the case when Heaviside-
type pressure is applied. Scaling relations are given between quantities such
as average displacement, velocities, and forces acting on the planes. The
results of that analysis helped understand the dynamics of stress propagation
into the bulk of crystal consisting of many (N) crystallographic planes.

We used analytical Pater’s model of the chain of masses and springs to find
out that it describes perfectly well the data obtained by MD simulations in
the case of harmonic inter-atomic potential. However, we found also that in
nonlinear cases Pater’s model remains qualitatively valid as well, since a
quasi-linear relation between the wave vector and frequency of oscillations
remains preserved. That allowed us to determine that the speed of sound is a
strongly decreasing function of $n$. In case however, when nonlinearity
parameter $n$ becomes larger than 2 and/or in case of large applied pressures,
oscillations may become unstable and chaotic. The subject of stability of
propagating waves of stress requires separate systematic research. The jumps
in the velocity of dislocations (observed by us as well, in certain
conditions) may be caused by incoming front waves of increased stress. These
waves are present in the case of the harmonic medium as well, and in the case
when dislocations are present their strength may only be stronger and more
chaotic.

Oscillations like those of $v_{x}(t)$, $S_{xy}(t)$, etc, can easily be
misinterpreted as a result of typical time dependencies found when pressure is
applied abruptly (Heaviside-type) and energy losses are present. Oscillations
in such cases can be observed in $P_{xy}(t)$ as well (some of the results on
$P_{xy}(t)$ presented in [2] have been misinterpreted by us).

Similar oscillations of various quantities are often observed in MD
simulations, either due to some justified physical origin or caused by an
unintentional mistake in the methodology of performing the simulations.
Therefore, the significance of Pater’s model is in explaining, that in a case
of perfect harmonic inter-atomic potential and no energy losses, we still
ought to observe oscillations in certain conditions that are related to the
discrete nature of atoms and their periodic arrangement in crystals.

An unresolved question is how adding more degrees of freedom (in particular,
allowing the movement of atoms in Z-direction) will change the dynamics of
stress propagation? What will be the effect of chaotic movement due to the
finite temperature on the propagation of waves and stability of crystal
lattice? Will finite temperature change the effective inter-atomic potential?
Would it make sense to consider fluctuation-dissipation theory as a
description of mechanisms introducing energy losses and leading to the damping
of oscillations?

## References

  * [1] T. S. Holstad, L. E. Dresselhaus-Marais, T. M. Ræder, B. Kozioziemski, T. van Driel, M. Seaberg, E. Folsom, J. H. Eggert, E. B. Knudsen, M. M. Nielsen, H. Simons, K. Haldrup, H. F. Poulsen, Real-time imaging of acoustic waves in bulk materials with x-ray microscopy, PNAS 120(39) (2023) 1–4.  doi:10.1073/pnas.2307049120.
  * [2] Z. Kozioł, Size effects in stress propagation and dynamics of dislocations: Fe–ni–cr steel, Modelling Simul. Mater. Sci. Eng. 30 (2022) 065010.  doi:10.1088/1361-651X/ac83de.
  * [3] Z. Kozioł, The role of stress propagation in dynamics of dislocations in crystals, to be published (2023).
  * [4] C. D., S. R., F. R., Theories of complexity, Complexity 8(3) (2003) 19–30.  doi:10.1002/cplx.10059.
  * [5] M. Herrera, T. M. Antonsen, E. Ott, S. Fishman, Dynamic localization of a weakly interacting bose-einstein condensate in an anharmonic potential, Phys. Rev. A 87 (2013) 041603(R).  doi:10.1103/PhysRevA.87.041603.
  * [6] K. N. Pathak, Theory of anharmonic crystals, Phys. Rev. 139(5A) (1965) A1569.  doi:10.1103/PhysRev.139.A1569.
  * [7] S. C. Lowell, Wave propagation in monatomic lattices with anharmonic potential, Proc. R. Soc. Lond. A 318 (1970) 93–106.  doi:10.1098/rspa.1970.0135.
  * [8] M. Hutereau, P. A. Banks, B. Slater, J. A. Zeitler, A. D. Bond, M. T. Ruggiero, Resolving anharmonic lattice dynamics in molecular crystals with x-ray diffraction and terahertz spectroscopy, Phys. Rev. Lett. 125 (2020) 103001.  doi:10.1103/PhysRevLett.125.103001.
  * [9] C. Nucera, F. L. di Scalea, Nonlinear wave propagation in constrained solids subjected to thermal loads, Journal of Sound and Vibration 333 (2014) 541–554.  doi:10.1016/j.jsv.2013.09.018.
  * [10] S.-Q. Hu, D.-Q. Chen, S.-J. Zhang, X.-B. Liu, S. Meng, Probing precise interatomic potentials by nonadiabatic nonlinear phonons, Materials Today Physics 27 (2022) 100790.  doi:10.1016/j.mtphys.2022.100790.
  * [11] A. von Hoegen, R. Mankowsky, M. Fechner, M. Forst, A. Cavalleri, Probing the interatomic potential of solids with strong-field nonlinear phononics, Nature 555 (2018) 79.  doi:10.1038/nature25484.
  * [12] G. Bonny, N. Castin, D. Terentyev, Interatomic potential for studying ageing under irradiation in stainless steels: the FeNiCr model alloy, Modelling and Simulation in Materials Science and Engineering 21 (8) (2013) 085004.  doi:10.1088/0965-0393/21/8/085004.
  * [13] L. Béland, A. Tamm, S. Mu, G. Samolyuk, Y. Osetsky, A. Aabloo, M. Klintenberg, A. Caro, R. Stoller, Accurate classical short-range forces for the study of collision cascades in Fe-Ni-Cr, Comput. Phys. Comm. 219 (2017) 11–19.  doi:10.1016/j.cpc.2017.05.001.
  * [14] M. Grado-Caffaro, M. Grado-Caffaro, Phonon energy in an anharmonic quasi-one-dimensional solid, Journal of Physics and Chemistry of Solids 72 (2011) 957–960.  doi:10.1016/j.jpcs.2011.05.002.
  * [15] A. Carré, S. Ispas, J. Horbach, W. Kob, Developing empirical potentials from ab initio simulations: The case of amorphous silica, Computational Materials Science 124 (2016) 323–334.  doi:10.1016/j.commatsci.2016.07.041.
  * [16] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comp. Phys. 117 (1995) 1.  doi:10.1006/jcph.1995.1039.
  * [17] O. Hellman, I. A. Abrikosov, S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84 (2011) 180301(R).  doi:10.1103/PhysRevB.84.180301.
  * [18] S. R. E., T. A., B. L. K., S. G. D., S. G. M., C. A., S. L. V., O. Y. N., A. A., K. M., W. Y., Impact of short-range forces on defect production from high-energy collisions, J. Chem. Theory Comput. 12(6) (2016) 2871 – 2879.  doi:10.1021/acs.jctc.5b01194.
  * [19] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, New York, 2008.
  * [20] N. Euler, Transformation properties of x” + f1(t)x’ + f2(t)x + f3(t)xn = 0, Journal of Nonlinear Mathematical Physics 4(3-4) (1997) 310–337.  doi:10.2991/jnmp.1997.4.3-4.7.
  * [21] P. Amore, F. M. Fernández, Exact and approximate expressions for the period of anharmonic oscillators, Eur. J. Phys. 26(4) (2005) 589.  doi:10.1088/0143-0807/26/4/004.
  * [22] T. Harko, F. S. N. Lobo, M. K. Mak, Integrability cases for the anharmonic oscillator equation, Journal of Pure and Applied Mathematics: Advances and Applications 10(1) (2013) 115–129.  arXiv:1304.1468.
  * [23] M. Znojil, I. Semorádová, Log-anharmonic oscillator and its large-n solution, Modern Physics Letters A 33(38) (2018) 1850223.  doi:10.1142/S0217732318502231.
  * [24] R. Barreto, M. F. Carusela, A. G. Monastra, Thermal conductance of suspended nanoribbons: interplay between strain and interatomic potential nonlinearity, J. Stat. Mech. (2017) 103201doi:10.1088/1742-5468/aa8c18.
  * [25] A. Pankov, Traveling waves in Fermi-Pasta-Ulam chains with nonlocal interaction, Discrete and Continuous Dynamical Systems - S 12(7) (2019) 2097–2113.  doi:10.3934/dcdss.2019135.
  * [26] S. Kashchenko, The interaction of waves in the Fermi–Pasta–Ulam model, Communications in Nonlinear Science and Numerical Simulation 91 (2020) 105436\.  doi:10.1016/j.cnsns.2020.105436.
  * [27] E. Fermi, P. Pasta, S. Ulam, M. Tsingou, Studies of the nonlinear problems (5 1955).  doi:10.2172/4376203.
  * [28] L. D. Landau, E. M. Lifshitz, Mechanics: Volume 1 (Course of Theoretical Physics), 3rd Edition, Butterworth-Heinennan, Oxford, 1976.
  * [29] M. Abramowitz, I. Stegun (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th Edition, New York: Dover, 1970\.
  * [30] D. F. Lawden, Elliptic Functions and Applications. Applied Mathematical Sciences, Vol. 80, Springer-Verlag, New York, 1989.
  * [31] E. T. Whittaker, A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, 4th Edition, Cambridge University Press, Cambridge, 1964.
  * [32] A. D. de Pater, The propagation of waves in a semi-infinite chain of material points and springs representing a long train, Vehicle System Dynamics: International Journal of Vehicle Mechanics and Mobility 3:3 (1974) 123–140.  doi:10.1080/00423117408968452.