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.

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 crust1 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 properties3–7 and surface defects,8,9 signal-processing10–12 and acousto-optical13 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 formalism27–32,44 applied either in the frequency domain using a system of equations for interacting harmonics27 or in the form of an evolution equation for the SAW profile in time domain29,30 have been developed to describe the variation of different harmonic constituents and the overall evolution of SAW profiles in isotropic28–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 system56–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 regime10) 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.

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” technique70 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 

FIG. 1.

Schematic representation of the interactions between different components of the proposed computational “synchrotron” approach to simulation of the generation of SAWs, their nonlinear propagation, and structural modification of the surface region (see description in the text).

FIG. 1.

Schematic representation of the interactions between different components of the proposed computational “synchrotron” approach to simulation of the generation of SAWs, their nonlinear propagation, and structural modification of the surface region (see description in the text).

Close modal

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 conditions75–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 interactions76–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.

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 nonlinearity80 (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 function81 used to ensure that the interactions vanish at a cutoff distance of rc = 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.

FIG. 2.

The stress–strain curves calculated for an fcc crystal described by the Lennard–Jones (green curve) and harmonic (red curve) potentials and loaded under the uniaxial strain conditions. The dashed line shows a linear stress–strain dependence.

FIG. 2.

The stress–strain curves calculated for an fcc crystal described by the Lennard–Jones (green curve) and harmonic (red curve) potentials and loaded under the uniaxial strain conditions. The dashed line shows a linear stress–strain dependence.

Close modal

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 τ=(mσ2/ε)1/2, υ=σ/τ=(ε/m)1/2, and T = ɛ/kB, respectively, where kB 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.

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, Lx, 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., Ly> 2rc. To minimize computational time, Ly is chosen to be close to the minimum value, Ly = 7.8σ, which corresponds to 5 fcc unit cells.

FIG. 3.

Schematic representation of the computational system used in simulations of free propagation of SAWs. The color shows an instantaneous distribution of kinetic energy in the field of an initial sinusoidal SAW propagating in the direction of the x axis. The blue and red colors correspond to low and high values of the kinetic energy, respectively.

FIG. 3.

Schematic representation of the computational system used in simulations of free propagation of SAWs. The color shows an instantaneous distribution of kinetic energy in the field of an initial sinusoidal SAW propagating in the direction of the x axis. The blue and red colors correspond to low and high values of the kinetic energy, respectively.

Close modal

The choice of the depth of the computational system, Lz, 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 at82,83 or very close to47 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 Lz 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−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 Lz = 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 Lz 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 Lz = 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.

FIG. 4.

Contour plots showing the evolution of the kinetic energy in MD simulations of SAWs propagating in two systems of the same depth of Lz = 2λ but different boundary conditions applied at the bottom of the computational system, free boundary in (a) and rigid boundary in (b). The horizontal axis shows the distance passed by the wave, which is proportional to the time of the simulation. The vertical axis represents the depth under the surface where the initial wave is generated at the start of the simulation. The color corresponds to the kinetic energy per atom averaged over one SAW wave period and surface subdomains with a size of 3.8σ in the z direction (axes are defined in Fig. 3). The simulations are performed for SAWs with wavelength of λ = 155σ, strain amplitude of e0 = 8 × 10−5, and zero temperature of the substrate. Periodic energy exchange between the dynamically coupled SAWs propagating along the two free surfaces of the system is observed in panel (a), where the direction of the energy transfer is shown by white dashed arrows. No such energy redistribution is seen in panel (b), where the rigid boundary condition is applied at the bottom.

FIG. 4.

Contour plots showing the evolution of the kinetic energy in MD simulations of SAWs propagating in two systems of the same depth of Lz = 2λ but different boundary conditions applied at the bottom of the computational system, free boundary in (a) and rigid boundary in (b). The horizontal axis shows the distance passed by the wave, which is proportional to the time of the simulation. The vertical axis represents the depth under the surface where the initial wave is generated at the start of the simulation. The color corresponds to the kinetic energy per atom averaged over one SAW wave period and surface subdomains with a size of 3.8σ in the z direction (axes are defined in Fig. 3). The simulations are performed for SAWs with wavelength of λ = 155σ, strain amplitude of e0 = 8 × 10−5, and zero temperature of the substrate. Periodic energy exchange between the dynamically coupled SAWs propagating along the two free surfaces of the system is observed in panel (a), where the direction of the energy transfer is shown by white dashed arrows. No such energy redistribution is seen in panel (b), where the rigid boundary condition is applied at the bottom.

Close modal
FIG. 5.

The dependences of the amplitudes of the SAWs induced near the free bottom surfaces of computational systems with depths, Lz, equal to 2λ, 3λ, and 4λ. In all three simulations, the initial SAWs generated near the top surfaces of the computational systems have a wavelength of λ = 155σ and a strain amplitude of e0 = 8 × 10−5. The initial temperature of the systems is zero. The curve for Lz = 2λ corresponds to panel (a) in Fig. 4. These examples illustrate the reduction in the rate of energy exchange between the SAWs propagating along the top and bottom surfaces of the substrate with increasing Lz. They also show that even with the relatively weak coupling between the two waves at Lz = 4λ, the cumulative energy exchange may be rather efficient.

FIG. 5.

The dependences of the amplitudes of the SAWs induced near the free bottom surfaces of computational systems with depths, Lz, equal to 2λ, 3λ, and 4λ. In all three simulations, the initial SAWs generated near the top surfaces of the computational systems have a wavelength of λ = 155σ and a strain amplitude of e0 = 8 × 10−5. The initial temperature of the systems is zero. The curve for Lz = 2λ corresponds to panel (a) in Fig. 4. These examples illustrate the reduction in the rate of energy exchange between the SAWs propagating along the top and bottom surfaces of the substrate with increasing Lz. They also show that even with the relatively weak coupling between the two waves at Lz = 4λ, the cumulative energy exchange may be rather efficient.

Close modal

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, Lz, 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.

FIG. 6.

(a) The evolution of the amplitude of the second harmonic of a SAW propagating in MD simulations performed for computational systems with depths, Lz, equal to 2λ (red dashed line), 3λ (green solid line), and 4λ (blue dash-dotted line) and rigid boundary conditions applied at the bottom of all systems. The wavelength of the SAWs generated at the start of the simulations is 155σ, and a strain amplitude is e0 = 2 × 10−5. The initial temperature of the systems is zero. The plots demonstrate a significant difference between the cases of Lz = 2λ and Lz = 3λ, while the evolution of the second harmonic is almost the same for Lz = 3λ and Lz = 4λ. (b) The Rayleigh wave solution for the normalized displacement amplitudes vs depth under the surface (normalized by the wavelength λ) at the source of the wave, x = 0.

FIG. 6.

(a) The evolution of the amplitude of the second harmonic of a SAW propagating in MD simulations performed for computational systems with depths, Lz, equal to 2λ (red dashed line), 3λ (green solid line), and 4λ (blue dash-dotted line) and rigid boundary conditions applied at the bottom of all systems. The wavelength of the SAWs generated at the start of the simulations is 155σ, and a strain amplitude is e0 = 2 × 10−5. The initial temperature of the systems is zero. The plots demonstrate a significant difference between the cases of Lz = 2λ and Lz = 3λ, while the evolution of the second harmonic is almost the same for Lz = 3λ and Lz = 4λ. (b) The Rayleigh wave solution for the normalized displacement amplitudes vs depth under the surface (normalized by the wavelength λ) at the source of the wave, x = 0.

Close modal

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 Lx × Ly × Lz = 1860σ × 7.8σ × 5580σ (1200 × 5 × 3600 fcc unit cells) and consists of 86.4 × 106 atoms. It should be noted that the conclusion on the minimum depth requirement of Lz ≥ 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 x = 0 is shown in Fig. 6(b). The minimum depth Lz 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−2A, 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, Lx, should be a multiple of the wavelength of the fundamental harmonic of the generated wave, i.e., Lx = 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.

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 L12b defined by the difference between the wavevector of the fundamental harmonic, k1, and half the wavevector of the second harmonic, k2,

L12b=2π/|2k1k2|.
(1)

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 L12b/2 acquires a relative phase mismatch equal to π due to the difference between the phase velocity of the fundamental harmonic, c1=ω1/k1, and that of the second harmonic, c2=2ω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.

FIG. 7.

The evolution of the amplitudes of second harmonics of initially sinusoidal SAWs shown as a function of distance of free wave propagation, as predicted in MD simulations (a). The SAWs are generated with the same initial displacement amplitude of ∼0.019σ and wavelengths λ of 155σ, 77.5σ, and 38.8σ (the corresponding curves are marked on the plot). The initial temperature of the system is equal to zero. The spatial period of the beating in the amplitude of the second harmonic, L12b, decreases with decreasing λ due to the faster accumulation of the phase mismatch related to the spatial dispersion, as shown in (b), where the results for bulk longitudinal (compressive) and shear waves are also shown for comparison. The data presented in (b) are provided in Table I.

FIG. 7.

The evolution of the amplitudes of second harmonics of initially sinusoidal SAWs shown as a function of distance of free wave propagation, as predicted in MD simulations (a). The SAWs are generated with the same initial displacement amplitude of ∼0.019σ and wavelengths λ of 155σ, 77.5σ, and 38.8σ (the corresponding curves are marked on the plot). The initial temperature of the system is equal to zero. The spatial period of the beating in the amplitude of the second harmonic, L12b, decreases with decreasing λ due to the faster accumulation of the phase mismatch related to the spatial dispersion, as shown in (b), where the results for bulk longitudinal (compressive) and shear waves are also shown for comparison. The data presented in (b) are provided in Table I.

Close modal
TABLE I.

The distance L12b needed for the double-frequency wave with wavevector k2 to accumulate the phase mismatch |2k1k2|L12b=2π with the fundamental harmonic for different wavevectors k1, as determined in MD simulations of SAWs, bulk longitudinal (compressive) waves, and bulk shear waves of different wavelength λ. The corresponding plots of L12b(λ) are provided in Fig. 7(b).

Wavelength λ [σ] (unit cells) 19 (12) 39 (25) 78 (50) 155 (100) 
 Surface waves 620 3440 15 000 61 800 
Distance L12b [σ] 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 L12b [σ] 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 ω=ck is cubic in k, namely, ωckαk3, where normally α 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 ω(k) for the initial deviation from the non-dispersive propagation [which is equivalent to the cubic in ωck correction to the wavevector k(ω)], the spatial frequency 2π/L12b=|2k1k2| 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, L12bλ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 λ=77.5σ to λ=155σ leads to an increase in the beating period L12b 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 λk1. 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 k3 dispersion correction due to the scattering by periodically spaced “screens” should be multiplied by a factor of k1, reflecting the decrease in the effective number of scattering elements with increasing k. Consequently, for SAWs, the dispersion correction becomes quadratic, k1k3=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 k1. As a result, for such extremely short SAWs, the initial k2 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σ.

FIG. 8.

Schematic illustration of the origin of the k2 correction to the dispersion relation for SAWs. The elastic fields generated by SAWs of three different wavelength (black curves) extend into the bulk of the substrate to the depth on the order of λk1, leading to the corresponding scaling of size of the periodically spaced vertical “screens” (red dashed lines) composed of discrete scattering elements (atoms). As a result, the conventional k3 dispersion correction typical for bulk waves or 1D waves in chains of masses and springs transforms into quadratic dispersion correction, k1k3=k2, for SAWs.

FIG. 8.

Schematic illustration of the origin of the k2 correction to the dispersion relation for SAWs. The elastic fields generated by SAWs of three different wavelength (black curves) extend into the bulk of the substrate to the depth on the order of λk1, leading to the corresponding scaling of size of the periodically spaced vertical “screens” (red dashed lines) composed of discrete scattering elements (atoms). As a result, the conventional k3 dispersion correction typical for bulk waves or 1D waves in chains of masses and springs transforms into quadratic dispersion correction, k1k3=k2, for SAWs.

Close modal

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. VVII, 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.

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 harmonics25 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 

FIG. 9.

Contour plots showing the evolution of the kinetic energy distribution in a surface region of a substrate during the propagation of a SAW with an initial strain amplitude of 2% and the fundamental wavelength of 930σ, as predicted in an MD simulation. The initial temperature of the substrate is equal to zero. Color in the contour plot corresponds to the kinetic energy per atom averaged over the wave period, and the energy distribution is shown as a function of depth under the surface and distance traveled by the wave. The contour plot in (a) shows the gradual localization of the energy in the surface region during the nonlinear SAW propagation. Panel (b) shows the variation of the kinetic energy with a lower energy scale in a deeper part of the system, below 500σ, for the wave propagation distance exceeding 10 000σ. A part of the region shown in (b) is outlined by the dashed rectangle in (a). After traveling a distance of 10 000σ, the SAW develops a sharp shock front, leading to the onset of intense dissipation of the wave energy through an effective conversion of the high-frequency harmonics to thermal phonons. The propagation of the phonon wave packets carrying the energy into the bulk of the substrate is clearly visible in the right-hand part of panel (b).

FIG. 9.

Contour plots showing the evolution of the kinetic energy distribution in a surface region of a substrate during the propagation of a SAW with an initial strain amplitude of 2% and the fundamental wavelength of 930σ, as predicted in an MD simulation. The initial temperature of the substrate is equal to zero. Color in the contour plot corresponds to the kinetic energy per atom averaged over the wave period, and the energy distribution is shown as a function of depth under the surface and distance traveled by the wave. The contour plot in (a) shows the gradual localization of the energy in the surface region during the nonlinear SAW propagation. Panel (b) shows the variation of the kinetic energy with a lower energy scale in a deeper part of the system, below 500σ, for the wave propagation distance exceeding 10 000σ. A part of the region shown in (b) is outlined by the dashed rectangle in (a). After traveling a distance of 10 000σ, the SAW develops a sharp shock front, leading to the onset of intense dissipation of the wave energy through an effective conversion of the high-frequency harmonics to thermal phonons. The propagation of the phonon wave packets carrying the energy into the bulk of the substrate is clearly visible in the right-hand part of panel (b).

Close modal

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×103104σ, 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σ.

FIG. 10.

The sequence of snapshots of instantaneous distributions of kinetic energy in the substrate at the moments corresponding to the propagation distances of 9000σ, 11 000σ, and 13 000σ shown for a simulation illustrated by Fig. 9. The corresponding profiles of the surface strain exx are shown on the top of the contour plots. Wave packets corresponding to thermal phonons with the wavelengths of 2–4 fcc unit cells are clearly visible in the enlarged view provided as an inset in the right panel. The dissipation of the SAW through the emission of a cloud of phonons descending into the bulk of the substrate can be readily seen in the sequence of snapshots.

FIG. 10.

The sequence of snapshots of instantaneous distributions of kinetic energy in the substrate at the moments corresponding to the propagation distances of 9000σ, 11 000σ, and 13 000σ shown for a simulation illustrated by Fig. 9. The corresponding profiles of the surface strain exx are shown on the top of the contour plots. Wave packets corresponding to thermal phonons with the wavelengths of 2–4 fcc unit cells are clearly visible in the enlarged view provided as an inset in the right panel. The dissipation of the SAW through the emission of a cloud of phonons descending into the bulk of the substrate can be readily seen in the sequence of snapshots.

Close modal

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 conditions75–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

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ɛ/kB, 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.

FIG. 11.

Evolution of the horizontal surface strain exx profiles predicted in MD simulations of free propagation of SAWs with initial strain amplitudes of 2% (a), 3.4% (b), and 4.7% (c). The SAWs propagate in the direction from left to right, the fundamental wavelength of the waves is 937σ, the temperature of the substrate is 0.15ɛ/kB, and the velocity of the wave propagation at this temperature is about 4.87 σ/τ. The first panel in each row shows the initial sinusoidal wave. The moments of time for which the strain profiles are shown in the middle and right panels are marked on the plots and are chosen to represent the same stage of the nonlinear distortion of the SAW profiles in each simulation. Gradual steepening of the SAW profiles leading to the formation of shock fronts is observed in all simulations and is shown in the middle panels. Pronounced asymmetric increase in the strain magnitude at larger propagation distances (larger evolution times) is observed for sufficiently strong waves, as can be seen in right panels of (b) and (c). The shock fronts form at ascending sides of the SAW profiles, as expected for the positive value of the parameter of local nonlinearity intrinsic to the LJ potential. The red dashed–dotted line in the right panel of (b) shows the wave profile at a much longer time of 8936τ and illustrates wave attenuation due to the strong dissipation near the shock front. The movies showing the evolution of the strain profiles in simulations illustrated in (b) and (c) are provided in the supplementary material.

FIG. 11.

Evolution of the horizontal surface strain exx profiles predicted in MD simulations of free propagation of SAWs with initial strain amplitudes of 2% (a), 3.4% (b), and 4.7% (c). The SAWs propagate in the direction from left to right, the fundamental wavelength of the waves is 937σ, the temperature of the substrate is 0.15ɛ/kB, and the velocity of the wave propagation at this temperature is about 4.87 σ/τ. The first panel in each row shows the initial sinusoidal wave. The moments of time for which the strain profiles are shown in the middle and right panels are marked on the plots and are chosen to represent the same stage of the nonlinear distortion of the SAW profiles in each simulation. Gradual steepening of the SAW profiles leading to the formation of shock fronts is observed in all simulations and is shown in the middle panels. Pronounced asymmetric increase in the strain magnitude at larger propagation distances (larger evolution times) is observed for sufficiently strong waves, as can be seen in right panels of (b) and (c). The shock fronts form at ascending sides of the SAW profiles, as expected for the positive value of the parameter of local nonlinearity intrinsic to the LJ potential. The red dashed–dotted line in the right panel of (b) shows the wave profile at a much longer time of 8936τ and illustrates wave attenuation due to the strong dissipation near the shock front. The movies showing the evolution of the strain profiles in simulations illustrated in (b) and (c) are provided in the supplementary material.

Close modal

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ɛ/kB.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 nonlinearity80 (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.

FIG. 12.

The horizontal surface strain exx profile predicted in an MD simulation of free propagation of a SAW in a substrate described by the short-range harmonic (parabolic) pair potential and exhibiting a negative value of the parameter of local nonlinearity. The SAW propagates in the direction from left to right, the fundamental wavelength of the wave is 155σ, the initial strain is 8.5%, and the initial temperature of the substrate is zero. The strain profile is shown for a time of 197τ, when the shock front of the wave is well developed. A pronounced asymmetric increase in the strain magnitude and the formation of a cusp exceeding the level of the initial strain are clearly visible. In contrast to the simulations performed with the LJ potential (see Fig. 11), however, the cusp appears at the upper (tensile) part of the shock front, and the shock itself forms at the descending side in the SAW profile.

FIG. 12.

The horizontal surface strain exx profile predicted in an MD simulation of free propagation of a SAW in a substrate described by the short-range harmonic (parabolic) pair potential and exhibiting a negative value of the parameter of local nonlinearity. The SAW propagates in the direction from left to right, the fundamental wavelength of the wave is 155σ, the initial strain is 8.5%, and the initial temperature of the substrate is zero. The strain profile is shown for a time of 197τ, when the shock front of the wave is well developed. A pronounced asymmetric increase in the strain magnitude and the formation of a cusp exceeding the level of the initial strain are clearly visible. In contrast to the simulations performed with the LJ potential (see Fig. 11), however, the cusp appears at the upper (tensile) part of the shock front, and the shock itself forms at the descending side in the SAW profile.

Close modal

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 υ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 υ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.

FIG. 13.

The horizontal υx and vertical υz velocity profiles predicted in MD simulations of free propagation of SAWs in a substrate described by the LJ potential. The SAWs propagate in the direction from left to right, the fundamental wavelength of the waves is 937σ, the temperature of the substrate is 0.15ɛ/kB, and the initial strain amplitude is 2% in (a) and 4.7% in (b). The red dashed lines show the initial velocity profiles at the start of the simulations. The black lines are the velocity profiles plotted for times when the shock fronts of the waves are well developed, 2730τ in (a) and 945τ in (b). The profiles of the horizontal velocities exhibit shock fronts similar to those observed for the horizontal strain, i.e., υxexx (see right panels in Fig. 11), while the profiles of the vertical velocities exhibit sharp spikes.

FIG. 13.

The horizontal υx and vertical υz velocity profiles predicted in MD simulations of free propagation of SAWs in a substrate described by the LJ potential. The SAWs propagate in the direction from left to right, the fundamental wavelength of the waves is 937σ, the temperature of the substrate is 0.15ɛ/kB, and the initial strain amplitude is 2% in (a) and 4.7% in (b). The red dashed lines show the initial velocity profiles at the start of the simulations. The black lines are the velocity profiles plotted for times when the shock fronts of the waves are well developed, 2730τ in (a) and 945τ in (b). The profiles of the horizontal velocities exhibit shock fronts similar to those observed for the horizontal strain, i.e., υxexx (see right panels in Fig. 11), while the profiles of the vertical velocities exhibit sharp spikes.

Close modal

Taking into account the proportionality υxexx, 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 substrates31 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 υ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 suggested94 that it has the same functional dependence as the distance of the shock formation for the initially sinusoidal bulk waves,95,96

xsh=1/(βexx0k),
(2)

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

The shock front formation distance xsh predicted in MD simulations for an LJ substrate with a temperature of 0.15ɛ/kB is plotted in Fig. 14 as a function of the inverse strain 1/exx0 and wavelength λ=2πk1. The error bars in Fig. 14 show the uncertainty of the measured distances xsh 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 xsh 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,94 cannot be applied in the analysis of the MD results.

FIG. 14.

The dependence of the shock formation distance xsh on the initial strain amplitude exx0 (a) and fundamental wavelength λ (b) predicted in MD simulations of SAW propagating in a substrate described by the LJ potential. The temperature of the substrate is 0.15ɛ/kB. The dependence on the initial strain amplitude in (a) is evaluated for wavelength λ=937σ. The dashed lines show the prediction of Eq. (2) with β = 0.85. The insets in (a) show the horizontal strain profiles after the shock front formation that correspond to two data points. The dependence on wavelength in (b) is evaluated for two values of the initial stain amplitude, exx0=3.4% (red points) and exx0=4.7% (green points). The error bars show the uncertainty in the measured distances xsh introduced by the variability in the slopes of fully developed shock fronts.

FIG. 14.

The dependence of the shock formation distance xsh on the initial strain amplitude exx0 (a) and fundamental wavelength λ (b) predicted in MD simulations of SAW propagating in a substrate described by the LJ potential. The temperature of the substrate is 0.15ɛ/kB. The dependence on the initial strain amplitude in (a) is evaluated for wavelength λ=937σ. The dashed lines show the prediction of Eq. (2) with β = 0.85. The insets in (a) show the horizontal strain profiles after the shock front formation that correspond to two data points. The dependence on wavelength in (b) is evaluated for two values of the initial stain amplitude, exx0=3.4% (red points) and exx0=4.7% (green points). The error bars show the uncertainty in the measured distances xsh introduced by the variability in the slopes of fully developed shock fronts.

Close modal

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 β 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 β 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 xsh 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 (λ=26σandλ=39σ), 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 β at large strains approaching the material damage threshold.

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ɛ/kB. 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).

FIG. 15.

Snapshots of the atomic structure in the top region of the substrate and the corresponding strain profiles predicted in an MD simulation of free propagation of a SAW with an initial strain amplitude of 4.7% and a fundamental wavelength of 1874σ (1200 fcc unit cells). The temperature of the substrate is 0.15ɛ/kB. The first snapshot shown in (a) is for a wave that passed a distance of 7670σ from the source of the initial sinusoidal wave. The two subsequent snapshots in (b) and (c) are for two moments of time when the wave passes extra 510σ and 1020σ, respectively. The snapshots depict only the top ∼200σ deep parts of the computational system, and the corresponding profiles of the SAW strain are shown below the snapshots. The positions of the shock front are marked by black arrows. Sharpening of the wave leads to the generation of material damage in the form of localized zones of plastic deformation. The darker (blue online) color corresponds to intact fcc structure of the substrate atoms, and lighter dashed zones correspond to the regions affected by the plastic deformation (stacking faults left behind by Shockley partial dislocations are shown by green color).

FIG. 15.

Snapshots of the atomic structure in the top region of the substrate and the corresponding strain profiles predicted in an MD simulation of free propagation of a SAW with an initial strain amplitude of 4.7% and a fundamental wavelength of 1874σ (1200 fcc unit cells). The temperature of the substrate is 0.15ɛ/kB. The first snapshot shown in (a) is for a wave that passed a distance of 7670σ from the source of the initial sinusoidal wave. The two subsequent snapshots in (b) and (c) are for two moments of time when the wave passes extra 510σ and 1020σ, respectively. The snapshots depict only the top ∼200σ deep parts of the computational system, and the corresponding profiles of the SAW strain are shown below the snapshots. The positions of the shock front are marked by black arrows. Sharpening of the wave leads to the generation of material damage in the form of localized zones of plastic deformation. The darker (blue online) color corresponds to intact fcc structure of the substrate atoms, and lighter dashed zones correspond to the regions affected by the plastic deformation (stacking faults left behind by Shockley partial dislocations are shown by green color).

Close modal

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 Si5–7 and fused silica.7 

FIG. 16.

The abrupt change of the strain profile at the onset of plastic deformation in a simulation illustrated in Fig. 15. The strain profiles just before and just after the SAW-induced emission of dislocations are shown by the red dashed and black solid lines, respectively. The profiles are plotted in a frame of reference moving with the wave. The inset with a light blue background provides an enlarged view of a part of the strain profile near the sharp cusp generated at the SAW shock front. A movie showing the evolution of the strain profile and illustrating the effect of the onset of material modification on the shape of the shock front is provided in the supplementary material.

FIG. 16.

The abrupt change of the strain profile at the onset of plastic deformation in a simulation illustrated in Fig. 15. The strain profiles just before and just after the SAW-induced emission of dislocations are shown by the red dashed and black solid lines, respectively. The profiles are plotted in a frame of reference moving with the wave. The inset with a light blue background provides an enlarged view of a part of the strain profile near the sharp cusp generated at the SAW shock front. A movie showing the evolution of the strain profile and illustrating the effect of the onset of material modification on the shape of the shock front is provided in the supplementary material.

Close modal

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 studies51–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.

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.

  1. 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.

  2. 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.

  3. 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.

  4. 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.

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ɛ/kB.

Movie #1 (937sigma-3.4%.mp4): exx profile for a SAW with exx0=3.4% and λ=937σ. 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 λ=937σ. 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 λ=1874σ. 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.

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.

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

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:

u(x,z,t)=AReexp[q1kz]q1γ1q2γ2exp[q2kz]exp[ik(xct)],
(A1)
w(x,z,t)=AImγ1exp[q1kz]q1γ1q2γ2γ2exp[q2kz]exp[ik(xct)],
(A2)

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:

(C11ρc2C44q2)(C44ρc2C11q2)+q2(C44+C11)2=0,
(A3)

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

γi=qiC44+C12C44C11RC11qi2.
(A4)

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=ρc2/C11 by solving the following equation:

1C11C44R1C122C112R2=R2(1R).
(A5)

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.

1.
A.
Ben-Menahem
and
S. J.
Singh
,
Seismic Waves and Sources
(
Springer-Verlag
,
New York
,
1981
).
2.
C. E.
O'Connell-Rodwell
,
B. T.
Arnason
, and
L. A.
Hart
, “
Seismic properties of Asian elephant (Elephas maximus) vocalizations and locomotion
,”
J. Acoust. Soc. Am.
108
,
3066
3072
(
2000
).
3.
Modern Acoustical Techniques for the Measurement of Mechanical Properties
, edited by
M.
Levy
,
H. E.
Bass
, and
R.
Stern
(
Academic Press
,
New York
,
2001
).
4.
P.
Hess
, “
Surface acoustic waves in materials science
,”
Phys. Today
55
(
3
),
42
47
(
2002
).
5.
A. M.
Lomonosov
and
P.
Hess
, “
Impulsive fracture of silicon by elastic surface pulses with shocks
,”
Phys. Rev. Lett.
89
,
095501
(
2002
).
6.
V. V.
Kozhushko
,
A. M.
Lomonosov
, and
P.
Hess
, “
Intrinsic strength of silicon crystals in pure- and combined-mode fracture without precrack
,”
Phys. Rev. Lett.
98
,
195505
(
2007
).
7.
G.
Lehmann
,
A. M.
Lomonosov
,
P.
Hess
, and
P.
Gumbsch
, “
Impulsive fracture of fused quartz and silicon crystals by nonlinear surface acoustic waves
,”
J. Appl. Phys.
94
,
2907
2914
(
2003
).
8.
A. M.
Lomonosov
,
P. V.
Grigoriev
, and
P.
Hess
, “
Sizing of partially closed surface-breaking microcracks with broadband Rayleigh waves
,”
J. Appl. Phys.
105
,
084906
(
2009
).
9.
F.
Hofmann
,
M. P.
Short
, and
C. A.
Dennett
, “
Transient grating spectroscopy: An ultrarapid, nondestructive materials evaluation technique
,”
MRS Bull.
44
,
392
402
(
2019
).
10.
I. A.
Viktorov
,
Rayleigh and Lamb Waves Physical Theory and Applications
(
Plenum Press
,
New York
,
1967
).
11.
M. F.
Lewis
, “
Rayleigh waves—A progress report
,”
Eur. J. Phys.
16
,
1
7
(
1995
).
12.
K.-Y.
Hashimoto
,
Surface Acoustic Wave Devices in Telecommunications: Modeling and Simulation
(
Springer-Verlag
,
Berlin
,
2000
).
13.
Y. J.
Liu
,
X.
Ding
,
S.-C. S.
Lin
,
J.
Shi
,
I.-K.
Chiang
, and
T. J.
Huang
, “
Surface acoustic wave driven light shutters using polymer-dispersed liquid crystals
,”
Adv. Mater.
23
,
1656
1659
(
2011
).
14.
G.
Schmera
and
L. B.
Kish
, “
Surface diffusion enhanced chemical sensing by surface acoustic waves
,”
Sens. Actuators B Chem.
93
,
159
163
(
2003
).
15.
A.
Buvailo
,
Y.
Xing
,
J.
Hines
, and
E.
Borguet
, “
Thin polymer film based rapid surface acoustic wave humidity sensors
,”
Sens. Actuators B Chem.
156
,
444
449
(
2011
).
16.
L. Y.
Yeo
and
J. R.
Friend
, “
Surface acoustic wave microfluidics
,”
Annu. Rev. Fluid Mech.
46
,
379
406
(
2014
).
17.
X.
Ding
,
S.-C. S.
Lin
,
B.
Kiraly
,
H.
Yue
,
S.
Li
,
I.-K.
Chiang
,
J.
Shi
,
S. J.
Benkovic
, and
T. J.
Huang
, “
On-chip manipulation of single microparticles, cells, and organisms using surface acoustic waves
,”
Proc. Natl. Acad. Sci. U.S.A.
109
,
11105
11109
(
2012
).
18.
M. C.
Jo
and
R.
Guldiken
, “
Dual surface acoustic wave-based active mixing in a microfluidic channel
,”
Sens. Actuators A Phys.
196
,
1
7
(
2013
).
19.
C. D.
Wood
,
S. D.
Evans
,
J. E.
Cunningham
,
R.
O’Rorke
,
C.
Wälti
, and
A. G.
Davies
, “
Alignment of particles in microfluidic systems using standing surface acoustic waves
,”
Appl. Phys. Lett.
92
,
044104
(
2008
).
20.
M.
Hennig
,
J.
Neumann
,
A.
Wixforth
,
J. O.
Rädler
, and
M. F.
Schneider
, “
Dynamic patterns in a supported lipid bilayer driven by standing surface acoustic waves
,”
Lab Chip
9
,
3050
3053
(
2009
).
21.
S. R.
Heron
,
R.
Wilson
,
S. A.
Shaffer
,
D. R.
Goodlett
, and
J. M.
Cooper
, “
Surface acoustic wave nebulization of peptides as a microfluidic interface for mass spectrometry
,”
Anal. Chem.
82
,
3985
3989
(
2010
).
22.
J.
Ho
,
M. K.
Tan
,
D. B.
Go
,
L. Y.
Yeo
,
J. R.
Friend
, and
H.-C.
Chang
, “
Paper-based microfluidic surface acoustic wave sample delivery and ionization source for rapid and sensitive ambient mass spectrometry
,”
Anal. Chem.
83
,
3260
3266
(
2011
).
23.
Y.
Huang
,
S. H.
Yoon
,
S. R.
Heron
,
C. D.
Masselon
,
J. S.
Edgar
,
F.
Tureček
, and
D. R.
Goodlett
, “
Surface acoustic wave nebulization produces ions with lower internal energy than electrospray ionization
,”
J. Am. Soc. Mass Spectrom.
23
,
1062
1070
(
2012
).
24.
A. A.
Kolomenskii
,
H. A.
Schuessler
,
V. G.
Mikhalevich
, and
A. A.
Maznev
, “
Interaction of laser-generated surface acoustic pulses with fine particles: Surface cleaning and adhesion studies
,”
J. Appl. Phys.
84
,
2404
2410
(
1998
).
25.
M. V.
Shugaev
,
A. J.
Manzo
,
C.
Wu
,
V. Y.
Zaitsev
,
H.
Helvajian
, and
L. V.
Zhigilei
, “
Strong enhancement of surface diffusion by nonlinear surface acoustic waves
,”
Phys. Rev. B
91
,
235450
(
2015
).
26.
L. V.
Zhigilei
and
H.
Helvajian
, “
Acoustic processes in materials
,”
MRS Bull.
44
,
345
349
(
2019
).
27.
E. A.
Zabolotskaya
, “
Nonlinear propagation of plane and circular Rayleigh waves in isotropic solids
,”
J. Acoust. Soc. Am.
91
,
2569
2575
(
1992
).
28.
D. J.
Shull
,
M. F.
Hamilton
,
Y. A.
Il’insky
, and
E. A.
Zabolotskaya
, “
Harmonic generation in plane and cylindrical nonlinear Rayleigh waves
,”
J. Acoust. Soc. Am.
94
,
418
427
(
1993
).
29.
M. F.
Hamilton
,
Y. A.
Il’insky
, and
E. A.
Zabolotskaya
, “
Local and nonlocal nonlinearity in Rayleigh waves
,”
J. Acoust. Soc. Am.
97
,
882
890
(
1995
).
30.
M. F.
Hamilton
,
Y. A.
Il’insky
, and
E. A.
Zabolotskaya
, “
Evolution equations for nonlinear Rayleigh waves
,”
J. Acoust. Soc. Am.
97
,
891
897
(
1995
).
31.
M. F.
Hamilton
,
Y. A.
Il’inskii
, and
E. A.
Zabolotskaya
, “
Nonlinear surface acoustic waves in crystals
,”
J. Acoust. Soc. Am.
105
,
639
651
(
1999
).
32.
A.
Lomonosov
,
V. G.
Mikhalevich
,
P.
Hess
,
E. Y.
Knight
,
M. F.
Hamilton
, and
E. A.
Zabolotskaya
, “
Laser-generated nonlinear Rayleigh waves with shocks
,”
J. Acoust. Soc. Am.
105
,
2093
2096
(
1999
).
33.
N.
Kalyanasundaram
, “
Nonlinear surface acoustic waves on an isotropic solid
,”
Int. J. Eng. Sci.
19
,
279
286
(
1981
).
34.
N.
Kalayanasundaram
, “
Nonlinear mode coupling of surface acoustic waves on an isotropic solid
,”
Int. J. Eng. Sci.
19
,
435
441
(
1981
).
35.
N.
Kalyanasundaram
,
R.
Ravindran
, and
P.
Prasad
, “
Coupled amplitude theory of nonlinear surface acoustic waves
,”
J. Acoust. Soc. Am.
72
,
488
493
(
1982
).
36.
R. W.
Lardner
, “
Nonlinear surface waves on an elastic solid
,”
Int. J. Eng. Sci.
21
,
1331
1342
(
1983
).
37.
R. W.
Lardner
, “
Nonlinear Rayleigh waves: Harmonic generation, parametric amplification and thermoviscous damping
,”
J. Appl. Phys.
55
,
3251
3260
(
1984
).
38.
R. W.
Lardner
, “
Waveform distortion and shock development in nonlinear Rayleigh waves
,”
Int. J. Eng. Sci.
23
,
113
118
(
1985
).
39.
M.
Planat
, “
Multiple scale analysis of the nonlinear surface acoustic wave propagation in anisotropic crystals
,”
J. Appl. Phys.
57
,
4911
4915
(
1985
).
40.
R. W.
Lardner
, “
Nonlinear surface acoustic waves on an elastic solid for general anisotropy
,”
J. Elast.
16
,
63
73
(
1986
).
41.
R. W.
Lardner
and
G. E.
Tupholme
, “
Nonlinear surface waves on cubic materials
,”
J. Elast.
16
,
251
265
(
1986
).
42.
V. E.
Gusev
,
W.
Lauriks
, and
J.
Thoen
, “
Theory for the time evolution of nonlinear Rayleigh waves in an isotropic solid
,”
Phys. Rev. B
55
,
9344
9347
(
1997
).
43.
V. E.
Gusev
,
W.
Lauriks
, and
J.
Thoen
, “
New evolution equations for the nonlinear surface acoustic waves on an elastic solid of general anisotropy
,”
J. Acoust. Soc. Am.
103
,
3203
3215
(
1998
).
44.
V. P.
Reutov
, “
Use of the averaged variational principle for describing multiwave interactions of elastic surface waves
,”
Radiophys. Quant. Electron.
16
,
1307
1316
(
1973
).
45.
L. D.
Landau
and
E. M.
Lifshitz
,
Theory of Elasticity
(
Pergamon Press
,
New York
,
1986
).
46.
A. J.
Manzo
and
H.
Helvajian
, “
Demonstration of enhanced surface mobility of adsorbate cluster species by surface acoustic wave excitation induced by a pulsed laser
,”
Proc. SPIE
8969
,
896908
(
2014
).
47.
C.
Wu
,
V. Y.
Zaitsev
, and
L. V.
Zhigilei
, “
Acoustic enhancement of surface diffusion
,”
J. Phys. Chem. C
117
,
9252
9258
(
2013
).
48.
D. R.
Denison
, “
Phonic desorption
,”
J. Vac. Sci. Technol.
6
,
214
217
(
1969
).
49.
C.
Krischer
and
D.
Lichtman
, “
Observation of desorption from quartz induced by surface acoustic waves
,”
Phys. Lett. A
44
,
99
100
(
1973
).
50.
N. G.
Basov
,
E. M.
Belenov
,
M. A.
Gubin
,
M. S.
Kurdoglyan
,
V. V.
Nikitin
,
A. N.
Oraevskĭ
, and
B. N.
Chichkov
, “
New ways of obtaining cold atoms and molecules
,”
Sov. J. Quantum Electron.
17
,
919
922
(
1987
).
51.
H.
Nishiyama
,
N.
Saito
,
H.
Chou
,
K.
Sato
, and
Y.
Inoue
, “
Effects of surface acoustic waves on adsorptive properties of ZnO and NiO thin films deposited on ferroelectric substrates
,”
Surf. Sci.
433–435
,
525
528
(
1999
).
52.
Y.
Inoue
,
Y.
Matsukawa
, and
K.
Sato
, “
Effect of surface acoustic wave generated on ferroelectric support upon catalysis
,”
J. Am. Chem. Soc.
111
,
8965
8966
(
1989
).
53.
S.
Kelling
,
T.
Mitrelias
,
Y.
Matsumoto
,
V. P.
Ostanin
, and
D. A.
King
, “
Acoustic wave enhancement of the catalytic oxidation of carbon monoxide over Pt {110}
,”
J. Chem. Phys.
107
,
5609
5612
(
1997
).
54.
Y.
Inoue
, “
Effects of acoustic waves-induced dynamic lattice distortion on catalytic and adsorptive properties of metal, alloy and metal oxide surfaces
,”
Surf. Sci. Rep.
62
,
305
336
(
2007
).
55.
Y.
Inoue
, “
Acoustic enhancement of surface reactions
,”
MRS Bull.
44
,
361
371
(
2019
).
56.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon Press
,
Oxford
,
1987
).
57.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulations: From Algorithms to Applications
(
Academic Press
,
San Diego
,
1996
).
58.
L. V.
Zhigilei
,
Z.
Lin
,
D. S.
Ivanov
,
E.
Leveugle
,
W. H.
Duff
,
D.
Thomas
,
C.
Sevilla
, and
S. J.
Guy
, “
Atomic/molecular-level simulations of laser-materials interactions
,” in
Laser-Surface Interactions for New Materials Production: Tailoring Structure and Properties
, Springer Series in Materials Science Vol. 130, edited by
A.
Miotello
and
P. M.
Ossi
(
Springer Verlag
,
New York
,
2010
), pp.
43
79
.
59.
B. L.
Holian
and
G. K.
Straub
, “
Molecular dynamics of shock waves in one-dimensional chains
,”
Phys. Rev. B
18
,
1593
1608
(
1978
).
60.
V.
Zhakhovskiĭ
,
S.
Zybin
,
K.
Nishihara
, and
S.
Anisimov
, “
Shock wave structure in Lennard-Jones crystal via molecular dynamics
,”
Phys. Rev. Lett.
83
,
1175
1178
(
1999
).
61.
B. L.
Holian
and
P. S.
Lomdahl
, “
Plasticity induced by shock-waves in nonequilibrium molecular-dynamics simulations
,”
Science
280
,
2085
2088
(
1998
).
62.
K.
Kadau
,
T. C.
Germann
,
P. S.
Lomdahl
, and
B. L.
Holian
, “
Microscopic view of structural phase transitions induced by shock waves
,”
Science
296
,
1681
1684
(
2002
).
63.
V. V.
Zhakhovsky
,
M. M.
Budzevich
,
N. A.
Inogamov
,
I. I.
Oleynik
, and
C. T.
White
, “
Two-zone elastic-plastic single shock waves in solids
,”
Phys. Rev. Lett.
107
,
135502
(
2011
).
64.
L.
Koči
,
E. M.
Bringa
,
D. S.
Ivanov
,
J.
Hawreliak
,
J.
McNaney
,
A.
Higginbotham
,
L. V.
Zhigilei
,
A. B.
Belonoshko
,
B. A.
Remington
, and
R.
Ahuja
, “
Simulation of shock-induced melting of Ni using molecular dynamics coupled to a two-temperature model
,”
Phys. Rev. B
74
,
012101
(
2006
).
65.
D. S.
Ivanov
,
L. V.
Zhigilei
,
E. M.
Bringa
,
M.
De Koning
,
B. A.
Remington
,
M. J.
Caturla
, and
S. M.
Pollaine
, “
Molecular dynamics simulations of shocks including electronic heat conduction and electron-phonon coupling
,”
AIP Conf. Proc.
706
,
225
228
(
2004
).
66.
T.
Hatano
, “
Dislocation nucleation in shocked fcc solids: Effects of temperature and preexisting voids
,”
Phys. Rev. Lett.
93
,
085501
(
2004
).
67.
A. M.
Dongare
,
A. M.
Rajendran
,
B.
LaMattina
,
M. A.
Zikry
, and
D. W.
Brenner
, “
Atomic scale studies of spall behavior in nanocrystalline Cu
,”
J. Appl. Phys.
108
,
113518
(
2010
).
68.
S.-N.
Luo
,
Q.
An
,
T. C.
Germann
, and
L.-B.
Han
, “
Shock-induced spall in solid and liquid Cu at extreme strain rates
,”
J. Appl. Phys.
106
,
013502
(
2009
).
69.
A.
Strachan
,
T.
Çağın
, and
W. A.
Goddard
 III
, “
Critical behavior in spallation failure of metals
,”
Phys. Rev. B
63
,
060103
(
2001
).
70.
V. Y.
Klimenko
and
A. N.
Dremin
, “
Structure of a shock-wave front in a solid
,”
Sov. Phys. Dokl.
25
,
288
289
(
1980
).
71.
V. V.
Zhakhovskii
,
K.
Nishihara
, and
S. I.
Anisimov
, “
Shock wave structure in dense gases
,”
J. Exp. Theor. Phys. Lett.
66
,
99
105
(
1997
).
72.
S. I.
Anisimov
,
V. V.
Zhakhovskii
, and
V. E.
Fortov
, “
Shock wave structure in simple liquids
,”
J. Exp. Theor. Phys. Lett.
65
,
755
761
(
1997
).
73.
M. M.
Budzevich
,
V. V.
Zhakhovsky
,
C. T.
White
, and
I. I.
Oleynik
, “
Evolution of shock-induced orientation-dependent metastable states in crystalline aluminum
,”
Phys. Rev. Lett.
109
,
125505
(
2012
).
74.
V. V.
Zhakhovsky
,
M. M.
Budzevich
,
A. C.
Landerville
,
I. I.
Oleynik
, and
C. T.
White
, “
Laminar, cellular, transverse, and multiheaded pulsating detonations in condensed phase energetic materials from molecular dynamics simulations
,”
Phys. Rev. E
90
,
033312
(
2014
).
75.
C.
Schäfer
,
H. M.
Urbassek
,
L. V.
Zhigilei
, and
B. J.
Garrison
, “
Pressure-transmitting boundary conditions for molecular-dynamics simulations
,”
Comp. Mater. Sci.
24
,
421
429
(
2002
).
76.
E. T.
Karim
,
M.
Shugaev
,
C.
Wu
,
Z.
Lin
,
R. F.
Hainsey
, and
L. V.
Zhigilei
, “
Atomistic simulation study of short pulse laser interactions with a metal target under conditions of spatial confinement by a transparent overlayer
,”
J. Appl. Phys.
115
,
183501
(
2014
).
77.
M. V.
Shugaev
,
I.
Gnilitskyi
,
N. M.
Bulgakova
, and
L. V.
Zhigilei
, “
Mechanism of single-pulse ablative generation of laser-induced periodic surface structures
,”
Phys. Rev. B
96
,
205429
(
2017
).
78.
C.-Y.
Shih
,
M. V.
Shugaev
,
C.
Wu
, and
L. V.
Zhigilei
, “
Generation of subsurface voids, incubation effect, and formation of nanoparticles in short pulse laser interactions with bulk metal targets in liquid: Molecular dynamics study
,”
J. Phys. Chem. C
121
,
16549
16567
(
2017
).
79.
A. V.
Bolesta
,
L.
Zheng
,
D. L.
Thompson
, and
T. D.
Sewell
, “
Molecular dynamics simulations of shock waves using the absorbing boundary condition: A case study of methane
,”
Phys. Rev. B
76
,
224108
(
2007
).
80.
A. A.
Kolomenskii
,
V. A.
Lioubimov
,
S. N.
Jerebtsov
, and
H. A.
Schuessler
, “
Nonlinear surface acoustic wave pulses in solids: Laser excitation, propagation, interactions (invited)
,”
Rev. Sci. Instrum.
74
,
448
452
(
2003
).
81.
S. D.
Stoddard
and
J.
Ford
, “
Numerical experiments on the stochastic behavior of a Lennard-Jones gas system
,”
Phys. Rev. A
8
,
1504
1512
(
1973
).
82.
C.
Taillan
,
N.
Combe
, and
J.
Morillo
, “
Nanoscale self-organization using standing surface acoustic waves
,”
Phys. Rev. Lett.
106
,
076102
(
2011
).
83.
C.
Taillan
,
N.
Combe
, and
J.
Morillo
, “
Chladni figures at the nanoscale
,”
Eur. Phys. J. B
88
,
317
(
2015
).
84.
R.
Stoneley
, “
The propagation of surface elastic waves in a cubic crystal
,”
Proc. R. Soc. London Ser. A
232
,
447
458
(
1955
).
85.
J. R.
Ray
and
A.
Rahman
, “
Statistical ensembles and molecular dynamics studies of anisotropic solids
,”
J. Chem. Phys.
80
,
4423
4428
(
1984
).
86.
M. I.
Rabinovich
and
D. I.
Trubetskov
,
Oscillations and Waves in Linear and Nonlinear Systems
(
Kluver Academic Publishers
,
Dordrecht
,
1989
).
87.
N. M.
Ryskin
and
D. I.
Trubetskov
,
Nonlinear Waves
(
Fizmatlit
,
Moscow
,
2010
) (in Russian).
88.
A.
Lomonosov
and
P.
Hess
, “
Effects of nonlinear elastic surface pulses in anisotropic silicon crystals
,”
Phys. Rev. Lett.
83
,
3876
3879
(
1999
).
89.
C.
Zener
, “
Internal friction in solids II. General theory of thermoelastic internal friction
,”
Phys. Rev.
53
,
90
99
(
1938
).
90.
D. S.
Ivanov
and
L. V.
Zhigilei
, “
Combined atomistic-continuum modeling of short-pulse laser melting and disintegration of metal films
,”
Phys. Rev. B
68
,
064114
(
2003
).
91.
T.
Halicioğlu
and
G. M.
Pound
, “
Calculation of potential energy parameters from crystalline state properties
,”
Phys. Status Solidi A
30
,
619
623
(
1975
).
92.
S.
Zhen
and
G. J.
Davies
, “
Calculation of the Lennard-Jones n-m potential energy parameters for metals
,”
Phys. Status Solidi A
78
,
595
605
(
1983
).
93.
Y. A.
Pishchal'nikov
,
O. A.
Sapozhnikov
, and
V. A.
Khokhlova
, “
A modification of the spectral description of nonlinear acoustic waves with discontinuities
,”
Acoust. Phys.
42
,
362
367
(
1996
).
94.
E. Y.
Knight
,
M. F.
Hamilton
,
Y. A.
Il’inskii
, and
E. A.
Zabolotskaya
, “
On Rayleigh wave nonlinearity, and analytical approximation of the shock formation distance
,”
J. Acoust. Soc. Am.
102
,
2529
2535
(
1997
).
95.
O. V.
Rudenko
and
S. I.
Soluyan
,
Theoretical Foundations of Nonlinear Acoustics
(
Plenum Publishing Corporation
,
New York
,
1977
).
96.
D. T.
Blackstock
, “
Nonlinear acoustics (theoretical)
,” in
American Institute of Physics Handbook
, 3rd ed. (
McGraw-Hill
,
New York
,
1972
) Chap. 3n.
97.
L. K.
Zarembo
and
V. A.
Krasil’nikov
, “
Nonlinear phenomena in the propagation of elastic waves in solids
,”
Sov. Phys. Usp.
13
,
778
797
(
1971
).
98.
J. K.
Nørskov
,
T.
Bligaard
,
B.
Hvolbæk
,
F.
Abild-Pedersen
,
I.
Chorkendorff
, and
C. H.
Christensen
, “
The nature of the active site in heterogeneous metal catalysis
,”
Chem. Soc. Rev.
37
,
2163
2171
(
2008
).
99.
B.
von Boehn
,
M.
Foerster
,
M.
von Boehn
,
J.
Prat
,
F.
Macià
,
B.
Casals
,
M. W.
Khaliq
,
A.
Hernández-Mínguez
,
L.
Aballe
, and
R.
Imbihl
, “
On the promotion of catalytic reactions by surface acoustic waves
,”
Angew. Chem. Int. Ed.
(accepted) (
2020
).

Supplementary Material