Nonreciprocal wave propagation and parametric amplification of bulk elastic waves in nonlinear anisotropic materials

Parametric amplification of an elastic wave and a framework for using elastic waves that could enable a new generation of high performance, low noise acoustic amplifiers, mixers and circulators are presented. Using a novel approach with nonlinear materials produces highly desirable non-reciprocal characteristics. Parametric amplification of a weak elastic signal wave is achieved by an elastic pump wave of higher intensity. By careful selection of material orientation together with precise excitation of signal and pump waves, ‘up frequency conversion’ is suppressed and selective amplification of the elastic signal wave occurs at its original frequency. In addition, a general mathematical framework is developed and used for analytical studies of coupled wave equations in nonlinear anisotropic materials. The results obtained from the analytical studies are verified using a finite element implementation.


Introduction
Most electronic handheld consumer products today take advantage of surface and bulk acoustic wave (SAW, BAW) filters and delay lines [1,2]. These provide the advantage of small size due to slower wave velocity at radio frequency (RF) ranges relative to electromagnetic waves, and high efficiency due to reduction of resistive losses compared with similar purely electronic devices. They are compatible with existing integrated circuit (IC) fabrication techniques and can be integrated with other circuit elements [3]. However, as first stated by Helmholtz in 1859 and proved by Rayleigh in 1878, the propagation of acoustic waves in conventional linear media is reciprocal [4]; which is 'requiring the transmission of information or energy between any two points in space to be symmetric for opposite propagating directions' [5]. This reciprocal characteristic of acoustic waves limits their applications where directional dependency is desirable.
Creation of non-reciprocity has been mainly achieved with three approaches [6]. Spatiotemporal modulation of some elements of the system, applying an external symmetry breaking field such as an applied magnetic field, or utilizing nonlinear behavior of the system.
Time and space modulation can be applied to material properties or boundary conditions of the system to break reciprocity. In 2015, Swinteck et al [7] applied a light source with time and space variant intensity to a material with a large photo-elastic coupling to modulate its elastic properties. Spatiotemporal modulation of the elastic constants of the material produced a time dependent superlattice which demonstrated nonreciprocal propagation of a bulk elastic wave. Non-reciprocal propagation of an elastic wave in a beam with spatiotemporal modulation of its Young's modulus and density was investigated by Trainiti et al [8]. It was shown with both time and space modulation, the dispersion diagrams for this system were no longer symmetric with respect to the frequency axis and directional band gaps were created. Croenne et al [9] used spatiotemporal modulation of the electrical boundary condition applied to a periodically repeated assembly of piezoelectric material sandwiched between thin metallic electrode layers. They showed nonreciprocal transmission of an input longitudinal acoustic wave. Their results showed scattering effects such as frequency conversion and generation of harmonics. Several recent works have also reported experimental realization of reciprocity breaking in a time modulated system [5,10].
Another approach to produce non-reciprocal behavior is to apply an external symmetry breaking bias field. External magnetic field bias is commonly used [11][12][13][14] although similar symmetry breaking has been demonstrated in a linear acoustic device with a circulating fluid that creates an angular-momentum bias [15] and by using magneto-elastic coupling to create a gyrator [16]. Wang et al [17] demonstrated breaking time reversal symmetry using gyroscopic inertial effects that creates an apparent external force. Adding a spinning gyroscope to each lattice site, they showed additional topological bands are created that enables multimode propagation of an elastic wave on the edge of the material. While external field biasing has been theoretically and experimentally shown to be effective in some applications, it may not be desirable in terms of physical packaging, fabrication and increased dimensionality of the system.
Another prolific area, seen frequently in phononics and metamaterials research, uses material nonlinearity and asymmetry to break reciprocity [18][19][20][21][22]. Liang et al [18,19] demonstrated acoustic rectification by asymmetrically coupling a super lattice to a nonlinear medium. The nonlinear mechanism, however, do not break reciprocity at the fundamental frequency. Non-reciprocity is realized in the total acoustic flux at the boundaries. Non-reciprocal acoustic propagation in which the frequency of the incident wave was preserved has been demonstrated experimentally in a system composed of a granular chain and a conical rod at low frequencies [23]. Other works have investigated nonlinear material with hierarchal asymmetry [24][25][26]. Moore et al [25] showed breaking reciprocity within a unit cell featuring a hierarchical internal nonlinearity imposing one directional transfer of energy from larger to smaller scale. Fronk et al [26] extended this asymmetry in a lattice of non-reciprocal unit cells and showed the non-reciprocity at global scale. The general concept of using nonlinearity and asymmetry to break reciprocity is a common theme in the cited references; however, the source of nonlinearity and the asymmetry elements vary.
Each of the methods used to produce non-reciprocity has various strengths and weakness and therefore tend to be applied to specific application areas. For example, much of the work to date has focused on photonic devices or macroscale acoustic devices but has not addressed RF range applications. Other approaches rely on coupled resonators or construction techniques that are not compatible with current IC fabrication techniques. While many works have reported on the study of wave propagation in isotropic [27] or anisotropic nonlinear elastic materials [28][29][30], to the authors knowledge, none have shown potential for parametric amplification in a nonreciprocal RF application. Due to the high-quality factor of mechanical resonances, development of nonreciprocal amplification devices operating at RF frequencies based on elastic materials may prove to be a revolutionary concept and represents the area targeted by this paper. Here the framework to enable the investigation of such devices is derived.
This paper shows non-reciprocal parametric amplification and non-reciprocal propagation of bulk elastic waves in a homogenous anisotropic material. This is due to second order material and geometric nonlinearity combined with an elastic traveling pump wave introduced in the medium as a symmetry breaking element. In this system, non-reciprocal parametric amplification of a bulk elastic 'signal' wave is demonstrated. Further, it is shown that propagation of the 'signal' wave traveling through the system is non-reciprocal. The non-reciprocal propagation appears as a difference in the intensity of the 'signal' wave traveling with, versus traveling opposite, to the direction of the pump wave (i.e. energy exchange occurs preferentially). This difference in intensity is due to parametric amplification of the signal traveling in the direction of the pump wave but not in the opposite direction. Providing the required phase matching condition for parametric amplification process and eliminating the phase matching condition for the higher frequency component (sum frequency generation (SFG)), is the key to achieving parametric amplification.
A general derivation of the coupled-wave formalism for elastic waves of different frequencies is presented and used to identify the requirements for parametric amplification and the associated phase matching conditions. No restrictions are made on material anisotropy or direction of wave propagation. The methodologies used to derive the coupled wave equations are adopted from the field of nonlinear photonics [31]. A condition of slow variation of the elastic field amplitude over the distance of an elastic wavelength is assumed, aka slowly varying amplitude (SVA) assumption.
The coupled wave equations are simplified to provide an analytic solution for collinear propagation of pump and signal plane waves to investigate the process of parametric difference (down) frequency conversion. The theory predicts parametric amplification of an acoustic signal wave through difference frequency conversion. The numerical results and conclusions arrived at in the analytical theory are further investigated using a finite element method (FEM) model to simulate the propagation and interaction of bulk elastic waves in a nonlinear anisotropic medium. Parametric amplification and difference frequency conversion are observed in the computational implementation, in agreement with the analytic results.

Linear propagation of plane elastic waves and modal expansions
This section begins with a review of the linear theory of elastic plane waves propagating in a bulk anisotropic material. The modal expansion of plane wave solutions derived in this section is used in the development of the nonlinear theory.
The approach is to define the constitutive properties in the global, crystal, coordinate system, select a propagation direction, then determine the type of wave propagation that is possible in this propagation direction. In an isotropic material, this will lead to a bulk longitudinal wave propagating at bulk velocity and a shear wave propagating at shear wave velocity where the shear velocity is degenerate, i.e. independent of the displacement direction in the plane orthogonal to the propagation direction.
As described in detail in appendix A, in the absence of body forces and under a small strain-displacement assumption, a linear elastic wave with displacement field in its complex form can be written as equation (1.1) Following the methodology described in [30], solutions of the above equation can be expressed as 3 is a unit vector in the direction of propagation defined relative to the global coordinate system and considered to be fixed in space, and c is the wave speed (phase velocity). Substituting equation is the symmetric Christoffel acoustic tensor. Equation (1.3) is an eigenvalue problem inŪ and for non-trivial solutions the determinant of the multiplying matrix must equal to zero, i.e.
which gives a cubic polynomial in terms of rc . 2 Each eigenvalue, a c determined from the solutions of equation (1.5) has a corresponding eigenvector¯a U that is found from equation  Here and throughout the paper, parentheses on repeated indices, e.g. (a), are used to indicate no Einstein summation. The orthonormal eigenvectors¯a l form a 3D basis set and can be used for the expansion of any vector. Expressed in the eigenvector basis, the components of displacement are represented in equation (1.7)¯( The superscript on displacement indicate components referenced to the modal basis while the subscripted u refers to the global or material system. Unlike isotropic materials which have one distinct longitudinal and two degenerate shear speeds, anisotropic materials generally have three distinct eigenvalues or phase velocities, ( ) a c . Consequently, each shear mode can have a different wave velocity. This phenomenon is called 'Birefringence' [32]. Additionally, the modes¯a l are not necessarily purely longitudinal or purely transverse with respect to the propagation direction. However, one mode is predominately in the direction of propagation and is called the quasi-longitudinal mode and the two others are predominately normal to the direction of propagation and are called the quasi-shear modes.
Defining the scalar quantity x as the inner product of thep andx,¯( we can assume the modes are exponential functions of the scalar quantity x as equation (1.9) c is the magnitude of the wave vector for a wave in mode¯a l . A general wave with particle displacement not aligned with any of the eigen-basis vectors, does not have a single value for its propagation velocity but can be expanded in its modal form, equation (1.10), with each of its components having the velocity of the corresponding modal displacement The presumption that the displacement field of a wave is only in direction of one modal basis is called the mono-mode assumption. In this case only one phase velocity for the wave is excited.
Once a crystal orientation and a propagation direction relative to the crystal orientation has been specified, orthonormal eigenvectors¯a l are found that form a three-dimensional basis set for that propagation direction and can be used to expand any plane wave displacement vector propagating in that direction in the crystal. Consequently, the modal expansion of equation (1.7), while based on the eigen solution of the linear material constitutive properties, can be applied to both the linear and nonlinear systems. In the next section, we will apply this expansion to develop the nonlinear equations.

Second order nonlinear processes, parametric frequency conversion
In this section, the different second order frequency conversion processes are introduced and the coupled wave equations are represented for these processes. Only quadratic nonlinearity is investigated, since the order of the magnitude of the strain terms in typical BAW or SAW devices are small and adding higher order terms will not add a significant contribution to the behavior of the system. The final form of the equation is given for an arbitrary number of frequencies. Investigation of the parametric amplification process in an anisotropic nonlinear elastic material is based on the coupled wave equations discussed in this section.
A process is called parametric when there is no net exchange of energy between a traveling wave and the medium in which the wave is propagating. Consequently, the sum total energy of all the waves traveling in the medium is conserved. In addition, linear wave processes do not exchange energy between different frequency components. Hence, in a linear parametric process there is no exchange of energy between the waves and the medium or the wave components themselves. In a nonlinear medium, waves at different frequencies can exchange energy among themselves due to frequency mixing, or frequency conversion processes. In a nonlinear parametric process, even though there is no exchange of energy between the waves and the medium, different frequency components couple and energy transfers among them. For instance, in the presence of a weak signal wave and another wave at a different frequency propagating in a nonlinear medium, the energy transfer between the two waves can lead to parametric amplification of the weak signal wave with a corresponding reduction in the intensity of the other wave. Second order nonlinearities in elastic materials are described by nonlinear terms in the constitutive equations and by second order displacement gradients in the definition of the Lagrangian strain. These are the material and the geometric nonlinearities. As shown in appendix B, these two effects can be combined in the form of nonlinear elastic wave equations represented by equation (2.1) Complex displacement fields, useful for decomposing the frequency mixing equations, can be represented through the relation to the real fields as equation (2.2)  Since the relation in equation (2.7) must hold for all times, the exponential terms are collected so that the same frequency exists in each summand. This requires w w w = + q r s where r, q, and s can take on negative values. The negative values of r, q, and s are interpreted as appropriate to satisfy the conjugate frequency notation of equation (2.6). This implies equation The summation in equation (2.8) is over all frequency combinations that satisfy the constraint w w w = + q r s for each w q independently. The notation ( ) r s , here is used to represents a pair that satisfies this constraint and r is not independent of s. Equation (2.8) can be used to describe any plausible parametric frequency mixing process.
Parametric frequency mixing or frequency conversion, in a second order nonlinear process involving three frequencies can be classified into three categories. SFG or up conversion, second harmonic generation (SHG), and difference frequency generation (DFG) or down conversion. In SFG, two input waves, say w 1 and w , 2 generate a third higher frequency wave w w w = + .
In this equation, the component ( ) a u is the amplitude of the mode a,¯( ) a l , and defines the direction of 'particle' displacement (for example quasi-longitudinal, shear, etc). If we assume the displacement field at frequency w q is composed of only one mode, i.e. its particle displacements are aligned with only one of the eigen basis, we can simplify equations (2.9a)-(2.9b) This is referred to as the mono-mode assumption. If the variation of the amplitude of a wave over its wavelength is small, the SVA assumption applies. The SVA assumes the spatial dependence of the amplitude function for a wave can be decomposed into the product of a slowly varying envelope term, or amplitude, and a harmonic function capturing its oscillations. Under the SVA assumption the wave components at different frequencies, ( ) w a u , q propagating in directionp, i.e. along x, have the form of equation (2.10), is the SVA term and ( ) a k q is the magnitude of the wave vector for the wave component at frequency w q in mode a.
As it is described in detail in appendix B, with above assumptions, equation (2.8), leads to a system of coupled equations in which the amplitude of the displacement field at frequency w , q couples to the amplitude of the fields at frequency w r and w , s as long as they satisfy the condition w Here, C ijkl and C jiklmn are the components of the rank four and six stiffness tensors, respectively, and jinl km is the effective nonlinear rank six tensor (see appendix B). Equation (2.11) is the most general form of the coupled wave equations for a frequency mixing process in a nonlinear elastic material under the mono-mode and SVA assumptions, regardless of the number of frequencies involved. In general, analytic solutions of equation (2.11) do not exist and numerical methods must be applied for a specific frequency mixing process. In section 5 analytic solutions of equation (2.11) for the process of DFG involving three frequency components is given and its application to parametric amplification process is discussed. Before completing the analytic solutions, in the next section phase matching conditions are discussed to motivate the rationale for the DFG analytic solutions.

Phase matching condition
In general, when signal w sig and pump w p waves are supplied to the medium, conservation of elastic energy, with w w w = + , p SFG sig implies that with SFG, one phonon at signal frequency and another at pump frequency must annihilate simultaneously (combine) to generate a phonon at SFG. In DFG, w w w = -, p sig i a phonon at signal and idler frequencies are generated when a phonon at pump frequency annihilates. Both these processes can happen simultaneously when the pump and signal waves are input at the boundary. When parametric amplification of the signal wave is of interest, generation of a sum frequency wave reduces the intensity of both pump and signal waves, which results in an overall reduction in the intensity of the signal wave instead of amplifying it. Therefore, it is significant to know what are the conditions that allow a frequency mixing process to happen. This section provides an investigation of these effects.
In a process involving only three wave components such that w w w = + , q r s the change in intensity of a wave at frequencies w q with respect to x, has the form of equation (3.1) (see appendix C), , is a positive quantity that depends on the material properties and the wave intensities and ( ) y x In this equation, f c is a function of the effective material properties and is either 0 or p. Dk rsq is the phase mismatch in equation (2.11), and f , r f s and f q are the phases of the wave envelopes at frequenciesw , r w s and w , q respectively. The phasor of the wave envelope is represented by shows that the derivative of the intensity with respect to the propagation distance x changes sign periodically, with periodicity p y . While ( ) y > Cos 0, the frequency mixing process transfers energy from the wave components at frequencies w r and w s to the w q wave and reverses energy the flow when ( ) y Cos changes sign. Hence, elastic energy flows back and forth among different frequency components over a distance p y . The interaction distance before the frequency conversion process is reversed is called the coherence length.
In general, to maximize the coherence length, y should be minimized. y is composed of two parts, the sum of the wave envelope phases and the phase mismatch factor Dk .
rsq For the first part, r s q is determined by the material properties and the relative phases of the two input waves at frequency w r and w .
s In practice, these can be adjusted at the boundary, Consequently, the more significant factor to consider is Dk . rsq From equations (3.2) and (3.3), it is seen that with a smaller phase mismatch, larger coherence length is achieved. The ideal situation happens when the Dk rsq term goes to zero; this is called the phase matching condition.
In a second order, nonlinear frequency conversion process involving three frequencies such that w w w = + , q r s the phase matching condition is written as equation (3.4) Under the collinear propagation assumption, the phase matching condition, equation (3.4), turns into a scalar relation among the magnitude of the wave vectors as equation (3.5a), or in terms of phase velocities, = w c , k with mono-mode assumption, as equation (3.5b). ( Once the pump and signal waves are launched in a nonlinear media both SFG and DFG occur simultaneously. Considering w q to be the signal wave at frequency w , sig such that w w w = - where A DFG and A SFG are two positive numbers that can be obtained from equation (C. 22) in appendix C. Initially, at x = 0, the energy transfer from the signal to the sum frequency wave, reduces the intensity of the signal wave until ( ) y Cos SFG changes sign and the process is reversed. However, if the phase mismatch term is large for the process of SFG, the coherence length of this process will be small and its effect can be negligible on the propagation of the signal wave.
If we assume the effect of SFG is negligible, equation (3.6) reduces to (3.7) term changes sign at a larger distance from the origin, i.e. the amplification process has larger coherence length.
From the above discussion, it is concluded that parametric amplification is feasible when the associated phase matching condition is satisfied and the input waves have the appropriate initial phases.

Parametric amplification
Unlike electromagnetic waves, SAW and BAW are relatively non-dispersive over the typical frequency ranges of interest and, for purposes here, the wave velocity c is not a function of the frequency of the propagating wave. In an anisotropic material where the three waves have the same mode, i.e. all displacements are along the same¯a l vector, or in an isotropic material when the waves are all shear or all longitudinal mode, the velocities of the three traveling waves are equal, = = c c c , s r q and equation (3.4) is always satisfied. When this happens, both the DFG and SFG wave components are phase matched which leads to poor amplification or decay of the signal wave.
One way to overcome this obstacle is to utilize the birefringent property of anisotropic materials. In this case two waves propagating in the same direction but having different displacement modes can have different phase velocities. This way by choosing the proper direction of propagation and exciting the pump and signal waves in the desired modes, the phase mismatch term can be maximized for SFG and minimized for DFG. Consequently, the SFG wave has small coherence length and does not interfere significantly with the desired DFG process.
To find the best propagation direction and modal orientation for the pump and signal waves, we start with DFG. Substituting w w w = - Examining equation (4.1), for w p to be greater than w , sig the phase speeds either have to satisfy inequality (4.2a) or (4.2b) Going through the same procedure for SFG, inequalities (4.3a) and (4.3b) are achieved. The set of inequalities in equation (4.4a) always have smaller phase mismatch for the process of DFG. Therefore, solutions satisfying this set of inequalities are the ones sought here.
Assuming the eigen modes are ranked in order of decreasing phase velocity, the first pair of inequalities require the signal and pump waves to be in modesl 1 andl 2 respectively.
While the idler wave generates in all modes, for simplicity, under the mono-mode assumption, only the mode with largest coherence length is tracked and the other two are neglected. The largest coherence for this process occurs when the idler wave is in thel 3 mode. For the process of SFG, the smallest Dk , SFG occurs for itsl 2 mode, therefore, this mode has the largest coherence length and greatest effect on the system. The other modes are assumed to be negligible. With this, the phase mismatch term for the two processes are represented in equations (4.5) and (4.6) For situations when Dk SFG is large and Dk DFG is small, the SFG has negligible effect on the DFG and parametric amplification is significant. The accuracy of this assumption is further studied in the finite element simulation section. In these cases, equation (2.11) can be simplified to account only for signal ( ) w , sig pump( ) w p and idler ( ) w i frequencies resulting in equations (4.7)-(4.9), where all the terms in equations (4.7)-(4.9) are as defined in equation (2.11).
The above system of coupled equations represents three nonlinear equations with three unknowns. The amplitude of the pump wave is generally much greater than the amplitude of the signal wave in envisioned applications; consequently, we assume the magnitude of the pump wave to be almost constant throughout the interaction. This will allow us to treat U ( ) w p 2 as a know function that factors out of the equations for U ( ) w x , 1 sig The accuracy of this assumption is also discussed further in the finite element simulation results section. Under this assumption, equations (4.7) and (4.9) can be written as equations (4.10) and (4.11), respectively. Sinh

Breaking reciprocity
In this section, it is shown in a parametric amplification process when the required phase matching condition is satisfied for pump and signal waves traveling in the same direction, the signal wave is amplified and when the signal wave travels in a direction opposite the pump wave, the required phase matching condition cannot be satisfied. Therefore, the magnitude of the signal wave at the destination is different depending on its direction of propagation. This results in directional dependency of the propagation of the signal wave. It is concluded that in a nonlinear media with a traveling pump wave in one direction, non-reciprocal propagation of an elastic wave can be achieved.
The phase matching condition for the process of DFG for pump, signal and idler waves when, w w w -= ,  figure 1(b), the relationship among their wave vectors cannot satisfy the required phase matching condition and the frequency conversion process cannot take place in this situation. Therefore, without energy transfer between the pump and signal waves, the amplitude of the signal wave remains constant. This directional bias in amplitude of the signal wave is interpreted as nonreciprocity in its propagation. In addition, when the DFG does not happen, the idler frequency does not generate. That means comparing the Fourier spectrum of the elastic wave in the forward direction, figure 1(a), versus the backward direction, figure 1(b), the idler frequency is only observed for the forward direction and is absent in the backward direction.

Results
In this section, linear and nonlinear material properties of LiNbO 3 (appendix D), are used as an example of the concepts and equations derived. Although the piezoelectric properties of LiNbO3 are not considered in this paper, we envision the piezo coupling will be used to produce electromechanical devices similar to the SAW and BAW devices currently in commercial applications. This simplified system is the first step for designing practical devices. LiNbO 3 has a low loss and a high coupling factor that makes it easy to work with for applications at RF frequencies and LiNbO 3 has published data for the higher order coupling terms which are of sufficient magnitude to use for proof of concept.
In the first section, propagation directions with small phase mismatch for DFG and large mismatch for SFG are determined. Directions with a ratio of In the second part of this section, the nonlinear elastic wave equations are solved numerically with FEM simulations. In these simulations, the equations are solved considering both DFG and SFG processes (i.e. four frequency interaction). Further, in the numerical simulations, the SVA and mono-mode assumptions are not used. These results are compared with the analytical solutions, validating previous simplifications and assumptions. Figure 2 refers the direction of the wave propagation vectorp to the crystal axes in LiNbO 3 , with each propagation direction determined by two angles q and j with respect to the positive Z and X axes, respectively. In section 2, it is discussed that in a general anisotropic material, the wave speed for each of the three modes can be distinct. The phase velocity as a function of direction of propagation for each mode in LiNbO 3 is shown in figure 3. Figures 3(a)-(c) shows the wave speed for the quasi-longitudinal and the two quasi-shear modes, associated with the propagation direction in 3D space. In these plots, every point on the surface corresponds to a (q, j) pair that defines the propagation direction and the distance from the origin is the magnitude of the modes' wave speed. In addition, the speed values are color coded as indicated in the bar-legend attached. Figure 3(d) shows the wave speed of all three modes overlaid. In LiNbO 3 , for q =  0 , that is propagation in the Z direction, the two shear velocities are equal. This indicates the isotropic behavior of shear modes oriented in the X-Y plane. We also see from figure 3(d) the two quasi-shear modes have approximately half the longitudinal wave speed.

Analytical results
An optimum direction of propagation is the one with (1) minimum phase mismatch for the process of DFG,   An 'optimum' value for r is not shown here since the r value does not determine gain. Our analytic gain results, being based on three frequency interaction, do not consider the effects of SFG. Consequently, higher values of gain can be found in certain directions, but the assumptions needed for validity of the three frequency results are violated. By choosing low r values, we are considering the gain only in regions where the 3 frequency assumptions reasonably apply. Here we have considered the points with < r 0.2 as a valid metric for when SFG is negligible. This value is verified by the FEM results which show close agreement within this region. For the directions with < r 0.2, the signal gain for a DFG process, neglecting SFG, is plotted in figure 4(b) for the X-Y plane, and In the analytical results section the coordinate system which aligns with the crystal axes are assumed to be fixed in space and the direction of propagation directionp is defined with respect to this coordinate system, figure 2. For convenience in the numerical evaluation, a coordinate transformation is applied to align thep propagation direction with the Z′ axis, and the vectorsp×Z andp×Z×p are chosen to form the Y′ and X′ axes of the FEM coordinate system, respectively', see the figure 5(a).
, we apply these BC in COMSOL using a 3D model. The 3D geometry of the model consists of a bar that is several wavelengths long in the propagation direction z, and meshed with a single element across the width d in the transverse x and y directions. We impose periodic boundary conditions of continuity type to each side of the bar in the transverse ijkl mn and C , ijkl C ijklmn are the second-and third-order elastic constants.
x and y directions. This constrains the displacement components to be constant in x and y dimensions giving a 1D problem in 3D space.
To make sure the results are not influenced by reflection from ¢ Z boundaries, perfectly matched layers are used on both ends to simulate continuous propagation. The signal and pump waves are excited as prescribed displacements at the boundary with initial amplitudes U = 10 nm, Three simulations are presented to show the effect of different parameters on the system. In the first simulations, the importance of exciting the pump and signal waves in the directions suggested by equation (4.4a) are emphasized. In this simulation, pump and signal waves are sent inl 1 andl 3 directions, respectively, in contrast to what equation (4.4a) would suggest. It is shown with this excitation, the phase matching condition is not favorable for DFG and parametric amplification of the signal wave does not occur. In the second simulation, pump and signal waves are applied in thel 1 andl 2 directions, in agreement with equation (4.4a), and parametric amplification of the signal wave is demonstrated. The third simulation shows non-reciprocal propagation of the signal. In the third simulation, this is accomplished by taking the model from the second simulation and reversing the propagation direction of the signal wave. The resultant magnitude of the signal wave is then seen to be quite different from the second simulation.
In the first simulation, the waves propagate in the direction q =  82 and j =  70 and signal and pump waves are exited inl 1 andl 3 directions, respectively, which are not optimum directions to eliminate the generation of a sum frequency wave. Figure 5(b) shows in this case the signal wave does not amplify and its displacement magnitude reduces as it propagates over the first wavelength. This agrees with the results predicted in the analytical discussion.
The analytical coherence lengths, calculated from equation As expected, the plot shows that coherence length of the SFG is larger than the DFG since Dk DFG >Dk .

SFG
In the second simulation, the same geometry and material properties are used, however, a more favorable phase matching condition for DFG is applied. Graphical representation of the pump and signal waves in figure 6(a) shows the signal and pump waves are launched in the medium with modal displacements inl 1 andl 2 directions, respectively. This choice was made based on the inequalities derived in section 5, equation (4.4a). The eliminates its adverse effect on signal amplification and the red plot in figure 6(b) shows the increasing signal amplitude as the wave travels. This demonstrates the parametric amplification of the wave as the result of energy transfer from the pump wave. It is also seen that the amplitude of the pump wave is fairly constant, indicating the approximate, but reasonable, assumption of constant amplitude for the pump. Figure 6(c) plots the total intensity of the combined waves along with the intensity of each wave component as a function of position in the direction of propagation. It is seen that the total intensity of the waves is conserved and the intensity of the pump wave decreases while the intensity of the signal and idler waves increases. The conservation of the total intensity of the involved waves is due to nature of the problem being parametric. Figures 6(d) and (e) compare the plot of the signal and idler waves obtained analytically and numerically over a distance of four signal wavelengths in the material. This figure shows close agreement between analytic and numerical results for about three signal wavelengths, demonstrating the approximate range of validity for the SVA assumption.
The objective of this last simulation was to investigate non-reciprocity in propagation of the signal wave. The same model setup as the second simulation was used, except, the signal wave was sent in the opposite direction of the pump wave. The results plotted in figures 7(b) and (c) show that by reversing the direction of propagation of the signal wave, the pump and signal waves do not couple; consequently, the idler wave is not generated. Also seen, there is no energy transfer from or to the pump and signal waves.
Comparing the second and third simulations shows that the amplitude of the signal wave is quite different over equivalent distances of travel. For the nonlinear elastic media, the pump wave creates a symmetry breaking field and the signal wave propagation becomes non-reciprocal.

Conclusions
This work has shown the derivation and verification of a system of equations that can be used for determining amplification effects due to DFG associated with bulk elastic waves. A simplified version of these equations was solved quasi-analytically for the 'three wave' problem. These idealized equations indicate promising modes of operation for the non-degenerate eigenvalue case where amplification of the signal wave is achieved. This was verified numerically using a FEM model. Without the analytic equations to guide selection of the required system characteristics, i.e. which mode to operate the pump, signal and idler as well as frequency determination, it would be very difficult to create an amplifier due the large number of choices available in the parameter space and the relatively limited number of combinations that can produce parametric amplification. In this paper only the mechanical properties of the Lithium Niobate were considered and its piezoelectric effects were neglected. This simplified system is the first step for designing practical surface acoustic devices and the piezoelectric effects should be addressed in future studies.

Acknowledgments
This work is supported by the National Science Foundation's Emerging Frontiers in Research and Innovation (EFRI) program through award EFRI-1641128.

Appendix A. Linear propagation of plane elastic waves and modal expansions
In the absence of body forces, the equations of motion and small strain-displacement relations for a mechanical system are represented in equations (A.1) and (A.2) The mechanical fields defined in the above relations are real quantities.
=  is the symmetric Cauchy stress tensor, ρ is the material density, is the displacement vector and =  is the small strain tensor. Each field quantity can be represented by a function plus its complex conjugate, equation (A.3) Assuming material properties are instantaneous in time and local in space, the stiffness tensor C ijkl can be written as equation (A6) where d is the Dirac delta function. This allows the linear stress-strain constitutive relations to be expressed in time domain as equation (A.7) , .
Following the methodology described in [30], solutions of the above equation are expressed as equation where, = -I 1 , U k is the wave amplitude,¯( 3 is a unit vector in the direction of propagation defined relative to the fourth order constitutive tensor that is considered to be fixed in space, and c is the wave speed (phase velocity).
Substituting equation (A.11) into (A.10) results in (A.12) is the symmetric tensor called the Christoffel acoustic tensor and the components of p are defined in the crystal coordinate system. Equation (A.12) is an eigenvalue problem inŪ and for non-trivial solutions the determinant of the multiplying matrix must be zero This gives the eigenvalues as roots of a cubic polynomial in rc . 2 Each eigenvalue a c , determined from the solution of equation (A.14) has a corresponding eigenvector¯a U that is found from equation (A.12). Unit vectors This appendix ends with derivation of an identity that is used in appendix B to simplify the nonlinear equations developed. Once a crystal orientation and a propagation direction relative to the crystal orientation has been specified, orthonormal eigenvectors¯a l are found that form a 3D basis set for that propagation direction and can be used to expand any plane wave displacement vector propagating in that direction in the crystal. This relation will be used in appendix B to simplify the derived nonlinear equations.

Appendix B. Nonlinear propagation of plane elastic waves and coupled wave equations
When stress is applied to an element of the medium, the element undergoes a displacement given bȳ , , where the vectorsx andX describe the element position in the deformed and undeformed states, respectively.
Starting with the nonlinear constitutive equations, the second Piola-Kirchoff stress tensor is approximated as a series expansion of the strain energy with respect to the Lagrangian strain,  , ij in equation (B.1) where F is the elastic energy function,C ijkl is the rank four stiffness tensor, that is associated with linear behavior, C ijklmn is the rank six stiffness tensor that represents quadratic nonlinearity and the Lagrangian strain,  , ij which is given in equation (B.2), Here r 0 is the density of the matrial represented in the undeformed coordinate system.
The first Piola-Kirchoff stress is related to the second Piola Kirchoff stress by equation (B.8) As in appendix A, the definition of equation (A.3) is used to introduce the complex displacement field¯(¯) u x, t which transforms equation (B.12) into (B.13) and its complex conjugate The SVA assumption presumes if the variation of the amplitude of a wave over its wavelength is negligibly small, its mathematical representation can be composed of two terms. The envelope or the amplitude of the wave and a sinusoidal function capturing its oscillations. Under the SVA assumption the wave components at different frequencies,