Dynamics of stress propagation in anharmonic crystals: MD simulations

Anharmonic inter-atomic potential x to power of n, n>1, has been used in molecular dynamics (MD) simulations of stress dynamics of FCC oriented crystal. The model of the chain of masses and springs is found as a convenient and accurate description of simulation results, with masses representing the crystallographic planes. The dynamics of oscillations of two planes is found analytically to be given by Euler beta functions, and its scaling with non-linearity parameter and amplitude of oscillations, or applied shear pressure is discussed on examples of time dependencies of displacements, velocities, and forces acting on masses (planes). The dynamics of stress penetration from the surface of the sample with multiply-planes (an anharmonic crystal) towards its interior is confirmed to be given exactly as a series of Bessel functions, when n=2 (Schrodinger and Pater solutions). When n<>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.

Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

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 insight into sub-picosecond time scales and subatomic spatial resolution.
There are 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 extracting the linear response from simulation data.On another hand, how to find out the value of pressure acting on a dislocation in material, while the penetration process is a dynamic one of a comparable timescale as the time of simulation itself?Additionally, during the pressure front passage through the dislocation, 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 a related 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 provide a mathematical description of pressure penetration.At the same time, we excelled in our technical skills, and 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 the simulated physical phenomena, is of the crucial importance.
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 [3][4][5][6]: (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 knowing the dynamic and thermodynamic properties of real materials.
Non-linearity in the physics of crystals is omnipresent.Even small, non-linear contributions to inter-atomic potentials, may considerably influence the dynamics of particles, leading to chaotic motion and may cause seemingly unexpected energy losses [7].It may cause localization of vibrations, depending on the strength of the non-linear contribution [8].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 [9].Lattice dynamics and conditions of its stability have been addressed as well [10].
The presence of anharmonic potentials in materials is confirmed now by a broad spectrum of experimental techniques.X-ray diffraction and terahertz time-domain spectroscopy is used to resolve potential in molecular crystals [11], or higher harmonics detection in ultrasound propagation in steel [12] is used, or analysis of phonon oscillation spectra in graphene and other materials [13].Higher-order phonon non-linearities are observed by optical, infrared, and terahertz spectroscopy [14].
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 et al [15] and Béland et al [16] EAM potentials commonly used in simulations of steel have a small anharmonic contribution.In the case of Bonny's, an anharmonic addition may be approximated by |x| 2.5 .
Hence, in this article, we study the simplest possible non-linear inter-atomic potentials of the form V(x) = ϵ 0 x n , where n > 1, mostly it is between n = 1.5 and n = 3.60.Since closedform analytical solutions to equations of motion in that case are not known in general, we perform simulations of the dynamics of particles by using an MD simulation tool, LAMMPS [19].It is an overlooked opportunity that MD simulations may be used as a tool for inspiration and for guiding our understanding of physical processes, with the simulation results largely playing a role like that of experimental data.MD is often used in studies of lattice dynamics of solids [20].
In section 2 we describe the simulation setup and creation of non-linear 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, and pressure.
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 of the chain of masses and springs available for harmonic crystals.The model was first solved mathematically by Schrödinger in 1914 [21,22], and independently, by using another mathematical approach, by de Pater in 1974 [23].Both these fundamental works remain wildly unknown.Our simulation results provide the first evidence (to the best of our knowledge) for the validity of their solutions.The approach based on the model of masses connected by spring was found useful in case of explaining, for instance, Raman spectra in a few layers graphene [24,25], and in our analysis of Van der Waals interaction in graphene [26].
Finally, we discuss how the speed of sound is influenced by the non-linearity parameter n.

Simulations setup, potentials, and samples
Since we want to have the simulation conditions similar as much as possible to the situation we have already studied intensively and very accurately, we use a special configuration suitable when the dynamics of dislocations in FCC steel is studied.We assume the FCC structure of the material, which is oriented in a specific way [2]: Y direction is along [111] crystallographic axis, X is in [1 10] and Z in [ 11 2].In dynamical studies either an initial displacement of a layer of atoms is made in X direction, or pressure is applied to the upper sample layer in that direction, as shown by the label 'Force' and arrow in figure 1(c).Figures 1(a  In this orientation of the sample, it follows from symmetries that a possible 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 top (upper region in figure 1) to N for the lowermost layer (lower region).The number of atoms in every plane must be the same, and the uppermost region must contain one only layer, and also exactly the same number of atoms as all other layers.Otherwise, the simulation results obtained will be noisy and hard to interpretation.In particular, for instance, when the upper rigid region has a different than 1 number of layers, the time dependence of pressure on the sample surface as reported by LAMMPS will be entirely different, and it will depend also on the number of layers in that region.This fact is commonly miss-understood and ignored.
To achieve these conditions, the sample 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 Ycoordinate.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 the number of layers, N, ranges from 2 to 33.The size of samples in X-and Z-direction does not play a role in our simulations (at least in cases when it is larger than a few unit cells and as long as the number of atoms in a single layer allows to achieve good statistical averages over physical quantities).Simulations performed in this work have been conducted at temperature T = 0 (that allows us to avoid sample annealing).
We use the table pair style in LAMMPS to provide interatomic potential.A text file is created for that by using Perl scripts.It contains the 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.It means we are limited to some maximal values of the applied displacements or pressure.The cut-off on the lower side plays less role in simulations itself.We used 0.5 Å for that (such a short cut-off is not needed, as we found out later).For the convenience of numerical computation, we use the following formula for the tablestyle potential used in LAMMPS: 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.That potential is 0 at r − r 0 = 1 [Å] and it has value −ϵ 0 at r = r 0 , for any n.
We aim to reproduce possibly closely the results for n = 2 to these observed in steel in our previous work.Therefore, we choose ϵ 0 = 4.2 eV at the minimum of the potential well, as in the case of the potential of Stoller et al [27] characterized at the NIST 1 repository.For particles, we assign the mass of one atom of Fe.
The arrangement of atoms is shown in figures 1(a) and (b) (all atoms are in the equivalent surroundings), where crystal orientation in figure (a) is the same as in (c).
Let us concentrate on finding out the potential energy of the green atom.
If to assume that the position of the green atom is 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: When the green atom is moved in the X-direction for x, the distance between it and 3 blue atoms becomes: where r 0 = a/ √ 2. Let us write: r i /r 0 = √ 1 + y, where i is 1, 2 or 3. We have: Hence, to the lowest order of x: Since in (1) E p is a function of the absolute value of (r − r 0 ), by adding contributions from 3 atoms, we obtain as an approximation of E p (x), with respect to the potential minimum: 1 NIST, Interatomic Potentials Repository.Potential energy change Ep(x) of a single atom in a layer displaced for x with respect to the other layer, for several values of the non-linearity parameter n, from n = 1.5 to n = 3.5, as shown by labels.On the top of sets of data points for any given value of n, tiny green lines are drawn by using Ep(x) = 2ϵ 0 • (x/ξ) n , with one pair of parameters ϵ 0 = 4.240 eV, and ξ = 2.006 Å, common for all values of n.The inset illustrates the method of computation of Ep in LAMMPS: dots symbolize positions of atoms.There are three groups of atoms created: green, red and blue; each of them belonging to their own group.The blue group of atoms is moved in the X-direction for finding change in potential energy, which is computed for red atoms only.
Let us notice that (6) does not depend on r 0 in that lowest order approximation.
The potential has been checked in two ways.First, we calculated numerically E p (x) by using equations ( 1)-( 4), and we 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.Both methods give the same results.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.15%), with the overall potential described by this kind of function (figure 2): where ϵ 0 = 4.240 eV and ξ = 2.006 Å.When n decreases below 1.5, the accuracy of fitting the data by using equation (7) deteriorates quickly.Hence, assuming inter-atomic potential between two NN atoms in the form of (1) leads to the potential of interaction between any atom with all its NN on the crystallographic lattice in the form given by ( 6) as an analytical approximation, or as (7), as an approximation of MD simulation results.
The inset in figure 2 illustrates the method of computation of E p in LAMMPS, as we used it: dots symbolize positions of atoms.There are three groups of atoms created: green, red and blue; each of them belonging to their own group.The blue group of atoms is moved in the X-direction for finding change in potential energy, which is computed for red atoms only.
Attention must be paid to the proper determination of potential energy from MD simulations.Every atom has 12 NN, 6 of them on the same plane, and 3 on every of the two neighboring planes.Therefore, the total potential energy of an atom ought to be −12ϵ 0 .However, the potential energy per one atom, when x = 0, as reported by LAMMPS, is exactly −25.2 eV, i.e. it is −12ϵ 0 /2 (where ϵ 0 = 4.2 eV).When computing total potential energy of a system with a large number of atoms, and averaging it over all atoms, we obtain the same value, i.e. −12ϵ 0 /2.However, when we determine the potential energy change of a group of atoms in one layer, after it is moved (let say for red atoms in inset of figure 2), we ought to take into account that the total potential energy change of the entire system contains also contribution from the potential energy of green atoms as well, since these atoms are moved as well relative to red atoms.Change of their potential energy is exactly the same in value.It is easier to understand that by using the following example: let us assume that we have a system of two only atoms separated by distance r 0 , i.e. both are at the minimum of their potential energy well.In that case their total potential energy is −ϵ 0 .However, LAMMPS will report their potential energy per atom as −ϵ 0 /2, a value which is two times lower.Hence, we conclude that the potential energy of a single atom that is used in any later calculations of their dynamics must be two times larger than the reported by LAMMPS potential energy per atom.Not taking that into account results in a disagreement by a factor √ 2 between the observed in simulations and the computed frequency of oscillation of atoms (layers).The factor √ 2 is due to the scaling of oscillation frequency as square root of potential energy.
As a result of imposed conditions, the average X for the region lower does not change with time.Average Y for regions lower and upper remain the same.

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.
It is 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.[28]).
Finding exact, general solutions for potentials of that form, and approximate ones, have been discussed by Euler [29], Amore and Fernández [30], Harko et al [31].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 [31].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 → 0, [32].
Variants of the chain of masses and springs model have been used often to study aspects of non-linearities, chaos, and localization/de-localization of lattice vibrations [12,[33][34][35].The model gained a broad interest after the Fermi-Pasta-Ulam paradox was published in 1955 [36].
There is another way of gaining an understanding of solutions to the problem, and it is based on the use of molecular dynamics simulations.The method is flexible, and allows one to approach problems more and more accurately by introducing or removing some physical interactions and/or constraints, approximating them gradually to the real.
Before going into the analysis of the results of MD simulations, let us concentrate on presenting the known analytic solution to the problem of simple anharmonic oscillator with the power-law dependence of potential energy on the displacement between two particles.
We are interested in the following potentials: where ξ is a certain length scale (actually, already determined in the previous section).The absolute value of displacement, |x|, is there only for mathematical rigor.As we will see, in our case we will not need to integrate in a negative x-direction.
From the equation of motion, mẍ = −dV(x)/dx, we obtain an integral of the motion E, which is the total energy: It is well known [30,37] that the time of the motion may be given as an integral over the trajectory of motion C in the x − v phase space: In our case, the periodic motion of the particle is restricted to the interval x − < x < x + , where the turning points x ± satisfy V(x ± ) = E; that is to say, ẋ = 0 at those points, i.e. the kinetic energy becomes zero and E converts to the potential energy.Moreover, for any potential V(x) that is symmetric with respect to its minimum value at x = 0, the motion during one period consists of traversing the closed trajectory, for instance: x + → 0, 0 → x − , x − → 0, and 0 → x + .Equations of motion for any of these 4 trajectories are the same.Therefore, the period is 4 times the time of the movement of the particle between x = 0 and x = x + .Hence, since from (9) v 2 = (2/m)(E − V(x)), after substituting v(x) into (10), we have: We can introduce a dimensionless variable z representing the ratio of the potential energy to the total energy: Now, at the turning point x + potential energy V(x + ) equals to E, hence z = 1 in ( 12) at that point.Therefore the integration range in (11) changes from [0:x + ] to [0:1] when z is used instead of x as an integration variable.We can write (11) in the form: By introducing parameters p = 1/n and q = 1/2, we have: We recognize that the integral expression in equation ( 14) is the definition of Euler's β function [38]: Hence, we have: The introduced At the maximum displacement (at x = A, amplitude of oscillations), based on ( 8), the potential energy equals the total energy E, and we get A/ξ = (E/ϵ 0 ) 1/n .Therefore, we can relate the period of oscillations T to the amplitude A: Accordingly, the angular frequency can be given as: When n > 0, there are three cases known to us when close-form exact analytic solutions are available: (1) A trivial case when n = 1, resulting in a force not dependent on the distance between particles, describing the 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 [39,40].

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 figure 1(c) is used with only two layers, the upper and the lower.The upper layer is moved in X-direction for a certain distance ∆X (much less than the periodicity of the 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 a rigid body, while the lower layer is fixed.This is a realization of the idea of a classical oscillator where a particle of mass m enters a periodic intime motion in a potential well.The quantities determined in simulations are the coordinates of all atoms, their velocities, potential energy, components of virial stress tensor.We perform averaging of those quantities for all atoms in every layer.It ought to be noticed that the default resolution used in LAMMPS for reporting quantities such as potential energy, positions, and velocities of atoms ought to be increased for 2-3 orders to achieve the accuracy of simulations as reported by us.(18), , where parameters a and b are common for all curves.

Frequency of oscillations.
The period of oscillations T may be determined from the dependence on time of the surface pressure, P xy (t), from the displacement of the upper layer ∆X(t), or its velocity v x (t).Usually, we use the least-squares fitting method to find out the period by approximating the simulation data with the function sin(2π t/T).Albeit the fitted data in general, deviate strongly from that function, this is an accurate method for finding the oscillation period.
The scaling of the oscillation frequency on the exponent of potential is shown in figure 3, for a few values of the initial displacement, by using (18) in the form: Here, a = π/2ξ • 2ϵ 0 /m, and b = ξ.Notice that the parameter ϵ 0 in equation ( 18) must be multiplied by 2 when computing frequency in agreement with the potential experienced during MD simulations, as expressed by equation (7).The least squares fitting of all the data gives values of a = 42.35ps −1 and ξ = 2.020 Å.The value of ξ is in agreement with what is found from potential analysis E p (x, n) in the previous section.When ϵ 0 = 2 × 4.24 eV and ξ = 2.006 Å are used, a = (π/2ξ) • 2ϵ 0 /m = 42.4 ps −1 is obtained.

Fourier analysis.
Displacement curves as a function of time step are shown in figure 4.These data were prepared to find out how the components of a Fourier series describing X(t) evolve when n changes.We used fast Fourier transform (FFT) method.It requires a sequence of data points that are equally spaced in time and their number must be a power of 2, therefore we have 2 12 = 4196 points for one period of oscillations.
In the case of the displacement method, we performed FFT analysis of two kinds of data sets of X(t) sequences.First, we constructed X(t) dependencies numerically.Let us explain how that can be done.
When integration in ( 14) is done to z < 1, the resulting function is called the incomplete beta function: It allows us to find time as a function of z, t(z), (where z(x) is given by ( 12)).We obviously would like to have the reverse relation, x(t).For that, the inverse incomplete beta function is used, usually named as I(z, p, q) when it is regularized.We compute it numerically, for instance with betaincinv from Python module scipy.special.One ought to be aware that this function computes the 1/4th of the trajectory (similarly as β z in equation (19) gives us t(x) for 1/4th of T).Therefore, it is necessary to join properly together results of four similar calculations for obtaining the full oscillation period, as shown in figure 4 for a few values of n.
The second type of data sets analyzed by FFT are those obtained from MD simulations.MD results are practically identical to those obtained by the numerical method.However, obtaining reliable and accurate simulation data requires attention to certain details.First of all, we need to find out an exact value of the oscillation period T. In the next step, we ought to adjust the time step in LAMMPS script so that it equals exactly T/4196, in our case.
For FFT computations, we used Perl module Math::FFT.Instead, one can use for instance Python libraries through scipy.fft.
In our case, it is actually not necessary to use FFT for finding Fourier series amplitudes.An alternative method, giving a comparable accuracy, is to use fit command available in Gnuplot.One can approximate X(t) dependence by a finite sum of harmonics of sin(t) and/or cos(t), let say in the form: The fit command can find out coefficients a i for k equal up to a few tens.
Since functions shown in figure 4 are asymmetric with respect to t = 0, the Fourier decomposition must contain only sin(t) type of components.Moreover, due to the property X(t) = −X(T − t), only odd components will be present.Figure 5 shows the dependencies of amplitudes a k on the non-linearity parameter n when the overall functions are defined as values of all a k are positive, while for n > 2, there are values of both signs.Our data suggest that some a k (n) change signs more than once (as may be seen in this scale for a 5 (n) only).In the limit when n = 1, X(t) can be constructed from parabolic dependencies.It is easy to verify by computing proper Fourier integrals that in this case, all odd coefficients are given by a series a k = 32/π 2 • 1/k 3 (see also [41] and [42]).Red arrows on the left in figure 5 show these values for k = 1, 3, 5.In the limit of n → ∞, ∆X(t)/A tends to a triangle function of time.In that case, we found by integration that the Fourier series amplitudes are given as The arrows on the right show these limiting values.
It ought to be noticed that Fourier series amplitudes a k do not depend on the amplitude of oscillation A, when we normalize them by A, and time is normalized by T.

Relationships between displacement, velocity and pressure.
Since z in equation ( 19) is defined as (ϵ 0 /E) • (x/ξ) n , we can find out the dependence t(x).In particular, at z ≪ 1 we can use an approximation: β z (p, q) ≊ (1/p)z p (1 − z) q ≊ (1/p)z p .Therefore, based on equation ( 14), we can write for t ≪ T/4 (near the potential minimum): Using that equation, we obtain velocity at potential minimum, which is the maximal (in value), v max = 2E/m.It does not depend on n (as desired): it depends on the total energy E. Since E = ϵ 0 • (A/ξ) n , we have a scaling relation: That dependence is valid for data as included in figure 6, and it can be approximated (for amplitude A ranging from 0.005 Å to 0.2 Å) as v max (n) = 54.4 • (A/ξ) n/2 Å ps −1 , with ξ = 2.02 Å.Notice that the proportionality coefficient, 54.4Å ps −1 differs from that in ω scaling in figure 3 by a factor of 1.286 Å, which is exactly 2ξ/π, from the scaling factor in (18).
Figure 6 shows a relationship between velocity and displacement.Displacement is normalized by its amplitude (left vertical axis), and 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 ∆X(t) at ∆X = 0 and v x (t) at its minimum and maximum values.Velocity depends on kinetic energy, From that, the relation follows: 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 figure 1(c).This illustration helps us understand the physical significance of pressure in LAMMPS.Hence, while V(x) = ϵ 0 (x/ξ) n , force is given by the derivative of potential, F = dV/dx = nϵ 0 /ξ • (x/ξ) n−1 .Therefore, displacement raised to a power (n − 1) reproduces pressure curves, with a rescaling factor that is related in this case to nϵ 0 /ξ.In figure 7, 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 ∆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, (∆X/A) n−1 , normalized in such a way that the curves coincide with pressure values.Computed (∆X/A) n−1 curves are drawn on the top of pressure curves, by using tiny green lines that largely cover them up.

Heaviside method 3.3.1. Potential curves and frequency of oscillations.
The method described in the previous section, with two layers, where one of them is initially displaced, is suitable for determining the fundamental oscillation frequency related to the potential between layers, as we did for graphene [26] and steel [2].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 from those derived to this point.They would be the same only in the case of a perfect harmonic potential.If any non-linearities are present, the effective potential of interaction between layers after applying pressure becomes non-symmetric concerning the minimum of the potential well (8).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 ∆X/A (left vertical axis), pressure at the surface Pxy normalized by its maximum value (right vertical axis), and displacement raised to power n − 1, (∆X/A) n−1 , normalized in such a way that the curves coincide with pressure values.(∆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.
Hence, when pressure (force F) is applied in a Heaviside-type dependence on time, to a sample, as shown schematically in figure 1(c), the upper layer is moved away from the equilibrium position at x = 0 and tends to achieve another equilibrium position at a position x 0 in figure 8. Oscillations occur around x 0 .That equilibrium would be reached after a (possibly very long) simulation time 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 non-physical losses during oscillation that allow us to determine quickly the positions of equilibrium, x 02 .In the presence of a Heaviside type of force F, since force is a derivative of potential, we may modify potential given by (8) to the following form: Let us assume that the total energy of the particle is zero at x = 0 and x = A H (see figure 8).We find the amplitude of oscillations A H from the condition V(A H ) = 0, and x 0 is determined from the condition dV/dx = 0, resulting in: Figure 8. Potential energy curves without the force applied (the upper 4 curves) and with the force of F = 1 applied (the bottom 4 ones), for n = 1.2, 2,0, 3.0, and 10.0, when ϵ 0 = 1 and ξ = 1, according equation (22).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.
Notice that x 0 /A H tends to 0 when n ↘ 0, and it tends (very slowly) to 1 when n ↗ ∞, with x 0 /A H = 1/e ≈ 0.3679 for n = 1.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 ): We may assume now that the total energy of the particle is At the bottom of the potential well (at x = x 0 ), E converts into fully kinetic energy, while at x = 0 and x = A H the maxima of potential energy are reached during particle oscillation.
Let us notice that for n → ∞, V(x) is dominated by the Fx term when x ≲ x 0 , i.e. the potential there is well approximated by a linear one, as illustrated in figure 8 by the curve for n = 10.This time we are not able to provide exact expressions for the period of oscillations since we do not know how to perform integration analogous to that one in (11) for potential given by equation (22).However, we expect that some scaling properties derived so far will remain valid for the new potential.Hence, let us have approximate expressions for ω and T. By substituting |E 0 | instead of E in ( 16), we have a new period of oscillations: and the angular frequency: . (26)

Phase space portraits.
A good way of comparing the dynamics in the case of both methods (displacement, and Heaviside) is by drawing phase space portraits, as in figure 9.This is also a convenient way of testing the validity of ongoing simulations.The diagrams show the 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.In MD simulations, atoms movement for regions upper and lower must be restricted to X-direction only.Without these restrictions, minimal 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 ∆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.

Fourier series analysis.
Preparing quality MD simulation data for analysis of Fourier harmonics by using Heaviside method requires care in the proper choosing of the applied pressure value.At small values of n, close to 1, the amplitude of oscillations is very small and a large pressure must be applied and short timestep used.At large values of n, the situation is reversed.The displacement curves shown in figure 10 have been obtained for pressures in the range 10 −4 − 10 10 Pa and for timestep parameter between 10 −5 and 10 10 ps.The dependencies ∆X(t) are however not sensitive to the value of applied pressure when time is normalized by T and the displacement by A H . Let us pay attention that the normalization by A H differs from the normalization by the amplitude A used in the case of the displacement method.Due to the symmetries of X(t) we conclude right away that only harmonics of cos(t) will contribute to the displacement function.Moreover, since X(t) = X(T − t), even harmonics must be present.
When n → ∞, X(t) tends to a parabolic function, as shown by the two crossing blue curves in figure 10.The reason for that is that the potential well for x < x 0 may be approximated by a linear dependence on the displacement, as illustrated by the curve for n = 10 in figure 8, and x 0 → A H at the same time.This allows us to determine what are the limiting values of amplitudes of harmonics for large values of n.We performed proper calculation of Fourier integrals in that case, obtaining that ∆X(t)/A H can be represented by the Fourier series of coefficients a k = −(8/π 2 ) • (−1) k /k 2 (see also [41] and [42]).Arrows on the right in figure 11 are drawn at levels corresponding to values of these coefficients.
In the case of n ↘ 1, when n becomes close to 1, MD simulation results become unreliable, for reasons that are not fully understood by us, and therefore we show results for n starting from 1.2, and higher.
As it is seen in figure 11, for n < 2 all Fourier components are of the same sign, while for n > 2 the odd and even harmonics have the opposite sign.

Relationships between the displacement, velocity and pressure.
Asymmetry of potential well manifests strongly in the shape of all time-dependencies of the Heaviside method.Figure 12 shows how pressure reflects the force acting on the upper layer that comes from the potential derivative.
The red vertical lines intersect velocity curves at their maximum values.The green horizontal line is drawn at ∆X = x 0 /A H computed with (23).The point of crossing of red and green lines, labelled as (t 0 , x 0 ), is the coordinate of the potential minimum.Velocity there has an extreme value.t 0 is also the inflection point on displacement and pressure curves.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.
values of all a k are positive, while for n > 2 the odd and even harmonics have the opposite sign.When n → ∞, X(t) tends to a parabolic function, as shown by the two crossing blue curves.The calculated Fourier integrals in that limiting case give us: Arrows on the right are drawn at levels corresponding to values of these coefficients.The lines are for guiding eyes, only.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 x 0 /A H . On the top of pressure data, the computed (∆X(t)/A H ) n−1 is drawn with a tiny green line.
The pressure registered at the upper layer, P xy , is the result of forces acting on that layer from the lower layer.Instead of using the mass of the entire layer and its surface, it is convenient in any calculations to use the average mass of a single atom and its averaged surface area per atom, S. Hence, P xy related to the force (a derivative of potential), may be written as: The maximum pressure at the surface is when x = A H is reached, and it is given by: with what is obtained from computing b = ξ S/ϵ 0 = 0.0082GP −1 , with ξ = 2.02 Å and ϵ 0 = 8.48 eV.It agrees also with the inverse parameter, 122.4 GPa, determined already from P(A H ).
While equation ( 23) on amplitude scaling must be correct, the one on frequency ought to be treated as an approximation only, as mentioned already when deriving equations ( 25) and (26).
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 λ, F → λF, the new potential depth can be written as E 0,λ = λ n n−1 • E 0 .When simultaneous change is made to scales of distance and time, equations of motion remain the same for new quantities (equations ( 23)-( 25)): 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 λ → λ n 2(n−1) • v. Acceleration will scale up as a square of velocity: a λ → λ n n−1 • a.Since it is related to the force (pressure) acting on the particle, the same scaling will be valid for pressure as well: P λ,xy → λ n n−1 • P xy .

Anharmonic crystals
The description of the two-layer motion of crystallographic layers is the introduction only to the study of the dynamics of many-layer structures and will serve us as theoretical background for finding out how stress (sound) propagates in crystals.The model we use is essentially a 1-dimensional model of a chain of masses connected by springs, with masses representing crystallographic layers.That model has been found to describe well the stress penetration as analyzed by MD simulations in steel, it helped improve the accuracy of data analysis from simulations, as well it was helpful in developing a methodology for the simulation itself.
The first ideas about the usefulness of this model were proposed in [2], where also some evidence can be found for the existence of pressure oscillations within the sample interior.Details of this model and supporting data analysis will be published separately.

Dependencies on the number of layers N
So far, we introduced the period of oscillations (and angular frequency) for two interacting layers.However, in the case of multi-layer materials, 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 figure 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 the boundary conditions imposed on that lower layer).The reflected wave interferes with the incoming wave and reaches the upper sample surface.The process repeats.In cases when no losses are present, a perfect periodic dependence on time is found for quantities such as the displacement of layers, their velocity, pressure at the surface, etc.That period of oscillations depends, however, on the size of the sample in the 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 waves traversing the sample will be denoted as T N .
We will show some results obtained for samples of changing sizes, from N = 2 to N = 20, with most simulations obtained for a sample with N = 33.
Typical displacement versus time for one period of oscillations, where both quantities are normalized, ∆X(t) by its maximum value ∆X max and time by the period of oscillations T N , is shown in figure 14.Pressure at the surface, P xy (t), follows closely ∆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, ∆X(t) ∼ t 2 is found for any n.Moreover, ∆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 being caused by constant force.This is the same as observed in figure 12 for t ≪ T. The red dots in the inset of figure 14 have been calculated for n = 2, by using (33), as will be explained.The velocity of the uppermost layer is shown in figure 15 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 by the wave the sample length in Y direction.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 the applied pressure is high.
The points when v x (t) passes through a value of 0 can be used for an accurate determination of the half of the period of oscillations, T N /2 (15).We have verified that T N determined in this way depends linearly on N, T N = α • (N − 1), where α depends on the parameter of potential non-linearity 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 the applied pressure, though, when n ̸ = 2.
In figure 16 velocity of the uppermost layer is shown for the first half of the oscillation period for samples with different number 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 range of the horizontal axis, in every case of different n, has been chosen in such 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 the non-linearity parameter n.Let us notice that for N = 2 the curves v x (t) look quite the same, regardless of n.For N > 2 the curves are nearly identical as well (the difference is in their amplitude).This suggests that the observed oscillations and their similarity is not related to the non-linearity parameter n.

Pater's model of the chain of masses and springs for n = 2.0
The analytical model of a chain of masses and springs that correctly describes our results for n = 2, is provided in a little-known, old article by de Pater [23].We will use equations from that paper in the analysis of our simulation results.However, the model of a chain of masses and springs dates back as far as the beginning of the XX century, to these times when the differences between continuous and discrete mechanics were formulated, with the first exact analytical results attributed to Schrödinger [21,22].The ingenious work of Schrödinger remained however not known as well.
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, is given by: where J k are Bessel functions of the first kind, and θ = Ωt.The Ω there, for n = 2, is 2ϵ 0 /m, F 1 is the force applied at the first layer, the surface one, and m is the mass of the layer.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 could still be used for longer times, however, we want to avoid complicating this description.
In MD simulations, pressure at the surface (as reported by LAMMPS), P xy , is a volume average of the internal virial stress, ⟨S xy ⟩/V, where V in this definition is an average volume occupied by an atom.Virial stress can be treated as a summation of forces exerted on the surface by all layers inside the sample volume.Some authors consider virial stress as a not well-defined quantity, with no clear physical meaning.It is, however, a statistical quantity, like temperature, valid as well for non-equilibrium systems, as long as these can be assumed as ergodic.The literature contains a reach discussion of the subject [43][44][45][46][47].
We have checked accurately and verified that S xy may be used consistently in MD simulations.In particular, for instance, the relation P xy = ⟨S xy ⟩/V provides a valuable test of the validity of simulations: any discontinuities or change of slope in that dependence strongly indicates problems with the simulations.Finding the proportionality coefficient between these quantities is the best way to determine the average volume occupied by one atom 3 .
Hence, we identify f m as given by (35) to be related to the virial stress averaged over all atoms of layer m, ⟨S xy,m ⟩.Therefore, we can write: where S is the surface area of the uppermost layer, and the summation is valid for times t < T N /4, i.e. for times less than the time needed for the front of the pressure wave to reach the opposite surface of the sample.The maximal allowed range of index m in the above summation is N − 1, since there are a maximum of N − 1 layers contributing to the pressure at the surface layer 4 .In a case when the equilibrium of forces (stresses) was reached in the sample volume, the surface pressure would equilibrate N − 1 contributions from the pressure of N − 1 interior layers.That is why we ought to include a weight 1/(N − 1) in the above summation.Therefore, the proper form of (36) for computing the pressure contribution of layers ought to be: where we substituted equation (35) as an expression on forces from kth layer and used there P 0 xy instead of F 1 /S, where P 0 xy is the amplitude of the applied surface Heaviside pressure.Computational results of P xy (t) by using this method are shown in figure 17, where comparison with MD simulations are included, for samples with the number of layers from N = 2 to N = 20.That figure illustrates two aspects of computing surface pressure by using (37).First, with green lines we illustrate how the slope of P xy (t) changes with N. It is governed by the 1/(N-1) factor in (37), only (neglecting the obvious dependence on the amplitude of applied pressure, P 0 xy , and the timescale of the horizontal axis, which is determined solely by Ω).The second aspect of calculations is illustrated by magenta lines, for the sample with N = 20 layers.We use successive approximations of P xy (t) based on (36), by summing a finite number of contributions from forces f m (t), where m denotes the number of layers (counting from the surface of the sample, where force is applied).The magenta labels in figure 17 indicate the range of layers over which the summation is performed.
The only free parameter that needs adjustment for the proper description of P xy (t) is Ω.We used for it the same value in numerical calculations as used elsewhere, i.e.Ω = 24.82ps −1 .

Virial stress propagation
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 equation ( 36) is sufficient for the approximation of surface pressure.In that time range, MD results on ∆X(t) have a parabolic t-dependence, in perfect agreement with the Pater analytical results.That is shown by the red dots in the inset of figure 14, computed for n = 2 in the same way (by using (33) for the displacement).
The velocity of layers is computed with (34) and compared with the MD simulation results in figure 18(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 1.On the top of these lines, with points, the results of simulations are shown.There are only two parameters needed for obtaining a coincidence between analytical and MD simulation curves.These are Ω in (34) 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 (34): While in figure 17 we demonstrate that Pater's equations allow us to reconstruct surface pressure as a function of time, in figure 18(b) we show that virial stress for layers inside of Here Ω is 24.82 ps −1 , the same as in figure 17.The virial pressure (force) is computed with Pater's equation (35), and the velocity of layers with equation (34).The red horizontal line in (b) is drawn at 0.1 of the applied surface pressure.
the sample volume can also be computed with these equations, with a very good agreement between analytical results and MD simulation 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 layers, 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 Pater result for one layer only ought to be used.This explains why the virial stress for the layer labeled as 1 (the uppermost) in figure 18(b) is about 2 times smaller than that for the remaining layers.
The amplitude of velocity used in figure 18(a) that best fits the MD simulation results is 0.505 Å ps −1 .That amplitude is given by F 1 /mΩ (34).With the mass of the Fe atom 9.273 × 10 −26 kg, surface area per atom 5.49 Å 2 , and Ω = 24.82ps −1 , the amplitude of velocity is found to be 0.48 Å ps −1 .
Calculations like these, with the use of equations ( 33)- (35), 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 Math::Cephes Perl libraries.
While figure 18(b) shows the virial stress as a function of time for different layers, in figure 19 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 (available in Gnuplot) to a discrete data set of points (one point for one layer, with a 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 and 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.
The results shown here suggest that the amplitude of oscillations in quantities like v x (t) or S xy (t), or as a function of the Y-coordinate, grows up with time or with spacial distance.That is a typical result characteristic of solutions to the model of a chain of masses and springs.A natural question arises is the observed increase of amplitude of oscillations going to continue increase with time/distance?This question cannot be resolved by performing numerical modeling by using Bessel functions as provided by Pater equations-the range where accurate approximations of these functions is available is too narrow.We did preliminary MD simulations, performed on very long samples (up to 1000 layers; 0.2 µm in size), but still were not able to find clear conclusions.The nature of oscillations changes at long times (distances).The phenomenon is related to the harmonic case as well, when n = 2, but the effects are stronger for n increasing well over 2 and for increasing values of pressure at the surface.Some of our preliminary studies suggest, however, that under high pressures, a catastrophic, chaotic scenario of stress propagation is realized, i.e. the wave front changes abruptly its own shape.On another hand, there are theoretical arguments that the amplitude of similar oscillations decreases with time/distance [10].This subject requires separate systematic research, theoretical and based on MD simulations.Perhaps the new experimental techniques available [1] could also provide physical evidence for the nature of these phenomena.

Sound velocity
A peculiar property of Pater's solutions is an infinite speed of sound in the limit of very small displacements, as noticed by Schrödinger [21,22].The displacement of the mth layer, u m , is proportional to J 2m (2θ), which for small arguments can be approximated by u m ∼ 1/(2m)!• (Ωt) 2m .Hence, u m is non-zero (albeit it may remain very small) for any m and for any arbitrarily small time t > 0.
It is difficult to formulate a precise, rigorous criterion for determining the time when stress entering the volume starts to increase at a position of a given layer, or when displacement starts to change.That increase in time is initially slow.We are able to show that u m (t), as determined from MD simulations, can be accurately described by (Ωt) 2m functions of time, however for not more than 4 layers, only.
We can determine the speed of sound by using analytical equations of Pater and numerical methods of finding out when a certain threshold displacement of u m is reached for various layers, to determine from that the speed of sound c.That method gives results on c dependent on the value of the critical threshold assumed.However, that dependence is not strong and in the limit of large thresholds (when u m is of the order of 1) the well anticipated limit of a classical infinite chain of masses is approached.
The red horizontal line in figure 18(b) is drawn at 0.1 of the applied surface pressure and the points of crossing it by S xy (t) curves have been used to determine the speed of sound (this is a transverse sound component, with propagation in direction [111]).Full circles in figure 20(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 Ω, and to the distance d between layers, c = d • Ω.In our case, d = a/ √ 3, where lattice constant a is 3.56 Å.For comparison, empty circles in figure 20 That method gives us results on c(n) nearly identical to those determined with the first method when n < 2, with some small and interesting differences for larger n.It is simpler and faster than the first one.However, it must be used with caution in situations of non-linear materials, unless we understand well the mechanisms of wave propagation in a studied medium.
We think that the velocity of stress propagation is influenced largely by the asymmetry of the potential well, as shown in figure 8.The time needed for a particle to move on the left of the potential minimum is different from the time for traversing the right part: for n > 2 the time on the right is shorter than the time on the left.This asymmetry is well seen in figure 20(b), where velocity normalized by the period of oscillations T for a 2-layer case is shown for a range of n-values.Hence, we argue that the time on the right determines the speed of sound.
The point of potential minimum is marked as (x 0 , t 0 ) in 12. Position x 0 is known from (23).Since we have no analytical expression on x(t), we cannot 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 velocity, as in figure 20(b).The results on dependence of 4 • t 0 /T on n are shown in figure 20(c), together with x 0 /A H .
When we correct the sound speed computed from equation c = d • Ω by a factor 1/(2 − 4t 0 /T) (this is like assuming that the time T/2 − t 0 should replace T/4 on the period of harmonic oscillator), the results become very close to those determined from the measurements of the front wave, as shown by empty squares and the blue curve (for guiding eyes only), in figure 20(a).
The determined two-layer oscillation frequency ω in the displacement method for n = 2 is 26.9 ×10 12 s −1 , and the speed of the transverse sound is 55.404 Å ps −1 .Both these values are about 2 times higher than those found by us for steel 310S, when Artur EAM inter-atomic potential is used, where we have 14.6 ×10 12 s −1 and 30.29 Å ps −1 , respectively.The results of MD simulations for steel are in agreement with the ab − initio computations [48] and with experimentally determined phonon energies in Zarestky and Stassis [49] or steel [50,51].That means that the parameter ϵ 0 used as the amplitude of the potential energy well ought to be rescaled down about 4 times to reproduce properly the results found in steel.

Summary and conclusions
The method has been described for creating a table-style non-linear inter-atomic potentials between NN atoms in FCC crystals for the 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 us understand the dynamics of stress propagation into the bulk of crystal, consisting of many 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 interatomic potential.The model implies the existence of strong oscillations of stress, displacement and velocities as a function of time and position within the volume of elastic medium.We found also that in non-linear cases Pater's model remains qualitatively valid as well, since a quasi-linear relation between the wave vector and frequency of oscillations is preserved.That allowed us to determine that the speed of sound is a strongly decreasing function of n.In case, however, when non-linearity 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.
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 and energy losses are present.Oscillations in such cases can be observed in surface pressure 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.
The developed methodology of performing MD simulations and subsequent data analysis allowed us to achieve unprecedented accuracy, with results in very good agreement with the theoretical, analytical model for n = 2.
An unresolved question is how adding more degrees of freedom (in particular, allowing the movement of atoms in the 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 the stability of the crystal lattice?Will finite temperature change the effective inter-atomic potential, and how will it influence results of the model of chain of masses and springs?
) and (b) show the position of atoms from two perspectives.

Figure 1 .
Figure 1.The orientation of FCC sample and distance between NN atoms.(a) shows mutual atoms' arrangement as seen on x-y plane and (b) on x-z plane.(c) shows the sample, with the same orientation as in figure (a).The distance between green atom on the lower plane and 3 NN atoms on the next upper plane (blue atoms) is a/ √ 2, and it is the same as the distance between blue atoms, shown by the red arrow in (b).That distance is the translational period of potential in X-direction.

Figure 2 .
Figure 2.Potential energy change Ep(x) of a single atom in a layer displaced for x with respect to the other layer, for several values of the non-linearity parameter n, from n = 1.5 to n = 3.5, as shown by labels.On the top of sets of data points for any given value of n, tiny green lines are drawn by using Ep(x) = 2ϵ 0 • (x/ξ) n , with one pair of parameters ϵ 0 = 4.240 eV, and ξ = 2.006 Å, common for all values of n.The inset illustrates the method of computation of Ep in LAMMPS: dots symbolize positions of atoms.There are three groups of atoms created: green, red and blue; each of them belonging to their own group.The blue group of atoms is moved in the X-direction for finding change in potential energy, which is computed for red atoms only.

Figure 3 .
Figure 3. Dependence of oscillation frequency ω on the exponent of potential, for displacements from 0.025 to 0.2 Å, as shown in the legend.The fitting lines are functions based on equation (18), ω(n) = a/Ψ (n) • (∆X/b) (n−2)/2 , where parameters a and b are common for all curves.

Figure 4 .
Figure 4. Comparison between computed displacement curves (normalized by the amplitude of oscillations A), shown by different colors, and the data reconstructed from Fourier components (up to 13th harmonic, all odd), shown by tiny green lines.The results are shown for the following values of n: 1.0, 2.0, 3.0, 5.0, 10.0, 20.0, and 50.0.The period of oscillations is 2 12 = 4096 time steps.When n = 1, X(t) can be constructed exactly by parabolic curves.When n → ∞, ∆X(t)/A tends to a triangle-like function.

Figure 5 .
Figure 5. Fourier coefficients a k for periodic functions of ∆X(t)/A as shown in figure 4, with ∆X(t)/A = ∑ ∞ i =0 a 2i+1 • sin(2π(2i + 1)t/T).All even coefficients are equal to zero.Arrows on the left indicate values of a k for n = 1 when ∆X(t)/A becomes parabolic-like.Arrows on the right show values of a k in the limit of infinite n, when ∆X(t)/A becomes triangle-like function.The lines are for guiding eyes, only.

Figure 6 .
Figure 6.Displacement normalized by its initial value (amplitude, 0.2 Å in this case) is shown as a function of time normalized by T. Additionally, velocities vx as a function of normalized time (right vertical axis), are shown.Maxima or minima of vx(t) are found exactly at t = T/4 or at t = 3T/4, as indicated by vertical red lines.On the top of vx(t) curves denoted in the legend, we draw thin green lines vx(t) computed from ∆X(t) data by using the formula: vx(t) = vmax(n) • √ 1 − |∆X(t)/A| n (taking care of proper sign, depending on time range).

Figure 7 .
Figure 7.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 ∆X/A (left vertical axis), pressure at the surface Pxy normalized by its maximum value (right vertical axis), and displacement raised to power n − 1, (∆X/A) n−1 , normalized in such a way that the curves coincide with pressure values.(∆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 9 .
Figure 9. 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 .

Figure 10 .
Figure 10.Displacement curves (normalized by the amplitude of oscillations A H ) in the Heaviside obtained by using MD simulations are shown by different colors.Tiny green lines represent the data obtained by reconstruction from Fourier components (up to 11th harmonic), all odd).The results are shown for the following values of n: 1.2, 1.5, 2.0, 2.5, 3.0, 5.0, 10.0, 20.0, and 50.0.The period of oscillations is 4096 time steps.When n → ∞, X(t) tends to a parabolic function, shown by blue curves.

Figure 11 .
Figure 11.Fourier coefficients a k for periodic functions of ∆X(t)/A H as shown in figure 10, with∆X(t)/A H = − ∑ ∞ k=0 a k • cos(2π • k • t/T).When n < 2,values of all a k are positive, while for n > 2 the odd and even harmonics have the opposite sign.When n → ∞, X(t) tends to a parabolic function, as shown by the two crossing blue curves.The calculated Fourier integrals in that limiting case give us:a k = −(8/π 2 ) • (−1) k /k 2 .Arrows on the right are drawn at levels corresponding to values of these coefficients.The lines are for guiding eyes, only.

Figure 12 .
Figure12.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 x 0 /A H . On the top of pressure data, the computed (∆X(t)/A H ) n−1 is drawn with a tiny green line.

Figure 13 .
Figure 13.Scaling of oscillation frequency ω with the 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 ω(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 ω(n) and for A H (n).

Figure 14 .
Figure 14.Displacement of the upper layer normalized by its maximum value as a function of time normalized by the 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 ∆X(t) are observed for any values of n.The inset shows expanded displacement (not normalized by ∆Xmax) 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 (33).

Figure 15 .
Figure 15.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 the first 2 oscillation periods.

Figure 16 .
Figure 16.Velocity vx 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.

Figure 17 .
Figure17.Pressure computed with Pater equations is compared with MD simulation results of Pxy(t), under 1 GPa Heaviside pressure impulse, applied to samples with N layers, where N is between 2 and 20 (as shown by horizontally aligned labels), when non-linearity parameter is n = 2.00.The MD results are drawn with (densely packed) dots over green and magenta curves representing results obtained with Pater equations.Green curves show surface pressure based on(37), computed for N = 2… 20.The horizontal red broken line at pressure −1 GPa is crossed by lines when t = T N /4, and the numerical results from Pater equations (green lines) are supposed to be valid only for the data lying above that line.Magenta lines are computed for the sample with N = 20 layers.The labels (vertically aligned) denote the number of layers over which the summation of the pressure contribution was taken into account.
2θ) + . ... However, for a sample with the number of layers equal to N (N ⩾ 2) summation must be limited to J k , where k = 4(N − 1) − 1.Hence, the oscillations of the first layer are shown by green lines for samples with N from 2 to 20.

Figure 18 .
Figure 18.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 1.The surface Heaviside impulse pressure applied is 2000 MPa, the sample is with N = 33 layers, and n = 2.00.Here Ω is 24.82 ps −1 , the same as in figure17.The virial pressure (force) is computed with Pater's equation(35), and the velocity of layers with equation(34).The red horizontal line in (b) is drawn at 0.1 of the applied surface pressure.

Figure 19 .
Figure 19.Comparison of the dynamics of stress penetration into the sample interior for three values of non-linearity 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.
(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 the non-linearity parameter deviations are significant.

Figure 20 .
Figure 20.(a) Speed of sound is determined as a ratio of the distance between stress curves for layers 2nd and 33rd, like those in figure 18(b), and the time needed for the travel of pressure wave between these layers.It 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.