A new computational setup suitable for the exploration of nonlinear effects in free propagation and dissipation of surface acoustic waves (SAWs) is developed based on the molecular dynamics (MD) simulation method. First applications of the computational model demonstrate the ability of atomistic simulations to reproduce the key features of the nonlinear SAW evolution, which are distinct from their well-known counterparts in bulk wave propagation. In particular, the MD simulations predict the increasing localization of the acoustic energy near the surface of the substrate during the nonlinear sharpening of the wave profile, which leads to the formation of a shock front with characteristic cusps in the horizontal strain and velocity profiles. The peak values of surface strain and velocity associated with the cusps can significantly exceed those of the initial wave. Some of the effects revealed in the MD simulations are outside the capabilities of continuum-level models and have not been explored so far. These include the observation of an unusual quadratic correction to the dispersion relation at short wavelengths defined by the frequency-dependent localization of SAWs near the surface of the substrate, the establishment of a new mechanism of the energy dissipation at the SAW shock front, where SAW harmonics generated at the limit of frequency up-conversion transform very effectively into clouds of phonon wave packets descending into the substrate bulk, and the generation of localized zones of plastic deformation at a substantial distance from the wave source. Overall, the MD methodology developed for atomistic modeling of free SAW propagation not only enables detailed analysis of the intrinsic properties of nonlinear SAWs and verification of theoretical models but also opens up a broad range of opportunities for investigation of acoustically induced surface processes, material modification by SAWs, and the interaction of SAWs with preexisting crystal defects and other material heterogeneities.

## I. INTRODUCTION

Surface acoustic waves (SAWs) are elastic waves that propagate along free surfaces of solid materials and are characterized by strain profiles exponentially decaying with depth under the surface. The research interest in SAWs is defined by their relevance to various natural phenomena, such as the propagation of seismic waves generated by earthquakes in the Earth's crust^{1} or long-range communication among elephants,^{2} as well as by the rapid expansion of the area of practical applications of SAWs, currently encompassing nondestructive probing of mechanical properties^{3–7} and surface defects,^{8,9} signal-processing^{10–12} and acousto-optical^{13} devices, chemical sensing,^{14,15} micro-scale manipulation of fluid flow, biomolecules, and nanoparticles in microfluidics devices,^{16–20} and “atomization” of liquid samples for mass spectrometry analysis.^{21–23}

Some of the applications, such as the analysis of crack nucleation and evaluation of fracture strength of brittle materials,^{4–7} rely on nonlinear evolution of SAW profiles leading to the formation of shock fronts and the associated increase in the peak stress values. The nonlinear transformation of laser-generated SAW pulses has been shown to facilitate the removal of sub-micrometer particles in surface cleaning,^{24} whereas the enrichment of the SAW spectra with high frequency harmonics due to the nonlinear acoustic effects has recently been demonstrated to enable dynamic coupling and acoustic activation of molecular level processes, such as surface diffusion of small atomic clusters.^{25,26}

In general, the nonlinear sharpening of the SAW profiles and shock formation are proceeding through up-conversion of the initially excited frequency components, which has much in common with similar nonlinear effects observed for bulk waves. The nonlinear transformation of SAWs, however, has certain features that are unique and are not typical of bulk waves. In particular, due to the direct relation between the SAW wavelength and the characteristic depth of the strain field localization, the nonlinear frequency up-conversion results in increasing the strain energy concentration near the surface. This, in turn, can produce a noticeable raise in the strain magnitude in the vicinity of the developing shock fronts,^{27–32} which is not observed in nonlinear propagation of bulk elastic waves. Therefore, the propagation of strong SAWs exhibits a number of physical features that cannot be elucidated based on straightforward analogies with well-developed theoretical descriptions of nonlinear propagation of bulk waves.

In view of the specific features of SAWs, a number of continuum-level theoretical and numerical approaches have been developed for treating the nonlinear effects in SAW propagation. In particular, the approaches based on multiple scale perturbation method,^{33–41} solution of the evolution equation with slowly varying wave profile method,^{42,43} and Hamiltonian formalism^{27–32,44} applied either in the frequency domain using a system of equations for interacting harmonics^{27} or in the form of an evolution equation for the SAW profile in time domain^{29,30} have been developed to describe the variation of different harmonic constituents and the overall evolution of SAW profiles in isotropic^{28–30,32–38,42,44} and anisotropic elastic solids.^{31,39–41,43}

The continuum-level numerical and theoretical investigations have significantly advanced the understanding of the nonlinear evolution of SAWs. It has been demonstrated, in particular, how the variation of different harmonic constituents leads to distortion of the wave profile and formation of a shock front in the horizontal (parallel to the surface) component of the particle velocity and, correspondingly, horizontal strain. A number of peculiarities of the nonlinear sharpening of SAWs have been revealed, including a significant difference between the profiles for the horizontal velocity (exhibiting shocks) and vertical velocity (that forms spikes rather than shocks),^{27} as well as cusping of the wave profile near the shock front, with the peak values of the horizontal velocity/strain exceeding those of the initial wave before the shock formation.^{27–32} The latter effect can be attributed to an increasing localization of the strain energy of the wave in the surface region as the wave profile sharpens and the energy is transferred from the fundamental to higher-harmonic components of the wave. The effective nonlinearity in this case becomes nonlocal even for a homogeneous material, i.e., in contrast to bulk waves in a homogeneous media, the nonlinear distortion of a particular part of the SAW profile is defined by processes occurring in other parts of the wave, at different depths under the surface.^{27,29}

Despite the demonstrated ability of the continuum models to provide important insights into the evolution of a SAW profile during the wave propagation in a nonlinear medium, the limitations of these models with respect to the analysis of the interaction of acoustic waves with any material heterogeneities (crystal defects, precipitates, surface roughness, and adsorbates) are quite evident and difficult to overcome through straightforward extension of the existing theoretical approaches. In particular, the wave dissipation due to the thermal conduction and internal friction is often introduced in the equations on an *ad hoc* basis, to ensure numerical stability of the solution after the shock front formation.^{27,28,30–32} The absorption of acoustic energy in real materials, however, may be sensitive to the microstructure, leading to deviations from the typically assumed quadratic dependence of the attenuation coefficients on the wave frequency,^{27,28,30–32,37} which is characteristic for viscous dissipation in fluids.^{45} Moreover, while the quadratic in strain terms in the stress–strain relation (i.e., cubic in strain terms in the elastic energy) are commonly used in continuum models to account for the nonlinear elastic properties of the material, higher order nonlinearities may become important at large strains, close to the levels that cause damage or produce modification of material microstructure. Finally, the representation of the material modification itself, as well as the acoustic activation of various surface processes, such as diffusion,^{25,26,46,47} desorption,^{48–51} and heterogeneous catalysis,^{51–55} requires *a priori* knowledge of the acoustic coupling to different elements of material microstructure or surface species and is far beyond the current capabilities of continuum modeling.

Given the limitations of the continuum approaches discussed above, the classical atomistic molecular dynamics (MD) technique based on the numerical integration of the equations of motion for all atoms in the system^{56–58} presents an attractive alternative for investigation of free nonlinear SAW propagation, dissipation, and coupling to surface defects, heterogeneities, and adsorbates. The main strength of the MD method is that it does not require any assumption about the processes taking place in the systems that are investigated. The only input into the model is the interatomic potential that defines all the thermodynamic, structural, and mechanical properties of the material. This characteristic of the MD technique presents a significant advantage over the continuum-level methods where all relevant processes have to be known (and described mathematically) before the simulations can be performed.

The MD simulations of bulk shock waves have indeed been instrumental in providing valuable information on the shock front structure,^{59,60} shock-induced plasticity and phase transformations,^{61–64} and interaction of strong shock waves with voids and crystal defects.^{65–67} Beyond the analysis of the steady-state properties of the shock front and fast dynamic material response to the strong shock loading, however, there has been little progress in extending the MD simulations to the investigation of the nonlinear sharpening and dissipation of the bulk stress waves during their long-term propagation and interaction with material microstructure. Similarly, while the dynamic fracture (spallation) caused by the reflection of a strong pressure pulse from a free surface of a sample has been investigated in a number of studies,^{67–69} more “gentle” interactions of the bulk acoustic waves with surfaces causing atomic-level structural rearrangements, diffusion enhancement, or desorption have not yet been systematically explored. The main factors that have been limiting the application of the atomistic modeling in this area are the severe limitations on the time- and length-scales accessible for MD simulations that typically do not exceed several nanoseconds and hundreds of nanometers. These limitations are further magnified in the case of SAWs, where a realistic simulation of free nonlinear propagation of a SAW requires the dynamic representation of a several-wavelength-deep surface region of the substrate, while the size of the computational system along the direction of the SAW propagation should increase with increasing distance covered by the wave.

In this paper, we alleviate some of the limitations of the MD technique through the design of computational methodology suitable for the exploration of nonlinear effects in the propagation and dissipation of SAWs. The general idea of the computational method developed for the atomistic simulations of free propagation of nonlinear acoustic waves is briefly outlined below, in Sec. II. The constraints imposed on the choice of parameters of the computational setup by the peculiarities of the nonlinear SAWs are discussed and illustrated by simulation results in Sec. III. The conditions for the breakdown of the conventional continuum description of weak (linear regime^{10}) SAWs and the effects related to the discrete atomic nature of the material structure are discussed in Sec. IV. The nonlinear transformations of wave profiles, formation of shocks, localization of the acoustic energy near the surface, and rapid dissipation of high-frequency harmonics generated via nonlinear frequency up-conversion are discussed for strong (nonlinear regime) SAWs in Secs. V and VI. The ability of a strong SAW to produce quasi-periodic structural modification or “damage” in the surface region of the substrate along the path of its propagation is illustrated in Sec. VII. The results of the simulations are summarized and the opportunities provided by the computational model for the exploration of new applications of SAWs for surface characterization, modification, and acoustic processing are discussed in Sec. VIII.

## II. COMPUTATIONAL “SYNCHROTRON” FOR ATOMISTIC SIMULATIONS OF NONLINEAR ACOUSTIC WAVES

The general idea of the computational method developed for MD simulation of the evolution of an acoustic wave profile during the free propagation of the wave in a nonlinear medium is illustrated in the schematic shown in Fig. 1. At first glance, the schematic resembles a synchrotron for particle acceleration. Drawing on this analogy, the generation of the initial wave at the source in the computational “synchrotron” corresponds to the particle injection in a real synchrotron, the gradual evolution of the wave profile and shock formation during the propagation away from the source corresponds to the particle acceleration by the electric field, and the material modification by the shock wave generated in the computational “synchrotron” corresponds to the ejection and utilization of the accelerated particle in applications. While the general setup depicted in Fig. 1 can be used in simulations of both bulk and surface acoustic waves, in this paper, we focus our attention on SAWs only. The implementation of each of the three components of the computational “synchrotron” is briefly outlined below.

(1)

**“Injection”**: Depending on the purpose of the simulation, a continuous sinusoidal wave with a given frequency or a wave packet representing a SAW pulse generated by an external optical or electromechanical stimulus can be generated in the MD computational system at the start of the simulation. As explained in Sec. III, the generation of the wave packets or periodic waves can be done by assigning instantaneous displacements and velocities to all atoms in the system according to the corresponding solutions of continuum-level equations derived for the corresponding linear (harmonic) SAWs.(2)

**“Synchrotron”**: Once the initial wave is generated, the gradual evolution of the wave profile during free propagation in a nonlinear medium can be simulated in a computational setup adopting periodic boundary conditions in the direction of the wave propagation. In the case of a periodic wave, the length of the computational system in this direction must be a multiple of the wavelength. This approach makes it possible to perform analysis of the long-term evolution of the wave profile without the need for extending the size of the computational system to match the distance covered by the wave. Since the periodic boundary conditions are applied along the direction of the wave propagation, the increasing time in the simulation corresponds to the linearly increasing distance from the source of the wave under typical experimental conditions. The utility of the periodic boundary conditions for simulation of the long-term wave propagation is conditioned by the absence (or negligible extent) of heating or structural modification of the substrate by the propagating wave, so that the material can be effectively “reused” for simulation of further propagation of the wave circulating in the system. The nonlinear sharpening of the wave and the generation of a shock front can lead to the onset of a rapid SAW energy dissipation, as discussed in Sec. V. While the energy deposition due to the wave dissipation may be offset by maintaining a constant temperature throughout the substrate with a specially designed spatially resolved thermostating procedure, any irreversible structural material modification or “damage” caused by the propagating wave would immediately render the “reuse” of the damaged material in the “synchrotron” mode of the simulation impossible.(3)

**“Ejection”**: When the strain magnitude at the shock front, increasing due to the nonlocal nonlinearity of SAWs,^{27}reaches the levels sufficient for generation of irreversible material modification, a realistic modeling of the wave propagation through the pristine material away from the source can no longer be done with periodic boundary conditions applied in the direction of wave propagation. In this case, the periodic boundary conditions should be switched off and, if needed, new material could be added and equilibrated in front of the propagating wave. In order to extend the time scale of the simulations and avoid the linear increase of the number of atoms in the system with simulation time, the material behind the propagating acoustic pulse can be removed, so that the MD computational system of a fixed size remains stationary in the reference frame of the moving front of the wave. This approach is similar to the “moving window” technique^{70}that has been successfully applied in investigations of bulk shock waves in gases,^{71}liquids,^{72}and crystals,^{60,63,73}as well as in the analysis of detonation front instabilities.^{74}

Note that the addition of new material and the “moving window” approach are only needed when the long-term evolution of the wave after the generation of damage is of interest. In studies focused on the conditions for the onset of surface damage and the nature of surface modification, the simulations can be completed before the wave front reaches the outer boundary of the computational system. This scenario is exemplified by a simulation discussed in Sec. VII, where a gradual evolution of a SAW in the computational “synchrotron” produces a shock wave that is sufficiently strong to cause the generation of several zones of localized plastic deformation (emission of partial dislocations) in the surface region of the substrate. Alternatively, when long-term structural changes behind the front of a shock wave are of interest, but the evolution of the wave profile is not, the acoustic impedance matching boundary conditions^{75–77} may be used to simulate a non-reflective exit of the wave from the computational cell. This approach has been successfully applied to bulk waves generated in simulations of laser–material interactions^{76–78} and shock-induced chemical reactions,^{79} although its extension to SAWs still calls for an advancement of the computational methodology fully accounting for the specific properties of the surface waves.

## III. DETAILS OF MD COMPUTATIONAL SETUP

The main goal of the simulations reported in this paper is to explore, using the general computational approach outlined in Sec. II, the new opportunities provided by the atomistic modeling in the investigation of linear and nonlinear propagation of SAWs. The computational systems considered in the simulations, therefore, are chosen to be simple but representative of different classes of real materials. In particular, the simulations are performed with two pairwise interatomic potentials that reproduce materials with positive and negative values of the parameter of local nonlinearity^{80} (e.g., metals and fused silica, respectively), so that qualitatively different evolution of the wave profile (e.g., formation of shock front at either ascending or descending branches of the strain profile of initially harmonic SAWs) could be studied.

The material with a positive value of the parameter of local nonlinearity is described by the common Lennard–Jones (LJ) potential, with parameters σ and ɛ defining the length and energy scales of the interatomic interactions and a cutoff function^{81} used to ensure that the interactions vanish at a cutoff distance of *r _{c}* = 3σ. The negative value of the parameter of local nonlinearity is reproduced, somewhat counterintuitively, by a harmonic (parabolic) pair potential describing the interactions between nearest neighbor atoms. For a one-dimensional (1D) chain of particles interacting via the harmonic potential, the dependence of the elastic energy density on strain is also ideally parabolic. Consequently, the propagation of waves in the 1D chain is perfectly linear and does not result in the generation of higher frequency harmonics or other nonlinear effects. In the bulk of a three-dimensional (3D) crystal described by the harmonic pair potential, however, the anharmonic quadratic terms (i.e., second-order elastic constants) in the stress–strain dependence do appear due to purely structural reasons, and the sign of these terms is opposite to that in the model material described by the LJ potential. Specifically, the tensile deformation leads to the effective stiffening of the material described by the harmonic potential while the material softens if the LJ potential is used, as illustrated in Fig. 2. This contrasting behavior makes it possible to perform an initial exploration of the effect of the sign of the acoustic nonlinearity on the evolution of SAW profiles with simple and computationally efficient pairwise potentials.

Since most of the results discussed in this paper are obtained with the LJ potential, all the physical quantities are presented in reduced units scaled by the LJ length and energy parameters, σ and ɛ, and the atomic mass *m*. In particular, the time, velocity, and temperature are expressed in units of $\tau =(m\sigma 2/\epsilon )1/2$, $\upsilon =\sigma /\tau =(\epsilon /m)1/2$, and *T* = ɛ/*k _{B}*, respectively, where

*k*is the Boltzmann constant. For the harmonic potential, the force constant (second derivative of the potential function at the equilibrium interatomic distance) is chosen to match that of the LJ potential, while the equilibrium interatomic distance is chosen so that the zero temperature lattice parameter of the face centered cubic (fcc) crystal structure is the same as that predicted with the LJ potential.

_{B}The MD simulations are performed for a computational system shown schematically in Fig. 3. It represents a surface region of a substrate with a face centered cubic (fcc) crystal structure and (001) crystallographic orientation of the surface. The periodic boundary conditions are applied along the *x* and *y* directions, parallel to the free surface of the substrate. In the direction of the SAWs propagation, the minimum system size, *L _{x}*, is defined by the spatial periodicity of the simulated SAW and varies between 12 and 1200 fcc unit cells (corresponding to about 19σ and 1860σ in a substrate equilibrated at zero temperature) in simulations with corresponding values of SAW wavelength, λ. The constraint on the minimum size of the computational system in the

*y*direction, parallel to the front of the propagating SAW, is imposed by the use of the periodic boundary condition, which requires the size of the system to be at least twice larger than the cutoff distance of the interatomic interaction potential, i.e.,

*L*> 2

_{y}*r*. To minimize computational time,

_{c}*L*is chosen to be close to the minimum value,

_{y}*L*= 7.8σ, which corresponds to 5 fcc unit cells.

_{y}The choice of the depth of the computational system, *L _{z}*, and the type of the boundary conditions suitable for the bottom of the simulated system requires special attention. In contrast to earlier MD simulations of SAWs,

^{47,82,83}where close-to-sinusoidal wave profiles are generated and maintained through dynamic boundary conditions applied at

^{82,83}or very close to

^{47}the surface, the simulation of free nonlinear wave propagation requires the dynamic representation of a much deeper part of the substrate. The elastic field of a SAW decays exponentially into the bulk of the substrate (oscillations superposed on the exponential envelope may appear in anisotropic crystals), and the energy carried by the SAW is largely concentrated within a surface layer with a thickness on the order of the SAW wavelength, λ. One may expect, therefore, that a depth of the MD simulation system exceeding one or two wavelengths would be sufficient for reproducing the unperturbed propagation of the SAW in an infinitely thick substrate. The first MD simulations performed for systems with

*L*ranging from 2λ to 4λ, however, reveal a rather strong influence of the finite depth of the computational system on the wave propagation, even though the wave energy density decreases at these depths down to less than 10

_{z}^{−4}–10

^{−7}of that at the upper surface.

In the simulations where a free boundary condition is applied at the bottom of the computational system, a strong linear coupling between the top and bottom surfaces is found to create an effective channel of energy exchange between the SAWs propagating along the two surfaces. As a result, the wave energy is periodically transferred from one surface to another in a manner similar to the energy exchange between two coupled identical oscillators, as can be seen in Fig. 4(a). For a sample thickness of *L _{z}* = 2λ, the period of the energy exchange corresponds to the distance of about 39λ covered by the wave. Certainly, the efficiency of coupling between the SAWs propagating along the two surfaces can be reduced by increasing the sample thickness, leading to an increase in the energy exchange period. This effect is illustrated in Fig. 5, where different rates of the energy leakage from the original wave generated in the top surface region to the wave induced at the bottom of the computational system are shown for different values of

*L*equal to 2λ, 3λ, and 4λ. Despite the marked decrease in the rate of the energy exchange with increasing depth of the system, the increase of the bottom wave amplitude on the time scale of the simulation is still significant at

_{z}*L*= 4λ and, eventually, all the energy of the original wave will be fully transferred to the energy of the wave propagating along the opposite side of the system.

_{z}A simple and radical solution to the problem of the energy leakage from a SAW propagating along the top surface of a computational system of finite depth is to apply a rigid boundary condition at the bottom of the computational system, thus eliminating the possibility of the excitation of the second SAW. Indeed, as shown in Fig. 4(b), the application of the rigid boundary condition results in a complete disappearance of the energy leakage from the SAW propagating along the top surface region of the computational system. The elimination of the resonant coupling to the bottom surface of the system, however, does not mean that the computational setup with the rigid boundary condition is free of artifacts introduced by the finite depth of the system. In particular, one of the important consequences of the finite depth of the system (with either rigid or free boundary conditions applied at the bottom) is the appearance of an extra frequency dispersion. Indeed, since the characteristic depth of the SAW energy localization is proportional to the wavelength, the presence of the opposite boundary has a different effect on waves of different frequencies. It is well known that, for nonlinear wave transformation, the phase mismatch between the interacting harmonics in a dispersive medium can drastically change the resultant wave profile. It should be verified, therefore, that the depth of the system is sufficiently large for avoiding the introduction of a strong dispersion and ensuring a realistic representation of the propagation of SAW on a substrate of infinite depth. To mitigate the effect of the finite depth of the computational cell on the wave propagation and generation of higher-frequency harmonics, three simulations are performed for the same SAW with λ = 155σ and different depths of the computational system, *L _{z}*, equal to 2λ, 3λ, and 4λ. The evolution of the second harmonic of the wave, shown for the three simulations in Fig. 6(a), demonstrates a substantial effect of the depth for the transition from 2λ to 3λ, while the results predicted with 3λ and 4λ are almost identical.

Based on these observations, and given the fact that an increase in the depth of the simulated system leads to the corresponding increase in the computational cost of the atomistic simulations, the depth of 3λ and the rigid boundary condition at the bottom of the computational system are chosen in this work as a reasonable compromise between the minimization of the artifacts related to the finite depth of the system and ensuring computational efficiency of the simulations. The largest system simulated in this work has dimensions of *L _{x}* ×

*L*×

_{y}*L*= 1860σ × 7.8σ × 5580σ (1200 × 5 × 3600 fcc unit cells) and consists of 86.4 × 10

_{z}^{6}atoms. It should be noted that the conclusion on the minimum depth requirement of

*L*≥ 3λ is not general and is discussed above for a substrate described by the LJ interatomic potential and having the (001) surface orientation. For this substrate, the dependence of the displacement amplitude at

_{z}*x*= 0 is shown in Fig. 6(b). The minimum depth

*L*required for an adequate simulation of free SAW propagation can be determined from a condition that the maximum displacement amplitude at the bottom of the computational system stays below 10

_{z}^{−2}

*A*, where

*A*is the wave amplitude (see the Appendix for further details). The displacement profiles, defining the characteristic depth of the energy localization, have a strong dependence on the elastic properties of the substrate material and the crystallographic orientation of the surface. Consequently, the minimum depth requirement may change significantly for a different material system and surface orientation.

The initial SAWs generated in the present study are continuous periodic (sinusoidal) unidirectional Rayleigh waves propagating along the *x* axis. The investigation of the propagation of the periodic SAWs is a necessary first step in the analysis of the behavior of more complex wave profiles that can be represented as linear superpositions of harmonic waves. The consideration of continuous periodic waves has an additional advantage of enabling the analysis of the long-term evolution of the wave profile in a computationally efficient **“**synchrotron” approach described in Sec. II. This approach involves the application of periodic boundary conditions in the direction of the wave propagation, and the system size in this direction, *L _{x}*, should be a multiple of the wavelength of the fundamental harmonic of the generated wave, i.e.,

*L*=

_{x}*n*λ. For the sake of computational efficiency, we choose

*n*= 1 in all simulations reported in this paper.

The initial sinusoidal SAWs are generated by adding displacements and velocities to the instantaneous positions and velocities of all atoms in thermally equilibrated substrates according to the analytical expressions derived for Rayleigh waves. The extra velocities and displacements are added at the start of MD simulations, i.e., at time *t* = 0, and the waves are allowed to evolve freely during the simulations. In this work, the simulations are performed for SAWs propagating along the [100] direction of the (001) plane in a substrate with the fcc structure. The initial displacements and velocities are chosen according to the corresponding linear solution for a harmonic Rayleigh wave in a cubic crystal. The explicit equations based on the analysis reported in Ref. 84 are provided in the Appendix. The values of the elastic moduli of the LJ model material entering the equations for the Rayleigh wave are calculated based on the statistical fluctuation method.^{85} The value of SAW velocity predicted by the analytical equations can be verified in test simulations, adjusted in the analytical equations if needed, and used to represent the wave evolution in a reference frame moving with the running wave. This representation is convenient for comparison with analytical solutions, which are often formulated in a coordinate system moving with the wave. Note that it is critical to ensure an accurate match of the displacements and velocities added to the atomic motions at *t* = 0 and the actual distributions for a running harmonic wave of a given frequency and amplitude. Any deviation from the correct distributions may lead to the excitation of a wide spectrum of spatial and frequency harmonics propagating in various directions in the form of both SAWs and bulk waves.

## IV. PROPAGATION OF WEAK SAWs: SPATIAL DISPERSION

The implementation of the MD model for simulation of SAWs, described above, enables computational analysis of some of the properties of SAWs that are difficult to study analytically. In this section, we consider weak (low amplitude) SAWs, for which nonlinear distortions of the wave profiles are relatively small, and focus our attention on the conditions leading to the breakdown of the continuum description of wave propagation due to the effects related to the discrete atomic nature of the material structure. In particular, we consider manifestations of the spatial dispersion of high-frequency SAWs with wavelengths that exceed the lattice parameter of the crystalline structure of the substrate by only one to two orders of magnitude.

It is well known that, in the non-dispersive case, the amplitude of the second harmonic of a plane sinusoidal wave is initially proportional to the propagation distance. The presence of dispersion, however, introduces a phase mismatch between the interacting fundamental and higher harmonics, leading to a distinct periodic beating of the wave amplitude. For example, the superposition of fundamental and second harmonics should produce beating with a spatial period of $L1\u22122b$ defined by the difference between the wavevector of the fundamental harmonic, $k1$, and half the wavevector of the second harmonic, $k2$,

This periodic beating reflects the fact that a wave propagating from a nonlinear double-frequency source (the second harmonic appears naturally in the course of nonlinear SAW propagation) a distance of $L1\u22122b/2$ acquires a relative phase mismatch equal to *π* due to the difference between the phase velocity of the fundamental harmonic, $c1=\omega 1/k1$, and that of the second harmonic, $c2=2\omega 1/k2$. The beatings in the amplitudes of the second harmonics observed in MD simulations of freely propagating SAWs of three different wavelengths are shown in Fig. 7(a). In these simulations, the displacement amplitude of the initially sinusoidal waves is chosen to be ∼0.019σ, which is sufficiently small to prevent the formation of a shock front during the time of the simulations. As discussed below, in Sec. VI, the appearance of a shock front would drastically increase the wave dissipation rate, which would also affect the evolution of the amplitudes of the wave harmonics.

Wavelength λ [σ] (unit cells) | 19 (12) | 39 (25) | 78 (50) | 155 (100) | |

Surface waves | 620 | 3440 | 15 000 | 61 800 | |

Distance $L1\u22122b$ [σ] | Longitudinal waves | 3070 | 35 200 | 301 000 | |

Shear waves | 1030 | 9380 | 74 300 |

Wavelength λ [σ] (unit cells) | 19 (12) | 39 (25) | 78 (50) | 155 (100) | |

Surface waves | 620 | 3440 | 15 000 | 61 800 | |

Distance $L1\u22122b$ [σ] | Longitudinal waves | 3070 | 35 200 | 301 000 | |

Shear waves | 1030 | 9380 | 74 300 |

The strong decrease in the beating period with decreasing wavelength, apparent from Fig. 7(a), is related to an increasing manifestation of the spatial dispersion at shorter wavelengths. For the propagation of the second harmonic in systems with spatial dispersion related to the discrete atomic structure of the material, the common analytical solutions are based on the approximation that the lowest dispersion correction to the dispersionless dependence $\omega =ck$ is cubic in *k*, namely, $\omega \u2248ck\u2212\alpha k3$, where normally $\alpha $ is positive. A quadratic in *k* dispersion correction is absent, as the quadratic and other even-order terms correspond to dissipation rather than dispersion. In particular, a cubic in *k* dispersion term naturally appears for smaller wavelengths in the case of a periodic chain of masses interacting via elastic springs. By additionally taking into account the lowest (quadratic in amplitude) nonlinear term, the wave evolution in this system is described by the Korteweg–de Vries equation.^{86} Its solution (e.g., Ref. 87) predicts that, in the absence of dissipation, the maximum displacement amplitudes of the second harmonic should be the same for waves with fundamental harmonics of different $k1$ but identical amplitudes of the initial displacement. This prediction agrees fairly well with the results shown in Fig. 7(a).

In view of the usually cubic in *k* dispersion correction to $\omega (k)$ for the initial deviation from the non-dispersive propagation [which is equivalent to the cubic in $\omega \u2248ck$ correction to the wavevector $k(\omega )$], the spatial frequency $2\pi /L1\u22122b=|2k1\u2212k2|$ of the second harmonic beatings should also be cubic in $k1$. This cubic dependence is indeed found in MD simulations performed for bulk longitudinal (compressive) and shear waves of three different wavelengths λ. Note that, in the case of the shear waves, the second harmonic does not naturally appear during the simulation of the nonlinear wave propagation, and it is generated at the start of the simulation together with the fundamental harmonic. As can be seen from Fig. 7(b), the spatial period of beatings between the fundamental and double-frequency harmonics follow the expected cubic dependence, $L1\u22122b\u221d\lambda 3$, for the bulk waves with λ equal to 39σ and 78σ. An even stronger dependence observed for very large $k1$, when λ decreases down to just 12 fcc unit cells (see Table I), indicates that the contribution of higher order (beyond cubic) dispersion correction terms becomes significant. The results of MD simulations of SAWs, presented in Fig. 7 and summarized in Table I, however, demonstrate a behavior strikingly different from that of bulk waves. In particular, the doubling of the fundamental harmonic wavelength from $\lambda =77.5\sigma $ to $\lambda =155\sigma $ leads to an increase in the beating period $L1\u22122b$ by a factor of $22=4$ instead of a factor of $23=8$ expected for the cubic dispersion correction.

At first glance, this quadratic scaling looks counter-intuitive since, as mentioned above, the quadratic correction to the relationship between the frequency and wavenumber is normally associated with dissipation rather than dispersion. The unusual characteristic of the spatial dispersion of SAWs, however, can be attributed to the wavelength dependence of the degree of localization of a SAW in the surface region of the substrate. The dispersion of a SAW arises due to scattering not only from the atoms localized in the upper atomic layer but also from periodically spaced vertical “screens” with the effective vertical size defined by the depth of the SAW penetration, as shown schematically in Fig. 8. The elastic field of a SAW penetrates into the bulk of the substrate material to the depth on the order of λ. Thus, the size of the vertical “screens,” as well as the corresponding effective number of discrete scattering elements (atoms) in each screen, are proportional to $\lambda \u221dk\u22121$. As a result, a wave with twice shorter wavelength “feels” twice smaller number of discrete scattering elements. In order to account for this effect, the conventional $\u223ck3$ dispersion correction due to the scattering by periodically spaced “screens” should be multiplied by a factor of $k\u22121$, reflecting the decrease in the effective number of scattering elements with increasing *k*. Consequently, for SAWs, the dispersion correction becomes quadratic, $\u223ck\u22121k3=k2$, as indeed is observed in the MD simulations of SAWs with sufficiently large wavelengths (exceeding 20–30 fcc lattice parameters) [Fig. 7(b)]. For very large *k* (when λ decreases down to a dozen of lattice parameters), the role of the very upper atomic layers becomes dominant (closer to the 1D case), so that effective vertical aperture of the “screens” does not continue to decrease as $k\u22121$. As a result, for such extremely short SAWs, the initial $\u223ck2$ dependence of the dispersion correction gradually becomes steeper. Moreover, the contribution of higher-order dispersion correction terms may become non-negligible, as is also discussed above for the bulk waves with λ = 19σ.

This difference in the effect of the spatial dispersion observed for relatively weak bulk waves and SAWs discussed above provides a clear example illustrating the implications of an important intrinsic characteristic of SAWs, namely, the strong frequency-dependent localization of SAWs near the surface of the substrate. In Secs. V–VII, we will discuss a number of additional effects that stem from this unique characteristic of SAWs and manifest themselves at higher wave amplitudes, when nonlinear distortions of SAW shapes become significant.

## V. SAW ENERGY CONCENTRATION AND DISSIPATION AT THE LIMIT OF FREQUENCY UP-CONVERSION

One of the most attractive characteristics of the atomistic computational setup developed in this work for simulation of SAWs is its ability to naturally reproduce the generation of higher frequency harmonics during the propagation of an initially sinusoidal wave in a nonlinear medium, leading to the transformation of the wave profile and localization of the acoustic energy near the surface. The evolution of the wave profile leading to the formation of a shock front is discussed in detail below, in Sec. VI. In this section, we focus on the elastic energy localization during the nonlinear propagation of SAWs and on the related phenomenon of rapid dissipation of the wave due to the conversion of high-frequency SAW harmonics to thermal phonons.

The effect of the concentration of the acoustic energy of SAWs in the course of their nonlinear propagation is illustrated in Fig. 9, where the kinetic energy distribution is shown for an MD simulation of a SAW with an initially sinusoidal wave profile, a wavelength of 930σ (600 fcc unit cells), and a strain amplitude of 2%. The kinetic energy per atom is averaged over the wave period and is plotted as a function of depth under the surface and distance traveled by the wave. For clarity, the MD simulation is performed for an initially cold sample (zero temperature) to eliminate the need for filtering out the thermal noise and extracting the collective velocities and displacements of atoms related to the propagating SAW. In the absence of thermal noise, the redistribution of the elastic energy of the wave closer to the surface of the substrate is clearly visible in Fig. 9(a). This increasing energy concentration near the surface is directly related to the nonlinear frequency up-conversion during the propagation of a sufficiently strong SAW,^{27–29} which leads to an enrichment of the wave spectrum with higher harmonics^{25} localized closer to the surface. Indeed, the concentration of the kinetic energy density at the surface of the substrate has also been predicted in a theoretical analysis of the harmonic generation in nonlinear propagation of a Rayleigh wave in an isotropic solid.^{28}

Due to the large strain amplitude and short wavelength, a pronounced increase in the concentration of the SAW energy near the surface of the substrate is observed within a relatively short wave propagation distance of ∼10λ. As discussed in Sec. VI, the characteristic time (or propagation distance) required for the frequency up-conversion leading to a significant modification of the SAW profile, generation of a shock front, and wave energy concentration near the surface scale linearly with the fundamental wavelength of the wave and inversely with the initial acoustic strain amplitude. Thus, this time/distance can be significantly longer for SAWs of much larger wavelength and smaller amplitudes typically generated in experiments.^{4–23} Nevertheless, the nonlinear sharpening of the wave profile and shock front formation has been observed in experiments performed with strong optically generated SAWs characterized by much larger initial wavelengths.^{5,88} The values of strain and characteristic spatial scales of its variation in the vicinity of well-developed shock fronts are comparable to those produced in the MD simulations, and the important processes occurring at the shock fronts and discussed below are likely to be similar in the simulations and experiments where sufficiently strong SAWs are generated. The shorter wavelength of the simulated SAWs and the rapid emergence of the nonlinear effects makes it possible to explore the key mechanisms of the SAW dissipation and SAW-induced material modification in atomistic simulations at an acceptable computational cost.

The dissipation of a strong SAW propagating in a nonlinear medium is a prime example of a phenomenon that still lacks complete theoretical understanding. While the mechanisms and rates of SAW energy dissipation at different stages of the wave profile evolution can be assumed in continuum-level models,^{27,28,30–32,37} they cannot be studied, confirmed, or discovered with these models. In spite of its inherent time- and length-scale limitations, the atomistic modeling, on the other hand, can help to reveal and quantify the new mechanisms and channels of the SAW dissipation. As an illustration, we can consider the late stage of the kinetic energy evolution shown in Fig. 9(b) with a focus on a lower energy region extending down to 2000σ, i.e., more than 2λ. While the initial increase in the SAW energy concentration near the surface of the substrate [Fig. 9(a)] can be explained by continuum models accounting for the enrichment of the wave spectrum with higher harmonics, the onset of the active emission of wave packets descending into the bulk of the substrate after a wave propagation distance of about $8\xd7103\u2212104\sigma $, clearly visible in Fig. 9(b), cannot be described by any of the existing models. After covering this distance, the processes of the nonlinear frequency up-conversion and wave energy localization practically reach their saturation. Further nonlinear transformation results in the generation of extremely high-frequency components comparable to the Debye frequency. In other words, the newly generated harmonics get close to the range of thermal phonons, and their transformation into even higher-order acoustic vibrations becomes impossible due to the atomically discrete nature of the material.

The simulations predict that the highest frequency harmonics of the SAW can effectively transform into high-frequency bulk modes, leading to the generation of phonon wave packets that channel the energy of the SAW into the depth of the substrate. The phonon wave packets moving into the bulk of the substrate appear as stripes of elevated energy in Fig. 9(b), with the slopes suggesting that the depth of their descent is increasing about twice as slow as the distance covered by the SAW along the surfaces. Additional details of the generation of the wave packets can be seen in Fig. 10, where snapshots of instantaneous distributions of the kinetic energy in the substrate are shown for three moments of time that correspond to the final stage of the sharpening of the wave profile (shown in the top panel of Fig. 10 and discussed in more detail in Sec. VI), when the process of the nonlinear frequency up-conversion reaches its limit. The active generation of high frequency wave packets propagating into the bulk of the substrate can be seen as sudden appearance of fine wavy patterns with wavelengths ranged from 2 to 7 fcc unit cell sizes as the propagation distance exceeds $104\sigma $.

Given that the minimum wavelength of the propagating acoustic phonons is 2 unit cells, the shortest visible packets (phonons) evidently correspond to the smallest wavelength of elastic perturbations with the maximal Debye frequency. As expected, an analysis of a sequence of snapshots, such as the ones shown in Fig. 10, yields group velocities of high-frequency wave packets that are much smaller than the velocity of the SAW. Moreover, the direction of the wavevectors observed for the wave packets with the shortest wavelength is about 45°–50° with respect to the surface, which means that the on-surface projections of the phonon wavevectors are about 1.5 times smaller than the total wavevector. Since the phonon wave packets are generated synchronously by the SAW, the phase velocity of the phonon profiles is ∼1.5 times smaller than the phase velocity of the SAW, as expected for the ultimately high-frequency phonons. The relatively sharp threshold-like onset of the emission of the wave packets at the limit of frequency up-conversion, as well as the fact that the simulation is performed at zero initial temperature of the substrate, suggests that the existing theoretical models of the acoustic energy dissipation may not be adequate for the description of the dissipation at the SAW shock front revealed in the MD simulations. Further systematic atomistic studies are needed to evaluate the dependence of this dissipation mechanism on the substrate temperature, material properties, and parameters of SAWs, so that an analytical description suitable for incorporation into continuum-level models could be formulated.

Note that an accurate quantitative evaluation of the wave dissipation may require a few technical adjustments to the computational setup used in the present study and described in Secs. II and III. First, the bulk waves generated at the SAW shock front could reach the bottom of the computational cell and reflect from the rigid boundary applied there. The relatively large depth of the computational cell and small magnitude of the bulk waves [note that the scale of the contour plot in Fig. 9(b) is reduced by a factor of 20 with respect to the one used in Fig. 9(a)] precludes any significant effect of the bulk wave reflection on the evolution of the SAW profiles. Nevertheless, a detailed quantitative analysis of the wave dissipation may require the application of acoustic impedance matching boundary conditions^{75–77} capable of reproducing non-reflective exit of the bulk waves from the bottom of the computational cell while recording the information on the energy carried away by the waves. Second, an adequate computational description of the thermoelastic mechanism of SAW dissipation,^{45,89} related to the spatial variation of temperature caused by the strain field of the propagating wave,^{47} relies on the ability of the model to describe the heat flow down the transient temperature gradients created by the wave. For metals, the thermal conductivity is dominated by electronic contribution, which is not accounted for in the conventional MD method. Thus, a realistic modeling of SAW dissipation in a metal substrate can only be done with a model incorporating a description of the electronic contribution to thermal conductivity of the material.^{64,65,90}

## VI. SAW PROFILE TRANSFORMATION AND SHOCK FRONT FORMATION

The increased localization of the acoustic energy of a SAW near the surface of the substrate and the onset of the rapid dissipation of the wave upon the shock front formation, discussed in Sec. V, are caused by the transformation of the wave profile during its propagation in a nonlinear medium. In this section, we consider the main features of the SAW profile evolution during the nonlinear transformation. The profiles of the horizontal surface strain $exx$ are shown in Fig. 11 for three MD simulations of free propagation of SAWs with initial strain amplitudes of 2%, 3.4%, and 4.7%. The fundamental wavelength of the waves, λ = 937σ, is sufficiently large to ensure that the velocity dispersion related to the discrete atomic structure of the material, discussed in Sec. IV, remains fairly weak. The temperature of the substrate is 0.15ɛ/*k _{B}*, which exceeds the room temperature for values of ɛ that can be obtained by rough fitting to bonding energies of most metals.

^{91,92}In contrast to the simulation performed at zero initial temperature and illustrated in Fig. 10, thermal vibrations of atoms can introduce a substantial level of noise to the strain profiles calculated based on instantaneous atomic positions. The calculation of the SAW strain profiles in finite temperature simulations, therefore, involves time and spatial averaging, mimicking that involved in macroscopic acoustic probes in real physical experiments. In particular, the displacements of substrate atoms are averaged over surface subdomains with sizes of about 3.8σ in the

*x*and

*z*directions. In addition, the values of displacements were temporally averaged with exponential moving average of 50τ in a frame of reference moving with the wave velocity.

The common trend in the evolution of strain profiles observed in all three simulations illustrated in Fig. 11 is the gradual steepening of the ascending branches of the initially sinusoidal SAWs and the eventual formation of sharp shock fronts. As briefly discussed in Sec. III and illustrated by Fig. 2, the substrate material described by the LJ potential exhibits a positive value of the parameter of local nonlinearity, which means that the elastic modulus increases upon compression and decreases upon tension. As a result, the compressive phase of the wave (negative $exx$) propagates faster that the tensile phase (positive $exx$), leading to the formation of the shock front at the ascending branch of the wave profile (Fig. 11).

The processes of the wave profile steepening and shock front formation described above are similar to those well known for the nonlinear transformation of bulk waves. The nonlinear transformation of SAWs, however, demonstrates some peculiarities that do not have direct analogies with the bulk waves. In particular, in the course of nonlinear propagation of a bulk non-dispersive wave, the points of the wave profile are merely displaced horizontally with different velocities, so that the amplitude of the developed shock does not exceed that of the initial sinusoidal wave, and only the gradient of stress (or pressure) becomes much higher as compared to the initial value. For the bulk waves and SAWs alike, the nonlinear steepening of the wave profiles and the formation of shock fronts can be described in terms of the generation of higher frequency harmonics. In contrast to the bulk waves, however, the nonlinear frequency up-conversion can result in a pronounced increase in the magnitude of the elastic strain in the vicinity of the shock front. Indeed, the formation of characteristic cusps on the compressive sides of the wave profiles, extending by more than 50% above the initial strain amplitude, can be seen in the right panels of Figs. 11(b) and 11(c). The formation of these cusps can be explained by the wavelength dependence of the characteristic depth of the SAW strain field localization, namely, the shorter wavelength harmonics of the wave are propagating closer to the surface. As a result, as the wave profile sharpens and the energy is transferred to higher-frequency components of the wave, the strain energy is increasingly concentrated in the surface region of the substrate, and the surface strain can exceed the amplitude of the initial wave. One of the important implications of the increase in the level of surface strain during the nonlinear propagation of a sufficiently strong SAW is the possibility of causing material damage at a distance from the source of the wave, as has been demonstrated in experiments on SAW-induced surface fracture of Si and fused silica substrates.^{5–8} Indeed, the cusping at the shock front in the simulation illustrated in Fig. 11(c) creates the conditions for the onset of plastic deformation just after the moment of time for which the profile is shown in the right panel, leading to the emission of partial dislocations at a time of 1002τ. This process of the generation of material damage in the form of localized zones of plastics deformation is discussed below in Sec. VII.

Note that the formation of SAW profiles with pronounced asymmetry between the compressive and tensile parts of the wave can be suppressed by the rapid dissipation of the high frequency harmonics generated at the shock front. In particular, in the simulation illustrated in Fig. 11(a), the initial amplitude of the surface strain, $exx0=2%$, is lower than those in the strain profiles shown in Figs. 11(b) and 11(c), the time of the shock wave formation is longer (see discussion at the end of this section), and the dissipation of the highest frequency harmonics generated at the shock front is sufficiently fast for preventing the formation of the cusp on the compressive side of the wave profile. Consequently, even after the shock front is fully developed, the upper/tensile and lower/compressive parts of the wave profile remain fairly symmetric, and the maximum surface strain remains close to the amplitude of the initial sinusoidal wave. The rate of the dissipation of the high frequency harmonics is sensitive to the temperature of the substrate and, for a SAW of the same wavelength and initial strain amplitude of 2%, a sharp cusp with a maximum strain of more than 4% is observed to form at the compressive side of the shock front at a four times lower substrate temperature of 0.036ɛ/*k _{B}*.

^{25}The dissipation of the high frequency harmonics at the shock front results in an overall loss of elastic energy of the wave, as can be seen from the significantly attenuated wave profile shown by the red dashed–dotted line in the right panel of Fig. 11(b) for a time of 8936τ.

The nature of the distortion of the SAW profile leading to the asymmetry with respect to tension and compression has its root in the joint influence of the material nonlinearity leading to the generation of higher harmonics and the gradual phase shift of the harmonics with respect to each other. The latter effect, i.e., the relative phase shift of the harmonics, can be caused by the spatial dispersion emerging due to the discrete atomic nature of the material structure, as discussed in Sec. IV, or due to some other specific properties leading to the dispersion of wave velocities in the nonlinear medium, e.g., the occurrence of relaxation.^{93} Moreover, the peculiarities of the structure of surface waves in crystals, where the complex-valued vertical components of the wavevector result in appearance of different phase corrections for harmonics of different orders, may produce the phase shifts and corresponding asymmetry of the shock front that are generically similar to the dispersion-induced distortions of the wave profile.^{31} In the systems considered in the MD simulations, both mechanisms, i.e., the spatial dispersion and the anisotropy-related phase shifts, may contribute to the asymmetry between the positive and negative sides of the wave profiles. For a given origin of the phase-shift-induced asymmetry of the wave profile, however, it is the sign of nonlinearity that defines the side of the shock front (compressive or tensile) on which the high-amplitude cusps appear on the horizontal strain and velocity profiles.

As discussed in Sec. III and illustrated in Fig. 2, the type of nonlinearity introduced by the LJ potential is characterized by a positive value of the parameter of local nonlinearity^{80} (i.e., the material stiffens under compressive strain). This results in the formation of shock fronts at the ascending sides of SAW profiles and the appearance of sharp cusps on the negative (compressive) sides of the shock fronts in Fig. 11. It can be expected that the opposite sign of the nonlinearity parameter, i.e., the negative sign of the quadratic correction to the linear stress–strain dependence in a substrate material, can alter the character of the nonlinear distortions of the wave profile. In particular, the shock front should form at the descending side of the wave profile, and the sharp cusp should appear at the positive (tensile) part of the shock front. These conjectures are confirmed in an MD simulation performed for a system described by a short-range harmonic potential applied to the interaction between the nearest neighbor atoms only. Despite the absence of physical nonlinearity due to the shape of the potential itself, an fcc crystal described by the harmonic potential exhibits nonlinearity of geometrical origin that has a sign opposite to that of the LJ system (see Fig. 2). Since the material represented by the harmonic potential is “indestructible” (cannot be “damaged” by any level of strain), the MD simulation in this case was performed for a SAW with a large initial strain of $exx0=8.5%$ and short fundamental wavelength of λ = 155σ, which resulted in the formation of the shock wave on a short time scale of ∼200τ despite the relatively weak nonlinearity of the substrate material. The surface strain profile after the shock formation is shown in Fig. 12. As expected, the cusp exceeding the initial strain amplitude appears at the tensile side of the shock front, and the shock front itself forms at the descending branch in the SAW profile and has an inclination opposite to that observed in Fig. 11.

The nonlinear evolution of the SAW profiles predicted in the MD simulations can be related to the results of continuum-level simulations and analytical studies.^{27–44} To facilitate the comparison, the strain profiles discussed above are complemented by the profiles of horizontal and vertical velocities plotted in Fig. 13 for the same SAWs that are illustrated in Figs. 11(a) and 11(c). The velocity profiles are shown for times that approximately correspond to those of the strain profiles shown in the right panels of Figs. 11(a) and 11(c), when the shock fronts of the waves are well developed. In a good agreement with the theoretical predictions,^{27,28,30,31} the profiles of the horizontal velocity component $\upsilon x$ exhibit characteristic sawtooth shapes with shock fronts similar to those observed for the horizontal strain in Fig. 11, while the profiles of the vertical velocity component $\upsilon z$ exhibit narrow spikes in the negative direction, i.e., away from the substrate. When comparing the profiles shown in Figs. 10–13 with theoretical predictions, it should be noted that the velocity is recorded as a function of time at a particular position in Refs. 27, 28, 30, and 31, whereas the wave profiles are plotted as a function of position at a fixed moment of time, so that the direction of the *x* axis is effectively reversed.

Taking into account the proportionality $\upsilon x\u221dexx$, the left panels of Fig. 13 once again illustrate the dependence on the wave profile evolution on the amplitude of the initial wave that is already discussed above for the strain profiles shown in Fig. 11. The shock front formation in Fig. 13(a) is not followed by the emergence of a cusp above the initial level of the horizontal velocity amplitude, while the cusped profile is generated for a wave with a higher initial amplitude. This observation is consistent with conclusions of a theoretical analysis,^{28} where the corners of the sawtooth profiles are predicted to evolve into cusps above a certain level of the wave amplitude.

Another interesting feature of the profiles predicted in the MD simulations is the asymmetry of the shapes of the strain and velocity profiles that becomes particularly pronounced at high wave amplitudes [Figs. 11(b), 11(c), 12, and 13(b)]. In particular, in contrast to theoretical predictions for isotropic materials,^{27,28,30} where the formation of identical cusps is observed on both negative and positive sides of the horizontal velocity profiles, and the spike of the vertical velocity has a symmetric shape, in the MD simulations, the cusps are generated on one side of the wave profiles, producing distinct asymmetry in the shapes of the wave profiles. The theoretical analysis that accounts for anisotropy of crystalline substrates^{31} demonstrates that the asymmetry of the wave profiles can naturally appear in crystalline substrates, where the direction of the vertical velocity spike and the side of the $\upsilon xandexx$ profiles where the cusp appears are controlled by the sign of the effective nonlinearity parameter. Moreover, it is demonstrated that the asymmetry of the wave profile may disappear for the same crystal but a different surface orientation and a different direction of wave propagation. In MD simulations, all these parameters can be controlled by the choice of the initial system and interatomic interaction, thus enabling verification of theoretical assumptions involved in the description of the generation of high-frequency harmonics and wave dissipation at large strains, where higher order nonlinear effects may play an important role.

The MD simulations also provide an opportunity to investigate the dependence of the distance the initially sinusoidal wave has to travel for the formation of a well-developed shock front on the initial wave amplitude and the wavelength. While quantitative analytical prediction of the shock formation distance remains challenging for Rayleigh waves, it has been suggested^{94} that it has the same functional dependence as the distance of the shock formation for the initially sinusoidal bulk waves,^{95,96}

i.e., the shock front formation distance $xsh$ is inversely proportional to the initial surface strain amplitude $exx0$, wavevector *k*, and effective nonlinearity parameter $\beta $.

The shock front formation distance *x _{sh}* predicted in MD simulations for an LJ substrate with a temperature of 0.15ɛ/

*k*is plotted in Fig. 14 as a function of the inverse strain $1/exx0$ and wavelength $\lambda =2\pi k\u22121$. The error bars in Fig. 14 show the uncertainty of the measured distances

_{B}*x*introduced by the lack of strict definition of the end of the shock formation process. The calculated wave profiles are distorted by the thermal noise, and, similar to experimental measurements, the slopes of the shock fronts do not exhibit ideal vertical rise due to the influence of wave dissipation. As a result, the slopes of fully developed shock fronts in the simulations vary in the range of 86°–88°, and the strict definition of

_{sh}*x*as a location where the slope of the horizontal velocity wavefront becomes 90°, used in continuum-level simulations of the wave propagation in an ideal, lossless solid,

_{sh}^{94}cannot be applied in the analysis of the MD results.

The theoretical expectation given by Eq. (2) corresponds to straight lines intersecting at the origin of coordinates. The constrained linear fitting of the data points in Fig. 14(a) omitting the point corresponding to the lowest strain amplitude yields the value of *β* = 0.85. The same value of *β* is then used for plotting the two dashed lines in Fig. 14(b), providing a reasonable description of the computational results. The typical value of the effective nonlinearity parameter $\beta $ in the case of bulk sound waves is several times greater for many crystals and amorphous materials, in the range of 4–8.^{97} At a qualitative level, the substantially smaller effective value of $\beta $ predicted in the MD simulations of SAWs can be attributed to the intrinsic feature of SAWs, namely, the frequency-dependent localization of the elastic energy of the wave harmonics discussed above. The fields of different-order harmonics are spatially separated in depth, which results in a weaker coupling between the harmonics as compared to the bulk waves. Indeed, the values of the nonlinearity parameter obtained for several homogeneous materials through a numerical solution of analytical equations derived for Rayleigh waves,^{94} *β* = 0.84 for polystyrene and *β* = 0.85–2.03 for several iron-based alloys, are in a good agreement with *β* = 0.85 estimated in MD simulations of SAWs propagating in an LJ substrate.

Note that while Eq. (2) provides a reasonable description of the data points generated in a series of simulations performed with different values of $exx0$ and λ, a substantial deviation is observed for some of the data points, suggesting that not all of the effects present in the MD simulations are accounted for in the theoretical scaling law. In particular, the deviation of the rightmost point in Fig. 14(a) toward larger *x _{sh}* can be explained by the influence of linear losses that counteract the shock front formation and become more pronounced at small strain amplitudes. We can expect that for SAWs of the same wavelength but even smaller strains, a shock front does not form at all, and the wave decays due to the linear dissipation. Note that, since the wave energy is proportional to λ

^{2}, the rate of the wave decay in the course of the nonlinear sharpening of the wave profile decreases with increasing λ, so that a shock front can still form for SAWs with relatively low $exx0$ but large λ. The deviation of the leftmost points from the dashed lines in Fig. 14(b) can be related to both stronger dissipation and dispersion of SAWs with shorter wavelengths. Indeed, we verified that for very short wavelengths ($\lambda =26\sigma and\lambda =39\sigma $), the effect of dispersion is so strong that a shock front is not developed, and the nonlinear evolution of the SAWs results in the formation of a train of soliton-like pulses. Finally, the shorter, as compared to the theoretical scaling law, shock formation distance observed in the simulation performed at the longest wavelength of λ = 1874σ and the largest strain amplitude of $exx0=4.7%$ can be attributed to the contribution of higher order nonlinear effects leading to an increase in the effective nonlinearity parameter $\beta $ at large strains approaching the material damage threshold.

## VII. MATERIAL MODIFICATION AND DAMAGE BY STRONG SAWs

The nonlinear sharpening and cusping of the wave profiles leading to the generation of surface strain in excess of the initial strain amplitude is a unique feature of SAWs that have important practical implications. Even if the initial strain in a SAW is lower than the level required for material modification or damage, the gradually increasing near-shock strains may reach the damage threshold at some distance from the wave source, causing the generation of surface or subsurface defects or the onset of surface cracking. Indeed, the acoustically induced surface cracking with strong laser-generated SAWs has been experimentally observed for Si and fused silica substrates at some distance from the laser spot,^{5–8} suggesting that the unique nonlinear behavior of SAW can be utilized for probing the ultimate fracture strength and transient dynamics of crack nucleation and growth in brittle materials. The theoretical treatment or continuum-level modeling of nonlinear effects at large strains, close to and at the levels that cause damage or produce modification of material microstructure, remains challenging. The MD simulations of SAWs, on the other hand, naturally account for higher orders of nonlinearity up to strains causing irreversible changes in the material structure and can be used for the exploration of the effects involving material modification and damage.

The ability of MD simulations to address the SAW-induced material modification is illustrated in this section by the results of a simulation where the nonlinear evolution of the wave profile leads to the emergence of localized zones of plastic deformation along the path of the SAW propagation. The simulation is performed for a SAW with an initial strain amplitude of 4.7% and a fundamental wavelength of 1874σ (1200 fcc unit cells) propagating in a LJ substrate with an initial temperature of 0.15ɛ/*k _{B}*. The initial strain magnitude of 4.7% is smaller than the threshold level for the dislocation emission, and the wave propagates without causing any material modification for four cycles of the computational “synchrotron” (see Sec. II). The nonlinear wave transformation leads to the gradual sharpening of the initially sinusoidal wave profile and the formation of a shock front with a cusp on the compressive side of the wave, similar to the evolution of the wave profiles discussed in Sec. VI and illustrated by Fig. 11. During the fifth cycle of the wave propagation through the computational cell, the shear stress generated at the shock front reaches a critical level required for the emission of dislocations from the surface of the substrate. The emission of semicircular loops of Shockley partial dislocations can be seen in Fig. 15, where three snapshots of atomic structure of the surface region of the substrate are shown, along with the corresponding wave profiles, at three moments of time corresponding to the dislocation emission. The atoms are colored according to the local structure, so that the intact fcc structure is blue, and the stacking faults left behind by Shockley partial dislocations are green. Since the dislocations are emitted along two of the close packed {111} planes that are not parallel to the

*y*axis, along which the periodic boundary condition is applied, the dislocation loops fold into layered zones where the stacking faults outlined by the dislocation lines are separated from each other by a distance defined by the width of the computational cell in the

*y*direction (5 fcc unit cells).

The plastic deformation of the surface region partially adsorbs the SAW energy and reduces the maximum strain at the shock front, as can be seen from Fig. 16, where an enlarged view of the strain profiles is provided for two moments of time, just before and after the SAW-induced emission of dislocations. The wave blunted by the plastic deformation propagates some additional distance in the elastic regime. The continuing frequency up-conversion and wave sharpening, however, result in the formation of a new cusp at the shock front, which again produces a new zone of plastic deformation. This process could be repeated several times resulting in the formation of a series of plastic deformation zones, as can be seen in Fig. 15 and in a movie provided in the supplementary material. The computational prediction of the generation of a sequence of spatially separated damage zones at a distance from the wave source can be related to the experimental observations of quasi-periodic arrays of SAW-induced cracks in brittle substrates, such as Si^{5–7} and fused silica.^{7}

When the shock front leaves the plastically deformed region, partial dislocations are retreating back due to the attraction to the surface (image force) and surface tension of the stacking fault, as can be seen by comparing snapshots in Figs. 15(a) and 15(b). While the small system size in the *y* direction limits the accuracy of the description of dislocation evolution in the simulation described in this section, the dislocation cross-slip and interaction between dislocations propagating along different slip planes can be expected to produce dislocation segments with limited mobility and create stable dislocation configurations in the subsurface region of the substrate. Moreover, the emission of dislocations generates atomic-scale roughness at the surface of the substrate, which may have important implications for chemical reactivity of the surface.^{98} Indeed, the ability of SAWs to substantially enhance the rates and selectivity of heterogeneous catalytic reactions has been demonstrated in a number of studies^{51–55} but still awaits clear mechanistic understanding.^{55,99} The roughening of the surface due to the SAW-induced generation of dislocations is a possible mechanism that may be responsible for the promotion of surface catalytic activity. The generation of dislocations and related surface modification by strong nonlinear SAWs are currently investigated for metal substrates described by realistic embedded atom method interatomic potentials, with the results to be reported elsewhere.

## VIII. SUMMARY

A computational “synchrotron” approach for atomistic MD modeling of long-term evolution of a surface acoustic wave profile during the free propagation of the wave in a nonlinear medium is developed in this work and applied for the initial exploration of nonlinear effects in the propagation and dissipation of SAWs. The computational approach takes advantage of the periodic boundary conditions applied in the direction of the wave propagation, which enables the reuse of the same material for the simulation of the nonlinear evolution of the wave profile without the need for extending the size of the computational system to match the distance covered by the wave. The circulation of the wave can be continued up to the onset of the irreversible material modification produced by the shock wave formed due to the nonlinear sharpening of the wave profile. The constraints on the choice of computational parameters suitable for the simulation of free propagation of SAWs are evaluated based on the consideration of the wavelength-dependent depth, the elastic field produced by SAWs in the bulk of the target, and the strength of the dynamic coupling between the top and bottom surfaces of the computational system.

First applications of the computational model demonstrate the ability of atomistic simulations to reproduce the key characteristics of the nonlinear transformations of SAW profiles predicted in earlier continuum-level numerical simulations and theoretical investigations. In particular, the MD simulations reproduce the generation of higher frequency harmonics and the formation of a sharp shock front during the propagation of an initially sinusoidal wave in a nonlinear medium, the increasing localization of the acoustic energy near the surface of the substrate due to the frequency up-conversion, and the rapid wave dissipation at the shock front. The MD simulations also predict the emergence of characteristic cusps in the horizontal strain and velocity profiles, as well as sharp spikes in the vertical velocity profiles, with the peak values of the velocity/strain significantly exceeding those of the initial wave before the shock formation. Notably, the complex nonlinear evolution of the wave profiles and surface energy localization are reproduced in MD simulations without any assumptions on the mechanisms and rates of energy redistribution among the harmonics, and all the nonlinear effects, including the sign of nonlinearity parameter, naturally emerge from atomic-scale dynamics defined by the interatomic potential and crystal structure of the substrate. As a result, the MD simulations enable verification of theoretical assumptions involved in the description of the harmonic generation and wave dissipation, particularly at large strains, where higher order nonlinear effects may play an important role.

More importantly, MD simulations of SAWs make it possible to explore the effects that are beyond the capabilities of the continuum-level models. Several examples of such effects revealed in the first MD simulations are listed below.

*An unusual quadratic in k correction to the dispersion relation at short wavelengths*, where the discrete atomic structure of the material starts to play a role, has been observed in simulations of relatively weak SAWs. This observation is attributed to the frequency-dependent localization of SAWs near the surface of the substrate, which leads to the decrease in the effective number of discrete scattering elements (atoms) experienced by SAWs of shorter wavelength.*A new mechanism of the energy dissipation at the SAW shock front*is observed to activate at the final stage of the sharpening of the wave profile, when the nonlinear transformation of a SAW wave profile naturally generates a cascade of extremely high-frequency harmonics. The energy channeling from the highest frequency harmonics of the SAW to the emission of a cloud of phonon wave packets slowly descending into the bulk of the substrate is clearly visualized by snapshots from MD simulations. Further investigations are needed to quantify this dissipation mechanism and formulate it in a form suitable for incorporation into continuum-level models.The dependence of the

*shock front formation distance*on the initial amplitude and wavelength of a SAW predicted in MD simulations can only be described by the commonly accepted theoretical scaling law at a semiquantitative level. The deviations between the theory and MD results suggest the need for refinement of the theoretical description, particularly at short SAW wavelengths and large strain amplitudes. By enabling controlled variation of the nonlinear material properties and the parameters of SAWs, the MD simulations can play an important role in the advancement of the shock formation theory.*The emission of dislocations*at a substantial distance from the source of an initially sinusoidal SAW is observed in MD simulations and explained by the nonlinear sharpening and cusping of the wave profile, leading to the generation of surface strain significantly exceeding the initial strain amplitude of the wave. These results demonstrate the ability of MD simulations to provide insights into the mechanisms of SAW-induced localized plastic deformation and suggest future applications in the exploration of various acoustic effects involving material modification and damage.

Overall, the computational methodology developed in this work for atomistic MD modeling of free nonlinear SAW propagation and dissipation provides a unified description of weakly nonlinear transformations and ultimate manifestations of the material nonlinearity in the form of material damage, thus enabling detailed analysis of the intrinsic properties of nonlinear SAWs. Moreover, the MD simulation opens up a broad range of opportunities for investigation of material modification by SAWs, acoustically induced surface processes, such as desorption, diffusion, and chemical reactions, as well as the interaction of SAWs with pre-existing crystal defects and other material heterogeneities. At a more general level, the atomistic modeling of SAWs may provide a solid scientific footing for the development of new material/surface processing techniques that fully utilize the benefits of the acoustic stimulus.

## SUPPLEMENTARY MATERIAL

The supplementary material includes three mp4 movies showing the evolution of the horizontal surface strain $exx$ profiles predicted in MD simulations illustrated in Figs. 11(b), 11(c), 15, and 16. All simulations are performed at a substrate temperature of 0.15ɛ/*k _{B}*.

**Movie #1** (937sigma-3.4%.mp4): $exx$ profile for a SAW with $exx0=3.4%$ and $\lambda =937\sigma $. The movie covers the total time of 8936τ from the start of the simulation and illustrates a gradual sharpening of the wave, shock front formation, cusping of the wave profile on the compressive side of the shock front, and attenuation of the wave due to the energy dissipation at the shock front.

**Movie #2** (937sigma-4.7%.mp4): $exx$ profile for a SAW with $exx0=4.7%$ and $\lambda =937\sigma $. The movie covers the total time of 1638τ from the start of the simulation and illustrates a gradual sharpening of the wave leading to the shock front formation and cusping of the wave profile, followed by several cycles of abrupt drop in the maximum strain upon the emission of dislocations and a partial recovery of the maximum strain through further sharpening of the wave and formation of a new cusp. The first dislocations are emitted at a time of 1002τ.

**Movie #3** (1874sigma-4.7%.mp4): $exx$ profile for a SAW with $exx0=4.7%$ and $\lambda =1874\sigma $. The movie covers the total time of 1926τ from the start of the simulation and illustrates the same processes as Movie #2, but for a wave with twice larger wavelength.

## ACKNOWLEDGMENTS

Financial support for this work was provided by the National Science Foundation (NSF) through Grant No. CMMI-1562929. Computational support was provided by the NSF through the Extreme Science and Engineering Discovery Environment (Project No. TGDMR110090). V.Y.Z. also acknowledges the support from the Russian Foundation for Basic Research (RFBR) under Grant No. 19-05-00536.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

### APPENDIX: ANALYTICAL SOLUTION FOR HARMONIC RAYLEIGH WAVE IN A CRYSTAL WITH CUBIC SYMMETRY

The analytical solution for displacements and velocities produced by a harmonic Rayleigh wave propagating in a crystal with cubic symmetry can be written based on the analysis provided in Ref. 84. Assuming that the *x* axis is directed along the wave propagation, the *z* axis is normal to the surface, and *z* = 0 is the position of the unperturbed surface, the equations for a wave propagating along the [100] direction of a (001) surface of a cubic crystal can be written in the following form:

where *u* and *w* are the displacements in the *x* and *z* directions, respectively, *A*, *k*, and *c* are the surface wave amplitude, wavenumber, and velocity, respectively. The parameters $q1$ and $q2$ are the two of the four roots of the following equation selected to produce a solution exponentially decaying with depth under the surface of the substrate:

where $\rho $ is the material density and $C11$, $C12$, and $C44$ are the elastic constants of the material. The parameters $\gamma i$ in Eqs. (A1) and (A2) can be expressed through the corresponding parameters $qi$ as

It should be noted that although Eqs. (A1) and (A2) are similar to the corresponding expressions formulated for Rayleigh waves in an isotropic material, in the case of the cubic crystals, the parameters $q1$ and $q2$ are complex-valued. These parameters define the in-depth variation of the SAW amplitude, and the fact that parameters $q1$ and $q2$ are complex numbers implies that the wave amplitude not only exponentially decays with depth but also experiences oscillations [see Fig. 6(b)].

Finally, the wave velocity *c* in Eqs. (A1) and (A2) can be determined in terms of the auxiliary parameter $R=\rho c2/C11$ by solving the following equation:

To generate a running SAW, both the displacements and velocities have to be added to the instantaneous atomic displacements and velocities arising from the thermal vibrations in an MD system equilibrated at a desired temperature. The velocity components can be calculated from Eqs. (A1) and (A2) as time derivatives of *u* and *w*.