Elastic Wave Eigenmode Solver for Acoustic Waveguides

A numerical solver for the elastic wave eigenmodes in acoustic waveguides of inhomogeneous cross-section is presented. Operating under the assumptions of linear, isotropic materials, it utilizes a finite-difference method on a staggered grid to solve for the acoustic eigenmodes of the vector-field elastic wave equation. Free, fixed, symmetry, and anti-symmetry boundary conditions are implemented, enabling efficient simulation of acoustic structures with geometrical symmetries and terminations. Perfectly matched layers are also implemented, allowing for the simulation of radiative (leaky) modes. The method is analogous to eigenmode solvers ubiquitously employed in electromagnetics to find waveguide modes, and enables design of acoustic waveguides as well as seamless integration with electromagnetic solvers for optomechanical device design. The accuracy of the solver is demonstrated by calculating eigenfrequencies and mode shapes for common acoustic modes in several simple geometries and comparing the results to analytical solutions where available or to numerical solvers based on more computationally expensive methods.


I. INTRODUCTION
Recent advances in several fields have attracted growing interest in the design of chip-scale acoustic devices that can interface with electrical and optical integrated components. Optomechanics is a prime example, where interacting acoustic and optical fields enable novel functionalities, such as ultra-sensitive quantum measurements 1,2 , narrowlinewidth lasers 3,4 , optomechanical memory 5,6 , non-reciprocity and optical diodes [7][8][9] , optical cooling 10-12 , phononic topological insulators 13,14 , optical amplifiers 15 , improved gravity wave detection 16,17 , microwave filters 18 , and quantum state transfer [19][20][21] . The field of RF micro-electromechanical systems is another important example where electroacoustic transduction of bulk and surface acoustic waves in acoustic resonators enables some devices which outperform conventional RF electronics. These include reconfigurable filters 22 , narrowband signal filtering 23 , and high quality factor (Q) resonators 24 . These devices are also being integrated into microelectronic systems for improved performance 25 . For all these acoustic wave based devices, good performance requires confining the acoustic energy to a small crosssectional area (waveguides) or volume (resonators), phase-matching the acoustic wave to transducer arrays and/or optical waves, and optimizing transduction efficiency. Numerical tools for designing and simulating acoustic waveguide modes are thus necessary to enable efficient device designs, intricate nanoscale coupling schemes, and novel device architectures.
Previous work has predominantly focused on full-wave simulation of 2D and 3D domains due to the importance of these problems in geophysics. 3D solvers have been developed for anisotropic, heterogeneous domains using both finitedifference 26 and finite-element 27 methods and are currently the predominant method for designing acoustic devices in the GHz frequency range. While many commercial software tools 28 allow the design of acoustic waveguides using a full 3D solver, sometimes more efficiently by reducing the volume using Floquet (periodic) boundary conditions, the most efficient approach is the maximally reduced, 2D problem formulation (disregarding cross-section symmetries). The 2D formulation returns an orthogonal, complete set of modes (field and frequency) at a specified propagation constant. In contrast, the solution of a 3D volume returns all resonant modes including those in higher order Brillouin zones, which are an artifact of the 3D formulation and usually undesired. This rigorous reduction of a 3D problem to a 2D simulation domain is referred to as a 2+1D simulation. Papers utilizing a 2+1D version of the finite-element method have been implemented to solve free, isotropic waveguide geometries in the ultrasound regime 29 and extended to axially symmetric waveguides 30 , embedded waveguides 31 , and viscoelastic materials 32 . It should be noted that generic FEM tools such as COMSOL, which provide an interface to solve arbitrary partial differential equations, are in principle capable of solving the 2+1D Cartesian acoustic waveguide simulation (by entering custom equations).
In this work, we demonstrate an acoustic waveguide mode solver based on the finite-difference method, analogous to electromagnetic (EM) mode solvers, which solves the linear isotropic elastic wave equation. The mode solver combines all capabilities of FEM mode solvers with the efficiency of finite-differencing, all physical boundary conditions as well as perfectly matched layers, directly overlaps with the optical Yee grid 33 , and is verified for acoustic frequencies ranging from 1 MHz to 10 GHz. We utilize a staggered-grid discretization that preserves second-order accuracy for all physical quantities of interest, by analogy with the Yee scheme first used in EM 33 and identical to the staggered grid of 34,35 . The solver finds the acoustic eigenmodes of a structure that is invariant along one Cartesian dimension. The source-free vectorial elastic wave equation in an inhomogeneous, isotropic, linear medium is formulated as an eigenvalue equation. Given a specific material configuration and the propagation constant, a unique set of eigenmodes with corresponding modal frequencies can be found. The acoustic problem is only a linear eigenvalue problem when formulated such that the propagation constant β is given and the eigenfrequency ω is the eigenvalue. This is notably different than the analogous electromagnetic problem which, due to Gauss's Law, can be formulated as a linear eigenvalue problem with either ω or β given and the other variable as the eigenvalue. A finite-difference method is used to sample the inhomogeneous medium and acoustic field over the computational domain, the cross-section of the structure. This operation transforms the continuous eigenvalue problem into a sparse matrix which can be solved by standard numerical sparse matrix eigen solver methods (e.g. shifted inverse method via Arnoldi iteration and a sparse linear solver). An example case is the 'eigs' function in MATLAB 36 . The solver calculates, to second order accuracy on the finite-difference grid, the acoustic eigenmodes of a straight acoustic waveguide with specified cross-section and propagation constant. The use of 2+1D mode solvers in both optical and acoustic domains allows for accurate and efficient calculation of propagation parameters and coupling terms, which can be input to the acoustooptic simulations for accuracy comparable to full 3D simulations. To this end, we have made our mode solver code, in MATLAB implementation, freely available 37 .
We present the theoretical and mathematical basis for the elastic wave equation eigenmode decomposition in Sec. II. This section also includes implementation details such as the finite-differencing scheme and boundary conditions. We then implement the method in MATLAB and validate its accuracy by analysis of the matrix construction, convergence tests, example cross-sections with analytical solutions, and other examples which can be solved numerically in 3D using a commercial FEM solver 28 in Sec. III.

II. THEORY, MATHEMATICAL FRAMEWORK, AND IMPLEMENTATION OF THE MODE SOLVER
The mathematical basis for the eigenmode solver and its implementation are described in this section. The sourcefree, linear, isotropic elastic wave equation is cast as an eigenmode problem with frequency eigenvalues and discretized. The only assumptions made are that the cross-section is invariant along one linear direction,ẑ; all materials are linear, isotropic, and time-invariant; and that a computational domain that is finite in cross-section represents well the modes of the structure (justified for 'confined' modes, such as in electromagnetics 38 ). The isotropic assumption is made for convenience and is not essential; the method can be applied to anisotropic media. This approach to solving for waveguide modes is analogous to electromagnetic mode solvers, where instead of an electric or magnetic field we solve for the elastic displacement field. Notably, in the elastic equation, there is no equivalent of Gauss's Law and thus we solve for three field components and have three polarization families. In EM, Gauss's Law reduces the eigen problem to two polarization mode families and specification of two field components fully defines the mode (e.g. E x , E y ).

A. Derivation of Isotropic Linear Elastic Wave Equation
We begin with Newton's 2 nd law written in a density formulation, the strain-displacement relation, and generalized Hooke's law where ρ(r, t) ≡ ρ(r) is the spatial distribution of material density, u(r, t) is the elastic displacement field, ∇· is the tensor divergence,σ(r, t) is the stress tensor,¯ (r, t) is the strain tensor, ∇ s is the symmetric spatial vector gradient, C(r, t) ≡ C(r) is the fourth order stiffness tensor, and : denotes a tensor inner product. Note that while stress and strain are second order tensors, due to symmetry considerations they can be unwrapped as six-vectors following the Voigt notation (σ = [σ xx σ yy σ zz σ yz σ xz σ xy ] T ) 39 , which is the form used in this paper. The strain tensor has an analogous form. The operator ∇ s , which acts on a vector to give a rank 2 tensor, is the adjoint of the tensor divergence operator (−∇·). C is a symmetric operator (even in the presence of loss) due to our assumption of linear, time-invariant materials 39 , and we use the notation ∂ i to refer to the partial derivative ∂/∂i. Substituting Eqs.
(2)-(3) into Eq. (1) in order to factor out the stress tensor yields the linear elastic wave equation written in terms of the displacement field We define a weighted displacement fieldũ ≡ √ ρu. This allows Eq. (4), a generalized eigenvalue problem, to be recast as an ordinary eigenvalue problem with a Hermitian operator in the absence of loss. Invariance of ρ and C with time means that the system has a spectrum defined by a linear eigenvalue problem by setting ∂ t → jω. The elastic wave equation can be written as This has the form of an eigenvalue equation (Λx =Hx), so we identify the eigenvalue as Λ ≡ ω 2 , the eigenvector as x ≡ũ, and the symmetrized elastic resonance operatorH as To clarify this expression, the modified differential operators (in the form that operate on the reduced, six-vector notation used forσ and¯ ) can be written in matrix form To this point, we have defined the acoustic resonator problem (in 3D). Next, when solving for eigenmodes of a structure with z-invariant geometry, by Fourier transformation along z the z-directed derivative becomes the propagation constant (∂ z = −jβ). This reduces the problem to one on the cross-sectional plane and ensures that only modes with the specified propagation constant will be returned by the solver.
Next, making the isotropic assumption, the stiffness tensor (in the form that operates on the six-vector notation) can be reduced to where λ(r, t) ≡ λ(r) and µ(r, t) ≡ µ(r) are the first and second Lamé parameters.

B. Discretization Scheme
In order to solve the eigenvalue problem in an arbitrary cross-section geometry (ρ, C) numerically, Eq. (5) is discretized to arrive at a form with a finite number of degrees of freedom, and is cast as a matrix eigenvalue problem. An appropriate 3D grid is first formulated which accurately captures the physics to second-order accuracy in the discretization and the grid is then collapsed to the 2D simulation domain. The collapse is performed such that the cross-sectional locations of all grid points are unchanged while the z-coordinates of all grid points are set to a single value. The zderivatives, which would be performed between two points in the 3D case, amount to multiplying a single grid point by −jβ in the 2D case. More formally, for a 3D grid with z-invariance, ∂ z → −jβsin(β∆z/2)/(β∆z/2) exp(−jβ∆z/2) 40 , but as ∆z → 0 then ∂ z → −jβ.
The discretization grid used is depicted in Fig. 1. This grid is based on physical intuition from solid mechanics: the state of an element cube is primarily described by the principal stresses (σ xx , σ yy , σ zz ), which we define to reside at the center of the cube. The principal stresses lead to the deformation of the faces of the cube, such that the center of each face is the sampling location of the associated normal displacement. If grid point (i, j, k) is associated with the principal stresses (and the center of the cubic volume elements), the corresponding displacements are u x : (i + 1/2, j, k), u y : (i, j + 1/2, k), and u z : (i, j, k + 1/2). An appropriate grid for the shear stresses is found to be a further half-step offset, such that shear stresses are located at σ xy : (i + 1/2, j + 1/2, k), σ xz : (i + 1/2, j, k + 1/2), and σ yz : (i, j + 1/2, k + 1/2). The strain grid is co-located with the stress grid. The finite-difference operators that approximate the spatial partial derivatives transform quantities from the stress/strain grid to the displacement grid and back, as expected and desired. Because isotropic materials are assumed, the material operator C does not induce a change of grid coordinates (the strain at a specific grid point is only related to the stresses at the same grid point). This method can also be extended to anisotropic materials at the cost of additional complexity in the C matrix, which must then have different coefficients sampled on different grids.
The staggered-grid scheme used here is analogous to the Yee grid 33 , commonly used in EM solvers, which preserves second-order accuracy for all fields due to centered differencing (when all materials/coefficients vary spatially on the scale of the discretization). Our discretization grid is slightly more complicated in that the elastic displacement field and shear stress directly replicate a Yee grid, while the principal stresses occupy an additional position at the center of each 'cube'.

C. Finite Difference Operators
Centered differences in 2 nd order differential equations can be formed by appropriate combinations of forward and backward differences on a staggered grid 40 . Denoting forward differences as∂ i and backward differences as∂ i , the differential operators can be rewritten in terms of forward and backward differences (where the propagation constant has been substituted for ∂ z ) as It should be noted that, because∂ i † = −∂ i , ∇ s and −∇· remain adjoints in discrete form (∇ † s = −∇·).

D. Matrix Operator Construction
Referring to Eq. (6), the resonance operator can be rewritten as a matrix. This matrix can be found by discretizing the displacement field and unwrapping the field into a single column vector. A convenient ordering was found as Thus, each vector component is unwrapped in the x − y plane along the x-dimension first and the components are then concatenated. This results in a vector of approximate length 3n x n y from an n x × n y × 3 matrix. A similar method is used to unwrap the stress tensor in a vector of approximate length 6n x n y of the form with the strain tensor having an equivalent form. Note that the different components σ ij and u i are of different lengths since they are sampled at different locations. Shown in Table I  Component ux uy uz σxx σyy σzz σyz σxz σxy n comp x nx + 1 nx nx nx nx nx nx nx + 1 nx + 1 n comp y ny ny + 1 ny ny ny ny ny + 1 ny ny + 1 The differential operator hereby becomes a matrix operator where each entry in Eqs. (9), (10), (11) is now a block matrix which relates one component to another. After populating the ∇·, C, and ∇ s matrices with the block matrices, they are multiplied together to givē where we have shortened λ + 2µ to the Voigt notation c 11 39 . We here note that the matrix, currently in a Hermitian form in the absence of loss, can be cast as a real symmetric matrixH (a special case of Hermitian matrices) by defining a modified eigenvector formũ ≡ [ũ x ,ũ y , −jũ z ], i.e. we expect that the u z component will be in quadrature with the u x , u y components. Hermitian matrices have real eigenvalues (energy conservation) and their eigenmodes form a complete, orthogonal set. Casting the operator matrix into a symmetric form, which has eigenvectors that can be chosen to be entirely real, additionally implies that u z is guaranteed to be in quadrature with u x and u y . If loss is present, then this modified matrix form will be complex symmetric rather than real symmetric and we would expect to obtain complex eigenvalues and non-orthogonal eigenmodes. Defining the modified operator matrixH such that ω 2ũ =H ũ , we can calculate this modified operator matrix by recognizing thatũ =Rũ andH =RHR −1 whereR is a diagonal matrix with diag(R) = [1, 1, −j]. We can then write the modified operator matrix as We remind the reader here that, even for complex finite difference operators,∂ T i = −∂ i so that this modified matrix H is symmetric. It is additionally real if both material parameters (λ, µ) and finite difference operators (∂ i ,∂ i ) are real, corresponding to no material loss and no complex coordinate stretching radiation (absorbing boundaries), respectively.

E. Modal Orthogonality
The orthogonality condition for vectors (in the sense of 1D tensors) is represents a 3-vector field in a 2D domain. Noting that thatũ * i ·ũ j = ρu * i · u j , we can then replace x →ũ and write the orthogonality condition and normalization as Noting that the particle velocity field is v i ≡ ∂ t u i = jω i u i , we can substitute u i = −jv i /ω i into Eqs. (17)-(18) to obtain measures of the kinetic energy E K = ρ|v| 2 /2 . Alternatively, because the potential energy V =¯ † C¯ /2 is equivalent to the kinetic energy when averaged 39 , this gives the orthogonality condition that would have been obtained if the wave equation were formulated with¯ as the free variable. This is analogous to the electromagnetic case, where resonator problems use an energy formulation to determine modal orthogonality 41 .

F. Boundary Conditions
Four boundary conditions are implemented: fixed boundary (u| Bound = 0), free boundary (σ ·n = 0), symmetry boundary (∂ n u n = 0, u t | Bound = 0), and anti-symmetry boundary (u n | Bound = 0, ∂ n u t = 0). Perfectly matched layers (PMLs), radiation-absorbing regions which permit simulation of radiating structures 42 , are also implemented. For all boundary conditions, we locate the boundary such that the boundary-normal displacements (u n ), corresponding out-of-plane shear stresses (σ nz ), and in-plane shear stresses (σ xy ) lie on the boundary. The choice of boundary location is depicted in Fig. 1(d) as the location of the red line. We additionally allow for vacuum (no material) in the simulation domain.
The fixed boundary condition requires the displacement field be set to zero on the boundary. An efficient and computationally stable implementation is to remove the corresponding boundary-normal displacements u n (since they are the only ones coincident with the boundary) from the matrix operator and solution vector, which is the method used here. We validate this boundary condition in Sec. III C by confirming that a beam with fixed boundary conditions returns the correct eigenmode shapes and frequencies (Fig. 4).
For the free surface boundary condition, the boundary-normal stress components must go to zero at the boundary. In a similar manner to the implementation of the fixed boundary condition, we remove the shear stresses σ nz , σ xy on the boundary. This boundary condition is also validated in Sec. III C by confirming that the eigenmode shapes and frequencies are accurately calculated (Fig. 4).
Symmetry and anti-symmetry boundary conditions are hybrids of free and fixed boundary conditions. Similar to the case in electromagnetics, symmetry of the normal displacement corresponds to anti-symmetry of the transverse displacement, while anti-symmetry of the normal displacement corresponds to symmetry of the transverse displacement. A symmetry (anti-symmetry) boundary therefore refers to symmetry (anti-symmetry) of the normal displacement and anti-symmetry (symmetry) of the transverse displacement.
The symmetry and anti-symmetry boundary conditions are implemented by combinations of the free and fixed boundary conditions. The proper implementation of these boundary conditions is validated in Fig. 5 by comparison of the free beam modes calculated with and without the symmetry/anti-symmetry boundary conditions. We enable vacuum within the simulation domain by removing both stress and displacement grid points within the vacuum region. We also remove the stress sampled on the vacuum-material boundary and double boundary-crossing derivatives as appropriate. This exactly emulates a free surface boundary condition along the interface, as desired.
We additionally implement PMLs in order to support simulation of leaky modes. The PML is not formally a boundary condition and is instead a region within the computational domain [42][43][44][45] . A desired use for this acoustic mode solver is to estimate propagation losses of a particular acoustic mode in a waveguide, and a key loss mechanism can be radiation loss (analogous to a distributed version of clamping loss in resonators, e.g. cantilevers). We refer to acoustic modes exhibiting radiation loss as a leaky modes, which can still be considered confined modes when the radiation loss rate is much smaller than the oscillation frequency.
The PML boundary condition can be implemented as a complex coordinate stretching near the boundary 45 . This can be achieved by making the grid discretization components ∆x and ∆y complex within the PML region. We choose to use a fixed boundary condition adjacent to the PML in order to force the acoustic field to go to zero at the boundary. In order to replicate the PML used by the commercial solver used for comparison, a linear coordinate stretching was used of the form where α is the scaling parameter of the PML. To demonstrate the accuracy of the PMLs implemented in this solver, we constructed a leaky acoustic waveguide that tunnels radiation across a barrier into the continuum on both sides (Fig. 6).

III. MODE SOLVER VALIDATION AND DEMONSTRATION
In order to demonstrate the accuracy of the proposed acoustic mode solver, we examine the resonance matrix to verify its Hermitian nature in Sec. III A, then compare several simple cases in our MATLAB implementation to both analytical solutions (Sec. III B) and a commercial, numerical 3D solver (Sec. III C). We verify the symmetry and anti-symmetry conditions in Sec. III D and demonstrate support of leaky modes in Sec. III E.

A. Mathematical Validation: Matrix Properties
To test our discretization scheme, we examine the properties of resonance operator matrices generated by the mode solver code. To provide an example of the resonance operator matrix, we implement a trial cross-section with two free and two fixed boundary conditions on a 4 × 4 pixel grid in Fig. 3(a). We additionally verify in Fig. 3(b) second order scaling of the error with discretization. For this case, we chose an anti-symmetric Lamb wave with a 10 µm wavelength in a 1 µm thick free slab and varied the discretization along the slab thickness direction. The calculated eigenfrequency and mode shape errors are relative to the analytical solutions, described further in Sec. III B. Our convergence plots indicate that the mode solver eigenfrequency error scales as ε ∝ (∆y) 2 , i.e. that our mode solver is second-order accurate as expected. The mode shape error scales as ε ∝ (∆y) 4 due to the squaring in the mode shape error formula. We additionally expect that, in the absence of material loss or PMLs, the resonance operator would be real and symmetric even after discretization which we verified using a 200 × 200 pixel grid simulation with both free and fixed boundary conditions. We also confirmed that the operator is complex symmetric when PMLs are used. Other hallmarks of Hermitian operators are real eigenvalues, orthogonal eigenvectors, and that the eigenvectors form a complete basis. We have shown an example demonstrating the orthogonality of the eigenvectors in both a free beam and a waveguide without PMLs, and lack thereof when PMLs are implemented, in Fig. 3(c).

B. Analytical Validation: Free Slab Waveguide
We choose a slab waveguide, which can be analytically solved, to compare the proposed mode solver with elastic wave theory. The simulation implementation and results are shown in Fig. 2. We investigate three fundamental acoustic waves in a thin plate [ Fig. 2(b)], namely pure shear waves, anti-symmetric Lamb waves, and symmetric Lamb waves 39 . We calculate the theoretical dispersion relations [ Fig. 2(c)] of these waves compared with the mode solver solutions as the wavelength is varied from 100 µm to 0.5 µm to achieve full coverage of the thin plate regime and into the thick plate (infinite half-space) regime. For the mode solver solutions, we use a plate which is much wider than its thickness (1 × 1000 µm) with the appropriate symmetry boundary conditions on the edges to approximate an infinitely long (1D cross-section) plate. We then calculate the error of the mode solver relative to theory, both in terms of mode shape and eigenfrequency, in Fig. 2(d). We sample the acoustic mode distribution at the center of the plate as an approximation to the theoretical 1+1D case. The mode shape error used here is formulated as where U o is the true mode shape and U M S is the mode shape found by the solver. For the slab case, we integrate over only a single dimension as the mode shapes are invariant along the second cross-sectional dimension. For the general case, the integration is carried out over the entire cross-section. The mode solver reproduces the key aspects of these acoustic waves with < 1% error across a range of wavelengths which captures the thick, thin, and wavelength-scale plate thickness regimes, indicating that the mode solver faithfully models linear elastic physics.

C. Numerical Validation: Free and Fixed Beam Waveguides
We now consider a rectangular beam waveguide (in both suspended and fixed, i.e. encased in an infinitely rigid medium, configurations) to compare the proposed mode solver with numerical solutions obtained from a commercial 3D solver 28 . We investigate the lowest order mode in each of four basic mode families: vertically (y) polarized shear waves, horizontally (x) polarized shear waves, torsional waves, and pressure waves. We choose as our cross-section a 2 × 1 µm silicon beam so as to avoid a degeneracy of the x-and y-shear modes. The same geometry is implemented in the commercial 3D solver with the beam length chosen equal to the simulated wavelength and Floquet periodic boundary conditions applied to the z-oriented boundaries. Both solvers are then used to find these four acoustic modes and corresponding eigenfrequencies across a wavelength range of 100 µm to 0.5 µm. These comparisons are made in Fig. 4 where we depict the comparative error of the mode solver relative to the 3D solver in terms of frequency and modeshape. For the majority of wavelengths, both the mode shape and frequency are very close (< 1%) to the 3D solver results which are defined here as the reference value.

D. Symmetry/Anti-symmetry Boundary Condition Validation
Having demonstrated the accuracy of the solver, we can use the free beam modes as a basis for comparison with modes computed in a reduced simulation domain that takes advantage of geometrical symmetries. For this case we halve (or quarter) the simulation domain and use a symmetry/anti-symmetry boundary condition on the cut plane to recreate the same mode. We have chosen to use both boundary conditions to recreate the same mode as a redundant test of the boundary conditions, and the pressure mode is chosen as a demonstration that two orthogonal boundary conditions can be used at the same time without loss of accuracy. This comparison is shown in Fig. 5, where examples of the modes and symmetry conditions are shown in Fig. 5(a) and the errors are shown in Fig. 5(b). The jump in error in both eigenfrequency and mode shape at small wavelengths is due to the simulated modes being nearly degenerate in frequency at these wavelengths, causing the solver to return mixed mode shapes and eigenfrequencies when the full simulation domain is used. The halved/quartered simulation domains did not have mixing due to the imposed symmetry requirements.

E. Solver Demonstration: Evanescent Waveguiding, Leaky Modes, and PMLs
We now construct an acoustic waveguide formed of two materials (a 'core' and a 'cladding') and find evanescently confined acoustic modes. Adding a radiation layer (formed of core material) to the waveguide cross-section creates a leaky acoustic waveguide that tunnels radiation across a (cladding) barrier into the radiation-mode continuum on both sides (see analogous optical slab waveguide 46 ). Implementation of PMLs on the simulation domain edges then allows outward radiating waves to be absorbed without reflection, enabling the simulation of leaky modes analogous to electromagnetic solvers 47 and calculation of waveguide mode radiation losses. The configuration, modal amplitude plots, radiation patterns, and calculated radiation losses are shown in Fig. 6. First, an evanescently guiding waveguide is constructed [ Fig. 6(a)] by embedding a strip of silica within a silicon substrate. This waveguide supports two 'polarizations' of acoustic modes, analogous to the Love and Rayleigh surface waves 39 . The former wave is chosen and the depicted log-scale plot of the horizontal displacement amplitude demonstrates that the wave is evanescently confined. We then add silica slabs to either side of the silica waveguide, separated by variable width silicon barriers, which cause the waveguide and guided acoustic modes to couple to the slab radiation continuum on both sides and become 'leaky' [Fig. 6(b)]. Oscillation of the field within the silica layer can be seen in the inset, and non-decaying intensity within the leaky region can be seen in the log-scale plot; both are indicative of the radiative loss of the guided mode into the adjacent silica layers.
An extended simulation domain with a longer silica layer further demonstrates the presence of this mode leakage in Fig. 6(c) where a cross-sectional slice plots both horizontal and vertical displacement amplitudes within the different material regions. The mode displacement is predominantly confined in the waveguide with an evanescent tail in the silicon barrier. Some of the acoustic energy tunnels across the barrier into the silica slab, where it becomes oscillatory, and then enters the PML and is attenuated as expected. A comparative plot of this propagation loss, as a function of silicon barrier width, is shown in Fig. 6(d) as numerical validation of the accuracy of the PML. By taking a Fourier transform of the mode leakage section we show that two primary radiation field components are excited in the leaky wave [ Fig. 6(e)], corresponding to quasi-Love and quasi-Rayleigh waves in the silica slab. Fig. 6(f) shows a cross-section in the x − z plane depicting the horizontal (left) and vertical (right) displacement amplitudes, where the wavefronts have been outlined for both polarizations and overlaid on both plots (gray for horizontal displacement, black for vertical).

IV. CONCLUSION
We have developed an eigenmode solver based on a finite-difference scheme on a staggered-grid for the linear, isotropic elastic wave equation. We have verified the accuracy of the solver by comparison with both theory and a 3D numerical solver (COMSOL). We have also used these comparisons to demonstrate correct implementation of fixed and free boundary conditions. We then used a leaky acoustic waveguide to demonstrate the solver's ability to simulate evanescent guiding, leaky modes, and numerically accurate PMLs. We expect that the mode solver will be of utility in the design of on-chip acoustic wave devices. It is particularly suited to interfacing with electromagnetic solvers on the Yee grid for applications based on acousto-optics and optomechanics.