Acoustic Waveguide Eigenmode Solver Based on a Staggered-Grid Finite-Difference Method

A numerical method of solving for the elastic wave eigenmodes in acoustic waveguides of arbitrary 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 (field and frequency) of the vector-field elastic wave equation with a given propagation constant. 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 that in 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 across four orders of magnitude in frequency in several simple geometries and comparing the results to analytical solutions where available or to numerical solvers based on more computationally expensive methods. The solver is utilized to demonstrate a novel type of leaky-guided acoustic wave that couples simultaneously to two independent radiation channels (directions) with different polarizations – a ‘bi-leaky’ mode.

In contrast to the 2D case, solution of this 'resonator' eigen problem on a 3D volume returns all resonant modes, including those in higher order Brillouin zones, which act as unwanted artifacts of the 3D formulation of the 2D waveguide problem. A 2 + 1D FEM has been previously implemented to solve waveguide geometries in the ultrasound regime (referred to as the Semi-Analytic Finite Element, or SAFE, method) [29][30][31][32] .
Although FEM is more sophisticated and general than FDM (which it includes as a subset), FDM on a uniform grid has a number of strengths when it comes to the design of nano-scale photonic, and we believe by extension also phononic, devices and circuits. While optical FDM solvers are dominant in academic and commercial tools, an FDM solver has not been previously demonstrated for acoustic 2 + 1D simulations. First, a key design parameter in coupled-waveguide and coupled-cavity circuits is frequency or propagation constant splitting (supermode frequencies/propagation constants), which measures the strength of coupling of two structures. A uniform-grid FDM can provide these frequency/propagation constant differences between resonant/guided modes of coupled cavities/waveguides accurately with a coarse grid because absolute errors are common-moded out, whereas irregular meshes used in a FEM can introduce 'grid induced' or numerical fictitious detunings between nominally degenerate elements. These fictitious detunings result in physically incorrect coupling dynamics, and thus require a much finer grid (to achieve high absolute accuracy). Second, design of coupled structures often involves generating design curves by sweeping dimensions (e.g. spacing between coupled components). In these situations an unstructured FEM grid does not sweep consistently because the mesh adapts to the geometry, and may result in a non-smooth result vs. sweep parameter. An FDM implementation on a uniform grid can provide smooth, physically accurate results with proper grid alignment. In addition, FDM is simple and the implementation compact, allowing for straightforward addition of new functionalities such as other coordinate systems (e.g. cylindrical or elliptical) and addition of support for anisotropic materials. Last, the staggered grid FDM we implement naturally mates to electromagnetic solvers implemented on a Yee grid 33 for optomechanical interaction simulations.
In this work, we develop and implement (and make the implementation freely available 34 ) an acoustic waveguide mode solver which solves the linear isotropic elastic wave equation based on FDM (analogous to the SAFE method) which is comparably accurate to 3D FEM and more computationally efficient. We implement all common boundary conditions including free, fixed, and symmetry/anti-symmetry, as well as perfectly matched layers (PMLs) which permit efficient modeling of outward radiation and leaky modes. We validate the implementation's accuracy from 1 MHz to 10 GHz to demonstrate its applicability in higher frequency regimes. A finite-difference method with staggered-grid discretization (identical to the staggered grid used in certain 3D acoustic solvers 35,36 ) is chosen to preserve second-order accuracy in primary and derivative physical variables and enable direct alignment with the optical Yee grid 33 for coupling to electromagnetics.
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. This β to ω 'resonator' formulation is the 'band structure' form, convenient for computing the band diagram, with complex ω representing loss. The interchanged case with β as the eigenvalue is, however, a quadratic eigenvalue problem. 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. The acoustic ω to β 'waveguide' formulation can be linearized 37,38 , but this second type of implementation is beyond the scope of this work. A finite-difference discretization approach is then applied and transforms the continuous eigenvalue problem into a sparse matrix which can be solved by standard numerical sparse matrix eigen solver methods (here, the 'eigs' function in MATLAB 39 ).
Notably, 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 used as inputs to coupled mode theory or acousto-optic simulations 40,41 for accuracy comparable to full 3D simulations but with orders of magnitude faster simulation times. Such fast simulation enables design approaches that are not practical with full 3D simulations, such as large parameter sweeps and optimizations. It also provides greater insight into the physics, and can be a building block enabling other physically insightful simulation techniques to be borrowed from optics, such as the film matching method 42,43 .

Theory, Mathematical Framework, and Implementation of the Mode Solver
The mathematical basis for the eigenmode solver and its implementation are described in this section. We assume that the structure's cross-section is invariant along one linear direction, ẑ, and that all materials are linear, isotropic, and time-invariant. 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, but instead of an electric or magnetic field we solve for the elastic displacement field. Notably, in the elastic equations, there is no equivalent of Gauss's Law and thus we solve for three field components and have three polarization families. In electromagnetics, Gauss's Law reduces the eigen problem to two polarization mode families and specification of two field polarization components fully specifies the mode (e.g. E x , E y ).

Derivation of Isotropic Linear Elastic Wave Equation.
We begin with Newton's 2nd law written in a density formulation, the strain-displacement relation, and generalized Hooke's law, respectively: where ρ(r, t) ≡ ρ(r) is the spatial distribution of material mass density, u(r, t) is the elastic displacement field, ∇⋅ is the tensor divergence, σ t r ( , ) is the stress tensor, t r ( , ) ε is the strain tensor, ∇ s is the symmetric spatial vector gradient, is the fourth order stiffness tensor, and : denotes a tensor inner product 44 . Note that while stress and strain are second order tensors, due to the underlying symmetries they can be unwrapped as six-vectors following the Voigt notation (σ σ σ σ σ σ σ = [ ] xx yy zz yz xz xy T ) 44 , which is the form used in this paper. The strain tensor has an analogous form. The operator ∇ s is the adjoint of the tensor divergence operator ( † ∇ = −∇⋅ s ).  is a symmetric operator (even in the presence of loss) due to our assumption of linear, time-invariant materials 44 , and we use the compact notation ∂ i to refer to the partial derivative ∂/∂i.
Substituting Eqs (2) and (3) into Eq. (1) in order to factor out the stress and strain tensors yields the linear elastic wave equation written in terms of the displacement field ρ∂ = ∇ ⋅ ∇ . u u : (4) t s 2  Invariance of ρ and  with time means that the system has a spectrum defined by a linear eigenvalue problem by setting ∂ t → jω. We define a weighted displacement field ρ ≡ u u  . This allows Eq. (4), a generalized eigenvalue problem with eigenvalue ω 2 when ∂ t → jω, to be recast as an ordinary eigenvalue problem with a Hermitian operator in the absence of loss. 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 u ≡ , and the symmetric elastic resonance operator To clarify this expression, the differential operators (in the form that operate on the reduced, six-vector notation used for σ and ε ) can be written in matrix form as 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 yields 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 c ij (r, t) ≡ c ij (r) are the Voigt notation stiffness coefficients which can be defined in terms of the Lamé parameters 44 as Discretization Scheme. In order to solve the eigenvalue problem in an arbitrary cross-section geometry (ρ, ) numerically, Eq. (5) is discretized over a finite 2D domain terminated by boundary conditions 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 but are all referenced to a single z position. The z-derivatives, which would be performed between two points in the 3D case, amount to multiplying a single grid point by −jβ in the 2D case. Formally, for a z-invariant structure on a 3D grid, , 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 [ Fig. 1(b)] 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 surface-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 (ε mn is co-located with σ mn ) and the material grid (, ρ) is located at the cube centers (i, j, k). 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 Scientific REpoRtS | 7: 17509 | DOI:10.1038/s41598-017-17511-x  does not induce a change of grid coordinates (the stress at a specific grid point is only related to the strains at the same grid point). This method can also be extended to anisotropic materials at the cost of additional complexity in the  matrix.
The staggered-grid scheme used here is analogous to the Yee grid 33 , commonly used in electromagnetic 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' .
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 45 . 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 ) aŝ^~j  Matrix Operator Construction. The resonance operator can be rewritten as a matrix, found by discretizing the displacement field and unwrapping the field into a single column vector. A convenient ordering is where {u i } is a vectorized form that is unwrapped such that adjacent points in the x-dimension are adjacent in the vector. 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 σ σ σ σ σ σ σ = [{ }{ }{ }{ }{ }{ }] xx yy zz yz xz xy T 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.
The differential operator hereby becomes a matrix operator where each entry in Eqs (7), (8), (9) is now a block matrix which relates one component to another. After populating the ∇⋅, , and ∇ s matrices with the block matrices, they are multiplied together to give the matrix operator H . We here note that the matrix, currently in a complex-valued Hermitian form in the absence of loss (due to jβ terms), can be cast as a real symmetric matrix ′ H (a special case of Hermitian matrices) by defining a modified eigenvector form u u ju u [ , , ] x y z    ′ ≡ − , i.e. the u z component will be in quadrature with the u x , u y components. Hermitian matrices have real eigenvalues (energy conservation) and their eigenvectors form a complete, orthogonal set. Symmetric matrices, as an additional trait, have eigenvectors that can be chosen to be entirely real. 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 eigenmodes with an unconjugated orthogonality condition that is not defined by energy conservation. Defining the modified operator matrix H ′ such that ω ′ = ′ ′   H u u 2 , we can calculate this modified operator matrix by recogniz- − . We can then write the modified operator matrix as^^~~^^~~~^ρ We remind the reader here that, even for complex finite difference operators and complex discretization Δx, ^∂ = −∂ i T i so that this modified matrix H ′ is symmetric. It is additionally real if both material parameters (c ij ) and finite difference operators ( , ) are real, corresponding to no material loss and no complex coordinate stretching (radiation absorbing boundaries), respectively.

Modal Orthogonality. The orthogonality condition for vectors (in the sense of 1D tensors) is
∬ ⁎ 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 (11) and (12) to obtain measures of the kinetic energy ∬ ( ) Alternatively, because the potential energy is equivalent to the kinetic energy when averaged 44 , 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 46 . 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) are radiation-absorbing regions which permit simulation of radiating structures 47 and 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. 2.3 by confirming that a beam with fixed boundary conditions returns the correct eigenmode shapes and frequencies.
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. 2.3 by confirming that the eigenmode shapes and frequencies are accurately calculated.
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 proper implementation of these boundary conditions is validated in Sec. 2.3 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 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 [47][48][49][50] . 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 50 . This can be achieved by making the grid discretization components Δx and Δy complex within the PML region. 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 in Sec. 2.5.

Mode Solver Validation and Demonstration
Mathematical Validation: Matrix Properties. To test our discretization scheme, we examine the properties of resonance operator matrices generated by the mode solver code. We implement a trial cross-section with two free and two fixed boundary conditions on a 4 × 4 pixel grid to allow the matrix structure to be visualized in Fig. 2(a). We additionally verify in Fig. 2(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. 2.2. 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 verified that, in the absence of material loss or PMLs, the resonance operator is real and symmetric after discretization using a 200 × 200 pixel simulation. 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, as well as the modified orthogonality when PMLs are implemented Fig. 2(c).
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. 3. We investigate three fundamental acoustic waves in a thin plate [ Fig. 3(b)], namely pure shear waves, anti-symmetric Lamb waves, and symmetric Lamb waves 44 . We calculate the theoretical dispersion relations [ Fig. 3(c)] of these waves compared with the mode solver solutions as the wavelength is varied from 100 μm to 0.5 μm to span 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) and sample the acoustic mode distribution at the center of the plate to approximate the theoretical 1 + 1D case. We then calculate the error of the mode solver relative to theory, both in terms of mode shape and eigenfrequency, in Fig. 3(d). The mode shape error used here is formulated as where U o is the true mode shape and U MS 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 <0.01% 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.
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 (which returns a considerably larger solution space and is correspondingly more computationally intensive). 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 frequencies 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 mode shape. For the majority of wavelengths, both the mode shape and frequency are very close (<1%) to the 3D solver results, defined as the reference value. The increased error at large wavelengths in the fixed beam, and lack of data at 100 μm, is due to difficulty in simulating these modes in the commercial solver.

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.
Solver Demonstration: Evanescent Waveguiding, Leaky Modes, and PMLs. We now construct an acoustic waveguide 7,51,52 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 waveguide pg. 93 53 ). 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 54 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 44 . 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 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 PML accuracy. 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 (u x ) and quasi-Rayleigh (u y ) waves in the silica slab. Figure 6(f) shows a cross-section in the x − z plane depicting the u x (left) and u y (right) displacement amplitudes, where the wavefronts have been outlined for both polarizations and overlaid on both plots (gray for horizontal displacement, black for vertical). The figure shows a waveguide mode that has two leakage channels -quasi-Love and quasi-Rayleigh -with different radiation angles. Multiple discrete radiative channels are not present in electromagnetic leaky modes. We suggest naming this type of mode 'bi-leaky' . It may find utility in design of phononic circuits (e.g. filters 55 ).

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. Potential future extensions to this solver include support for material anisotropy, a cylindrical coordinate system formulation for simulating ring and disk resonators, co-integration with an optical mode solver for optomechanics applications, and reformulation of the eigenvalue problem to solve for β as the eigenvalue with ω provided, similar to electromagnetic mode solvers for waveguides.
Data Availability. The datasets generated and analyzed during the current study are available from the corresponding author on reasonable request.