# Elastic Modes of an Anisotropic Ridge Waveguide

Authors: Ameya Galinde; Masoud Koochakzadeh; Abbas Abbaspour-Tamijani
Journal: Advances in Acoustics and Vibration
 (2012)

---

## Abstract

A semi-analytical method for finding the elastic modes propagating along the edge of an anisotropic semi-infinite plate is presented. Solutions are constructed as linear combinations of a finite number of the corresponding infinite plate modes with the constraint that they decay in the direction perpendicular to the edge and collectively satisfy the free boundary condition over the edge surface. Such modes that are confined to the edge can be used to approximate solutions of acoustic ridge waveguides whose supporting structures are sufficiently far away from the free edge. The semi-infinite plate or ridge is allowed to be oriented arbitrarily in the anisotropic crystal. Modifications to the theory to find symmetric and antisymmetric solutions for special crystal orientations are also presented. Accuracy of the solutions can be improved by including more plate modes in the series. Numerical techniques to find modal dispersion relations and orientation dependent modal behavior, are discussed. Results for ridges etched in single crystal Silicon are found to be in good agreement with Finite Element simulations. It is found that variations in modal phase velocity with respect to crystal orientation are not significant, suggesting that anisotropy may not be a critical issue while designing ridge waveguides in Silicon.


---

## Body

## 1. Introduction

The acoustic ridge waveguide has applications in many fields like microwave engineering and ultrasonics. Due to the low velocity of acoustic waves, compared to electromagnetic waves, miniature scale acoustic transmission lines can be fabricated at microwave frequencies for use in filters and resonators. An acoustic transmission line made from a vertical ridge etched in silicon is of special interest because of its ease of fabrication and small footprint. The modes of a ridge waveguide decay perpendicular to the edge surface and exhibit strong confinement to the free edge. This makes the ridge waveguide suitable for making low-loss acoustic transmission lines. The symmetric mode of the ridge also shows little dispersion over a wide frequency range and is similar to surface waves.A significant body of work can be found in the literature, addressing the modal analysis of acoustic ridge waveguides. However, most of the existing work focuses on isotropic ridge waveguides. The ridge as a topographic waveguide has been examined in [1]. Experimental results for a ridge made of the alloy duralumin are presented in [2]. A rectangular overlay waveguide over a rigid substrate has been analyzed in [3, 4]. A finite element method analysis of isotropic ridges has been presented in [5]. Using modes of an isotropic infinite plate, approximate modes of a semi-infinite thin ridge are obtained in [6]. Microwave circuit models for the solid-solid interface and free surface of a solid are given in [7, 8] and used in a transverse resonance formulation to derive approximate dispersion relations for an isotropic thin ridge [9, 10]. Flexural edge waves in isotropic semi-infinite plates have been studied in [11]. Existence of edge waves in orthotropic thin plates has been shown in [12]. Flexural edge waves on ridges in monoclinic crystals with the ridge lying in the plane of crystal mirror symmetry are examined in [13]. Their uniqueness and a method to obtain their velocity based on impedance calculations has also been presented. Using the modes of an isotropic infinite plate, 3-dimensional symmetric edge waves in semi-infinite plates have been analyzed in [14]. The only possible closed form expression obtained is for the case where the Poisson’s ratio is zero. Some recent work using finite element method for ridges on a rigid substrate has been reported in [15].An anisotropic ridge or a rectangular acoustic beam, in general, do not have exact analytical solutions. Some numerical methods must be used to arrive at an approximate solution. This paper aims to arrive at semianalytical solutions for the modes propagating along the edge of a semi-infinite plate. The basic idea is to express the modes of the semi-infinite plate in terms of the modes of an infinite plate. In the past, a similar semi-analytical approach has been used to find surface waves in infinite anisotropic half spaces and modes of anisotropic infinite plates (but not semi-infinite plates or ridge waveguides) [16–19], starting from an orthogonal and complete basis set of plane wave solutions in an infinite anisotropic crystal. Only, modes that decay in depth are considered in this paper. Normally a ridge waveguide is excited from the top, and so the modes of interest are only the ones that decay in depth. Thus, the solutions obtained for the semi-infinite plate can be taken as approximate solutions to the modes of an acoustic ridge waveguide with supporting structures that are sufficiently far away from the free edge such that the reflections from the interface can be neglected. The words ridge and semi-infinite plate are used interchangeably in this paper. A ridge waveguide having a width of wR is illustrated in Figure 1. The coordinate system is chosen such that the edge surface is in the y-z plane, the top surface normal is in x̂ direction, and propagation is along z-axis (lying on the middle of the top plane). This coordinate system will be used in all subsequent analysis. The ridge is made of a homogenous anisotropic material confined to the |y|≤wR/2, x≤0 region.(a) and (b) illustrate the geometry and coordinate system of the ridge waveguide. (c) shows the cross-section of a practical finite ridge waveguide etched in Silicon.
(a)(b)(c)Section2 outlines basic elasticity theory and a methodology for finding the modes of an anisotropic plate followed by the anisotropic ridge. Section 3 discusses the numerical tricks and issues involved in implementing the methods discussed in Section 2, along with a detailed example of a ridge etched in Silicon to illustrate the process step by step. In Section 4, results of the theory presented in this paper are verified with a finite element simulation in ANSYS, followed by dispersion curves and behavior of modes for ridges etched in Silicon at various angles. Finally, Section 5 summarizes and concludes the discussion. An Appendix is provided at the end, wherein some of the equations from Section 2 are expressed in greater detail.


## 2. Theory

A brief overview of elasticity theory and a mathematical description of the general approach are provided in this section. The methodology given here can, in general, be used to find propagating modes bound to the free edge of a ridge waveguide made of an anisotropic material with arbitrary crystal orientation. There are special cases where the plane of the ridge lies along one of the symmetry planes of the crystal. In such a situation, it is possible to find approximate symmetric and antisymmetric solutions. In this paper, all matrices are represented by bold symbols, unit vectors by capped symbols, vectors by symbols with an arrow on top, and matrix transpose with a superscriptT.


### 2.1. Elastic Wave Equation

The homogenous time harmonic acoustic field equations are [20](1)∇·T⃗=jρωv⃗,c∇sv⃗=jωT⃗,
where, ω is the frequency, ρ is the material density, v⃗=[vx,vy,vz]T=[v1,v2,v3]T is the particle velocity field, T⃗=[Txx,Tyy,Tzz,Tyz,Txz,Txy]T=[T1,T2,T3,T4,T5,T6]Tis the abbreviated notation for the second-order stress tensor:(2)T=[T1T6T5T6T2T4T5T4T3],
and c is the symmetric 6 × 6 elastic stiffness matrix in abbreviated notation after a bond transformation. The bond transformation expressed in (3) rotates the original stiffness matrix in crystal coordinate axes, cE, to the ridge coordinate axes using the bond transformation matrix, Mb, given in [20]:(3)c=MbcEMbT.
In a general triclinic crystal, which only has center point symmetry and no planes of mirror symmetry, cEhas 21 independent constants. For a monoclinic crystal, which has only one plane of mirror symmetry, the stiffness matrix simplifies to 13 independent constants. If the plane of mirror symmetry lies in the x-z plane, the corresponding monoclinic stiffness matrix can be used to obtain simplified expressions for symmetric and antisymmetric solutions to the plate or ridge. Expressions for ∇· and ∇s are given in (A.1) and (A.2), respectively, in the appendix.The field equations can be rewritten to give the following wave equation and stress field in terms of the velocity field:(4)[∇·(c∇s)+ρω2]v⃗=0,(5)T⃗=-j(c∇sv⃗)ω.
For uniform plane waves of the form v⃗upw=e-j(kzz+kxx+kyy)û, the elastic wave equation reduces to the matrix equation in (6), also known as the Christoffel  equation (for details see (A.3) to (A.6) in the appendix):(6)Γû=0.
Here, û=[ux,uy,uz]T is the velocity polarization, and k⃗=[kx,ky,kz]T is the propagation vector with k=|k⃗|. The solutions are obtained by setting the characteristic determinant equal to zero. For a given material and frequency, the surfaces in k⃗-space satisfying (7) correspond to the various slowness surfaces of the infinite crystal. For an isotropic medium, these surfaces are spheres, yielding the longitudinal and transverse solutions. In view of the ridge problem, kz is the propagation direction and has to be real, while kx must be complex to allow for decay in the transverse direction. By setting kz∈ℝ and kx∈ℂ, a 6th degree polynomial in ky is obtained, whose roots kym, in general, cannot be found analytically as stated by the Abel-Ruffini theorem [21] and must be found numerically as functions of kx and kz as kym=kym(kx,kz). If kym is not purely real, then it appears as a pair with its conjugate kym*. If the crystal has mirror symmetry in the x-z plane, and if kym is not purely imaginary, then it also appears as a pair with -kym. For every kym, the corresponding ûm is the eigenvector corresponding to the zero eigenvalue. In practice, where kym’s are found numerically, ûm may be found as the eigenvector corresponding to the smallest eigenvalue λΓ as described in (8):(7)det⁡(Γ(kx,kym,kz))=0,(8)Γ(kx,kym,kz)ûm=min⁡⁡(|λΓ|)ûm,m=1⋯6.


### 2.2. Free Anisotropic Plate

Since the wave equation is a linear equation, a linear combination of the plane wave solutions to the wave equation can be used to construct a solution for a free infinite plate, as shown in (9):(9)v⃗plate=∑m=16bmv⃗mupw.
If the plate is parallel to a plane of mirror symmetry of the crystal, then the roots of (7) appear in pairs of ±kym with polarizations [ux,±uy,uz]T. In such cases, the plate modes are found to be symmetric (even) or antisymmetric (odd) and can be expressed as shown in (10) and (11):(10)v⃗plate-s=e-jkzze-jkxx∑m=13bms[cos⁡⁡(kymy)uxm-jsin⁡⁡(kymy)uymcos⁡⁡(kymy)uzm],(11)v⃗plate-as=e-jkzze-jkxx∑m=13bmas[-jsin⁡⁡(kymy)uxmcos⁡⁡(kymy)uym-jsin⁡⁡(kymy)uzm].
The stress field can be calculated using (5) and made to satisfy the following free boundary condition on 3 traction components at each free surface with a surface normal n̂:(12)Tplaten̂|atfreeboundary=0.
The 2 free surfaces of the plate give rise to a matrix equation in 6 unknown coefficients represented by b̂=[b1,b2,b3,b4,b5,b6]T. Without loss of generality, the boundary conditions are applied at x=0, and z=0, as shown in (13) (see (A.7) in the appendix). In the symmetric or antisymmetric cases, the boundary condition needs to be applied only on one side of the plate with a coefficient vector b̂s/as=[b1,b2,b3]T to give a 3×3 matrix equation, as shown in (14) (see (A.8) and (A.9) in the appendix). Note that the symmetric cases have superscript s and antisymmetric cases have superscript as:(13)Bb̂=[Tplateŷ|x=0,y=wR/2,z=0-Tplateŷ|x=0,y=-wR/2,z=0]=0,(14)Bs/asb̂s/as=Tplate-s/asŷ|x=0,y=wR/2,z=0=0.
Exact solutions to these equations are obtained, where the characteristic determinant in (15) goes to zero. In general, it is not possible to find closed form analytical expressions for the solutions of this determinant. A practical approach is to perform a two-dimensional sweep in the kx-kz plane and search for the zeros of the determinant (note that kymn=kymn(kxn,kz) according to (7)). For the infinite plate, the search is conducted for real values of kx and kz, and results in closed contour in the kx-kz plane, or kxn=kxn(kz), corresponding to different propagating modes of the plate at the frequency of analysis. For the purpose of constructing the ridge modes, however, one needs to sweep kx in the upper half of the complex plane, Im(kx≥0), which corresponds to the edge-bound solutions. If kxn is not purely real, then it appears as a pair with its conjugate kxn*. If the crystal has mirror symmetry in the y-z plane, and if kxn is not purely imaginary, then it also appears as a pair with -kxn. For every kxn, the corresponding coefficient vector b̂n, theoretically representing the nullspace, is approximated as the eigenvector corresponding to the smallest eigenvalue λB as described in (16):(15)det⁡(B(kx,kz))=0,(16)B(kxn)b̂n=min⁡⁡(|λB|)b̂n,∀n.
In Section 3, issues involved with locating and tracking the zeros of det(B(kx,kz)) or the poles of log(|det(B(kx,kz))|) are discussed.


### 2.3. Anisotropic Semi-Infinite Plate or Ridge

General solutions to the ridge can be constructed as a linear combination ofN modes of the plate as shown in (17). Also, the symmetric and antisymmetric ridge modes can be constructed from the symmetric and antisymmetric modes of the plate. As mentioned above, only plate modes with Im(kx≥0) contribute to edge-bound solutions. If the crystal has mirror symmetry with respect to the y-z plane, and if kxn is not purely imaginary, then the modes should be chosen in pairs of [kxn,-kxn*] to enable standing waves that decay in depth. Accuracy of the ridge solution can be increased by including more plate modes in the series, that is, by increasing N. A simple convergence study to illustrate this has been presented in Section 3.3:(17)v⃗ridge-(s/as)=∑n=1Nanv⃗nplate-(s/as).
Without loss of generality, the boundary condition is applied at z=0 on the top surface of the ridge to yield the matrix equation (18), which should be satisfied over the range |y|≤wR/2. Since the boundary conditions are a function of y, in general it is not possible to obtain an exact solution to the ridge problem. There are two main steps in finding the ridge solution: (1) constructing an appropriate matrix A corresponding to the field quantities that must vanish on the top surface of the ridge and (2) finding the value of kz that makes this matrix singular and the corresponding nullspace â. One way of expressing the boundary conditions is by expanding them in terms of appropriate domain functions over the range of y. Point matching is a special case, wherein the domain functions are delta functions. Another approach is to expand the boundary conditions as a Taylor series about a fixed point in y0. For the case of free boundary conditions, setting the traction component and its higher order derivatives at a fixed point in y equal to zero, the Taylor series expansion about that point approaches zero, thereby forcing the traction in a neighborhood of that point to zero. Here, we have chosen the latter method to construct the approximate boundary conditions as shown in (19) (also see (A.12) to (A.16) in the appendix), where â=[a1,…,aN]T is the coefficient vector weighing each plate mode:(18)Tridgex̂|x=0,|y|≤wR/2,z=0=0,(19)dpdypTridgex̂∣x=0,y=y0,z=0=[A⃗1(p,y0)|⋯|A⃗N(p,y0)]â=0;p=0,1,…,P.
For 0<|y0|<wR/2, (19) gives 3 nontrivial equations for each value of p. For y=wR/2 or -wR/2, (19) only gives 2 nontrivial equations for p=0 because one traction component Txy is already zero by the plate boundary conditions. For y0=0, and the symmetric case, Txy, and it even, derivatives are zero at y0, and (19) produces only 2 non-trivial equations for the even values of p. Similarly, the odd derivatives of Txx and Txz are zero, yielding only 1 non-trivial equation for every odd value of p. For y0=0 and antisymmetric case, (19) gives 1 and 2 non-trivial equations for each even or odd value of p, respectively. The number of equations obtained for the pth derivative for different choices of y0 in the general, symmetric, and antisymmetric cases are summarized in Table 1. The non-trivial boundary condition equations are combined into a single matrix equation (20) as follows:(20)Aâ=0.
In order to have a square boundary condition matrix A(s/as), the values of p much be appropriately chosen such that the total number of equations is N.Table 1
Number of ridge boundary conditions obtained at each point on they-axis for the pth derivative for the general, symmetric, and the antisymmetric cases.
Number of equationsy0=00<|y0|<wR/2y0=±wR/2GeneralP=0332Odd p333Even p333SymmetricP=0232Odd p133Even p233AntisymmetricP=0132Odd p233Even p133After constructing a suitableA matrix encapsulating the boundary conditions of the ridge, the next problem is to find its nullspace â. If the matrix is not square, then the nullspace can be found either by using the method of Lagrange multipliers or by performing a search in the complex space ℂN of â. If the matrix is square, then the characteristic determinant in (21) is set to zero, and the values kzl corresponding to the propagating modes of the ridge are found numerically. The edge bound modes only correspond to those values of kzl for which none of the contributing kxnl’s with nonzero coefficients are purely real. For every kzl, the respective nullspace âl representing the coefficient vector is approximated as the eigenvector corresponding to the smallest eigenvalue λA as described in (22). As functions of frequency, the dispersion relations of the ridge modes can be expressed as  kzl=kzl(ω), kxnl=kxnl(kzl,ω) and, kymnl=kymnl(kxnl,kzl,ω):(21)det⁡(A(kz))=0,(22)A(kzl)âl=min⁡⁡(|λA|)âl,∀l.
For thicker ridges, in order to obtain higher order ridge modes, many plate modes must be used in the expansion. Since the solutions obtained are only approximates, the only way to ensure that the solutions are authentic is to check for convergence of kzl. If kzlN is the propagation constant obtained using N plate modes, then limN→∞kzlN=kzl.


## 2.1. Elastic Wave Equation

The homogenous time harmonic acoustic field equations are [20](1)∇·T⃗=jρωv⃗,c∇sv⃗=jωT⃗,
where, ω is the frequency, ρ is the material density, v⃗=[vx,vy,vz]T=[v1,v2,v3]T is the particle velocity field, T⃗=[Txx,Tyy,Tzz,Tyz,Txz,Txy]T=[T1,T2,T3,T4,T5,T6]Tis the abbreviated notation for the second-order stress tensor:(2)T=[T1T6T5T6T2T4T5T4T3],
and c is the symmetric 6 × 6 elastic stiffness matrix in abbreviated notation after a bond transformation. The bond transformation expressed in (3) rotates the original stiffness matrix in crystal coordinate axes, cE, to the ridge coordinate axes using the bond transformation matrix, Mb, given in [20]:(3)c=MbcEMbT.
In a general triclinic crystal, which only has center point symmetry and no planes of mirror symmetry, cEhas 21 independent constants. For a monoclinic crystal, which has only one plane of mirror symmetry, the stiffness matrix simplifies to 13 independent constants. If the plane of mirror symmetry lies in the x-z plane, the corresponding monoclinic stiffness matrix can be used to obtain simplified expressions for symmetric and antisymmetric solutions to the plate or ridge. Expressions for ∇· and ∇s are given in (A.1) and (A.2), respectively, in the appendix.The field equations can be rewritten to give the following wave equation and stress field in terms of the velocity field:(4)[∇·(c∇s)+ρω2]v⃗=0,(5)T⃗=-j(c∇sv⃗)ω.
For uniform plane waves of the form v⃗upw=e-j(kzz+kxx+kyy)û, the elastic wave equation reduces to the matrix equation in (6), also known as the Christoffel  equation (for details see (A.3) to (A.6) in the appendix):(6)Γû=0.
Here, û=[ux,uy,uz]T is the velocity polarization, and k⃗=[kx,ky,kz]T is the propagation vector with k=|k⃗|. The solutions are obtained by setting the characteristic determinant equal to zero. For a given material and frequency, the surfaces in k⃗-space satisfying (7) correspond to the various slowness surfaces of the infinite crystal. For an isotropic medium, these surfaces are spheres, yielding the longitudinal and transverse solutions. In view of the ridge problem, kz is the propagation direction and has to be real, while kx must be complex to allow for decay in the transverse direction. By setting kz∈ℝ and kx∈ℂ, a 6th degree polynomial in ky is obtained, whose roots kym, in general, cannot be found analytically as stated by the Abel-Ruffini theorem [21] and must be found numerically as functions of kx and kz as kym=kym(kx,kz). If kym is not purely real, then it appears as a pair with its conjugate kym*. If the crystal has mirror symmetry in the x-z plane, and if kym is not purely imaginary, then it also appears as a pair with -kym. For every kym, the corresponding ûm is the eigenvector corresponding to the zero eigenvalue. In practice, where kym’s are found numerically, ûm may be found as the eigenvector corresponding to the smallest eigenvalue λΓ as described in (8):(7)det⁡(Γ(kx,kym,kz))=0,(8)Γ(kx,kym,kz)ûm=min⁡⁡(|λΓ|)ûm,m=1⋯6.


## 2.2. Free Anisotropic Plate

Since the wave equation is a linear equation, a linear combination of the plane wave solutions to the wave equation can be used to construct a solution for a free infinite plate, as shown in (9):(9)v⃗plate=∑m=16bmv⃗mupw.
If the plate is parallel to a plane of mirror symmetry of the crystal, then the roots of (7) appear in pairs of ±kym with polarizations [ux,±uy,uz]T. In such cases, the plate modes are found to be symmetric (even) or antisymmetric (odd) and can be expressed as shown in (10) and (11):(10)v⃗plate-s=e-jkzze-jkxx∑m=13bms[cos⁡⁡(kymy)uxm-jsin⁡⁡(kymy)uymcos⁡⁡(kymy)uzm],(11)v⃗plate-as=e-jkzze-jkxx∑m=13bmas[-jsin⁡⁡(kymy)uxmcos⁡⁡(kymy)uym-jsin⁡⁡(kymy)uzm].
The stress field can be calculated using (5) and made to satisfy the following free boundary condition on 3 traction components at each free surface with a surface normal n̂:(12)Tplaten̂|atfreeboundary=0.
The 2 free surfaces of the plate give rise to a matrix equation in 6 unknown coefficients represented by b̂=[b1,b2,b3,b4,b5,b6]T. Without loss of generality, the boundary conditions are applied at x=0, and z=0, as shown in (13) (see (A.7) in the appendix). In the symmetric or antisymmetric cases, the boundary condition needs to be applied only on one side of the plate with a coefficient vector b̂s/as=[b1,b2,b3]T to give a 3×3 matrix equation, as shown in (14) (see (A.8) and (A.9) in the appendix). Note that the symmetric cases have superscript s and antisymmetric cases have superscript as:(13)Bb̂=[Tplateŷ|x=0,y=wR/2,z=0-Tplateŷ|x=0,y=-wR/2,z=0]=0,(14)Bs/asb̂s/as=Tplate-s/asŷ|x=0,y=wR/2,z=0=0.
Exact solutions to these equations are obtained, where the characteristic determinant in (15) goes to zero. In general, it is not possible to find closed form analytical expressions for the solutions of this determinant. A practical approach is to perform a two-dimensional sweep in the kx-kz plane and search for the zeros of the determinant (note that kymn=kymn(kxn,kz) according to (7)). For the infinite plate, the search is conducted for real values of kx and kz, and results in closed contour in the kx-kz plane, or kxn=kxn(kz), corresponding to different propagating modes of the plate at the frequency of analysis. For the purpose of constructing the ridge modes, however, one needs to sweep kx in the upper half of the complex plane, Im(kx≥0), which corresponds to the edge-bound solutions. If kxn is not purely real, then it appears as a pair with its conjugate kxn*. If the crystal has mirror symmetry in the y-z plane, and if kxn is not purely imaginary, then it also appears as a pair with -kxn. For every kxn, the corresponding coefficient vector b̂n, theoretically representing the nullspace, is approximated as the eigenvector corresponding to the smallest eigenvalue λB as described in (16):(15)det⁡(B(kx,kz))=0,(16)B(kxn)b̂n=min⁡⁡(|λB|)b̂n,∀n.
In Section 3, issues involved with locating and tracking the zeros of det(B(kx,kz)) or the poles of log(|det(B(kx,kz))|) are discussed.


## 2.3. Anisotropic Semi-Infinite Plate or Ridge

General solutions to the ridge can be constructed as a linear combination ofN modes of the plate as shown in (17). Also, the symmetric and antisymmetric ridge modes can be constructed from the symmetric and antisymmetric modes of the plate. As mentioned above, only plate modes with Im(kx≥0) contribute to edge-bound solutions. If the crystal has mirror symmetry with respect to the y-z plane, and if kxn is not purely imaginary, then the modes should be chosen in pairs of [kxn,-kxn*] to enable standing waves that decay in depth. Accuracy of the ridge solution can be increased by including more plate modes in the series, that is, by increasing N. A simple convergence study to illustrate this has been presented in Section 3.3:(17)v⃗ridge-(s/as)=∑n=1Nanv⃗nplate-(s/as).
Without loss of generality, the boundary condition is applied at z=0 on the top surface of the ridge to yield the matrix equation (18), which should be satisfied over the range |y|≤wR/2. Since the boundary conditions are a function of y, in general it is not possible to obtain an exact solution to the ridge problem. There are two main steps in finding the ridge solution: (1) constructing an appropriate matrix A corresponding to the field quantities that must vanish on the top surface of the ridge and (2) finding the value of kz that makes this matrix singular and the corresponding nullspace â. One way of expressing the boundary conditions is by expanding them in terms of appropriate domain functions over the range of y. Point matching is a special case, wherein the domain functions are delta functions. Another approach is to expand the boundary conditions as a Taylor series about a fixed point in y0. For the case of free boundary conditions, setting the traction component and its higher order derivatives at a fixed point in y equal to zero, the Taylor series expansion about that point approaches zero, thereby forcing the traction in a neighborhood of that point to zero. Here, we have chosen the latter method to construct the approximate boundary conditions as shown in (19) (also see (A.12) to (A.16) in the appendix), where â=[a1,…,aN]T is the coefficient vector weighing each plate mode:(18)Tridgex̂|x=0,|y|≤wR/2,z=0=0,(19)dpdypTridgex̂∣x=0,y=y0,z=0=[A⃗1(p,y0)|⋯|A⃗N(p,y0)]â=0;p=0,1,…,P.
For 0<|y0|<wR/2, (19) gives 3 nontrivial equations for each value of p. For y=wR/2 or -wR/2, (19) only gives 2 nontrivial equations for p=0 because one traction component Txy is already zero by the plate boundary conditions. For y0=0, and the symmetric case, Txy, and it even, derivatives are zero at y0, and (19) produces only 2 non-trivial equations for the even values of p. Similarly, the odd derivatives of Txx and Txz are zero, yielding only 1 non-trivial equation for every odd value of p. For y0=0 and antisymmetric case, (19) gives 1 and 2 non-trivial equations for each even or odd value of p, respectively. The number of equations obtained for the pth derivative for different choices of y0 in the general, symmetric, and antisymmetric cases are summarized in Table 1. The non-trivial boundary condition equations are combined into a single matrix equation (20) as follows:(20)Aâ=0.
In order to have a square boundary condition matrix A(s/as), the values of p much be appropriately chosen such that the total number of equations is N.Table 1
Number of ridge boundary conditions obtained at each point on they-axis for the pth derivative for the general, symmetric, and the antisymmetric cases.
Number of equationsy0=00<|y0|<wR/2y0=±wR/2GeneralP=0332Odd p333Even p333SymmetricP=0232Odd p133Even p233AntisymmetricP=0132Odd p233Even p133After constructing a suitableA matrix encapsulating the boundary conditions of the ridge, the next problem is to find its nullspace â. If the matrix is not square, then the nullspace can be found either by using the method of Lagrange multipliers or by performing a search in the complex space ℂN of â. If the matrix is square, then the characteristic determinant in (21) is set to zero, and the values kzl corresponding to the propagating modes of the ridge are found numerically. The edge bound modes only correspond to those values of kzl for which none of the contributing kxnl’s with nonzero coefficients are purely real. For every kzl, the respective nullspace âl representing the coefficient vector is approximated as the eigenvector corresponding to the smallest eigenvalue λA as described in (22). As functions of frequency, the dispersion relations of the ridge modes can be expressed as  kzl=kzl(ω), kxnl=kxnl(kzl,ω) and, kymnl=kymnl(kxnl,kzl,ω):(21)det⁡(A(kz))=0,(22)A(kzl)âl=min⁡⁡(|λA|)âl,∀l.
For thicker ridges, in order to obtain higher order ridge modes, many plate modes must be used in the expansion. Since the solutions obtained are only approximates, the only way to ensure that the solutions are authentic is to check for convergence of kzl. If kzlN is the propagation constant obtained using N plate modes, then limN→∞kzlN=kzl.


## 3. Numerical Calculations

### 3.1. Pole Tracking

As mentioned earlier, the solutions of (15), kxn(kz) cannot be found analytically; hence, the zeros of the determinant of B or the poles of log(|det(B)|) must be found numerically. Performing a raster scan in the complex kx plane for each kzto track the poles is a computationally intensive task. To expedite the process, some heuristic methods can be employed to predict and locate the poles. A high-resolution raster scan is performed for a few initial values of kz to identify the precise location of poles. The subset of the plate modes used for synthesizing the ridge solution must be carefully selected, and these are then tracked over successive values of kz to ensure the continuity of the determinant. Using a polynomial fit, the approximate location of poles or kxn’s for each mode are estimated for successive kz’s. A more precise location upto a desired degree of accuracy is then arrived at by using a combination of gradient descent and a 2-dimensional version of the direct line search algorithms. A good prediction helps to make the gradient descent approach converge faster and reduces the search area for the direct search algorithm. For the gradient descent algorithm, the estimates are advanced at each step such that the drop in the value of the function remains constant. Usually the gradient descent algorithm is used to track minimax and tends to become slower when it approaches a zero. In the present formulation, since the algorithm is operating in the vicinity of a pole, it tends to converge faster as the estimate advances towards the location of the pole. In some cases, when kz is scanned, it is observed that two poles come close and then recede away in opposite orthogonal directions, rendering the predicted values grossly off the mark. Pathological cases such as these can cause the starting point to lie on a saddle point or a flat area. These are taken care of by limiting the jump in the estimate to within a certain radius from the current position or by resorting to the direct search algorithm. The use of prediction and precision search algorithms makes the overall pole-tracking process hundreds of times faster than a brute force raster scan. This also helps in tracking the behavior of individual plate modes and gives an estimate for their computed dispersion relations.Once the required plate modes are tracked overkz to get kxn(kz), the appropriate boundary conditions (19) for the ridge top surface are applied on them to construct the A matrix. The zeros of the determinant of A or the poles of log(|det(A)|) represent the propagation constants of the ridge modes. Locations of these poles are kzl with the constituent plate modes, having kx=kxnl(kzl). This entire process of pole-tracking and application of ridge boundary condition yields the ridge modes only at a given frequency ω or ridge thickness, and at a given ridge orientation in the anisotropic crystal. In order to get the full dispersion relations or behavior of the ridge modes over frequency or orientation, this entire process has to be repeated for each frequency point or ridge orientation. This process can again be time-consuming and may need significant manual intervention. To automate this process and make it much faster, the plate modes are directly tracked over ridge modes. That is, instead of tracking kxn(kz), kzl(ω), and kxnl(kz,ω), or kzl(c) and kxnl(kz,c), are tracked using a combination of gradient descent and direct line search.


### 3.2. Spurious Poles

Whenlog(|det(B)|) is plotted, some spurious poles are observed that do not correspond to any plate modes, that is, they do not satisfy the plate boundary conditions. These poles appear at locations where kym’s become redundant with the same velocity polarization vectors. Such degeneracies occur when two kym’s meet and split in orthogonal directions. In order to remove these spurious peaks and ensure smooth functioning of the automatic pole tracking algorithm, the 6 × 6 matrix B̃ in (23), is defined for the general case with m=1⋯6. Correspondingly, a 3×3 matrix B̃(s/as) in (24) and (25), is defined for the symmetric and antisymmetric cases, respectively, with m=1…3. The determinant of the main plate boundary condition matrix B or Bs or Bas is divided by the determinant of this pole nullifying matrix B̃ or B̃s or B̃as, respectively, to cancel out the spurious poles in the logarithm of their absolute values, as shown in (26). The pole tracking algorithms use d(kx,kz) while tracking the plate modes:(23)B̃:[B̃1m,B̃2m,B̃3m,B̃4m,B̃5m,B̃6m]T≜[uxm,uym,uzm,kymuxm,kymuym,kymuzm]T,(24)B̃s:[B̃1ms,B̃2ms,B̃3ms]T≜[uxm,kymuym,uzm]T,(25)B̃as:[B̃1mas,B̃2mas,B̃3mas]T≜[kymuxm,uym,kymuzm]T,(26)⁡d(kx,kz)=log⁡(|det⁡(B(s/as)(kx,kz))det⁡(B̃(s/as)(kx,kz))|).


### 3.3. Example Case

An example is presented to illustrate the procedure for finding ridge modes starting from plate modes. All the numerical computations have been done in MATLAB. Symbolic math has been used to compute the long expressions resulting for the general anisotropic cases. Propagating modes at 350 MHz are calculated for a 2-micron wide ridge deeply etched in a〈100〉 Silicon wafer with the z-axis along the 〈001〉 crystal axis. Silicon is a cubic crystal with a density ρ=2300kg/m3 and 3 independent stiffness constants: c11=165.7GPa, c12=63.9GPa, and c44=79.56GPa.The first step is to identify the plate modes for a 2 micron wide plate in Silicon with the same orientation by doing a raster scan in the complexkx plane for kz=0, as shown in Figure 2. Plate boundary conditions for the general case have been used in this calculation; hence, m=1⋯6. The axes in this figure and subsequent k-space figures are normalized to the propagation constant of shear waves in bulk Silicon defined as ksh=ωρ/c44. The plate modes are then tracked over kz to obtain the mode constellation shown in Figure 3. Individual plate modes have been labeled in this figure. The apparent convoluted evolutions of the curves are manifestations of the inherent modal dispersion relations of the plate modes in the complex plane. The n=1, n=2, and n=3 plate modes are purely propagating (real kx) for small values of kz/ksh and can be identified as the fundamental longitudinal L1, fundamental shear SH0, and fundamental flexural F1 modes of a plate, respectively. Since the crystal orientation is such that it has mirror symmetry in the y-z plane, the plate modes are seen to evolve symmetrically about the imaginary kx-axis. At higher values of kz, the n=1 and n=2 combine to give a hybrid shear-longitudinal plate mode. This is in contrast to the isotropic plate, where the SH0, and L1 modes can be separately distinguished for all values of kz.Figure 2
Plate modes of a 2-micron wide plate withx-axis along the 〈001〉 crystal axis and z-axis along the 〈001〉 crystal axis at 350 MHz for kz=0.Figure 3
Plate modes of Figure2 tracked over kz at 350 MHz. The n=1, n=2, and n=3 plate modes are identified as the L1, SH0, and F1 modes of a plate, respectively.In order to find the modes of the ridge, only 3 symmetric and 3 antisymmetric plate modes withkxn (n=1⋯6) are used in the expansion, and the top boundary conditions are applied at y0=0 and p=0,1 (see Table 1), separately, for the symmetric and antisymmetric cases. This results in 3 × 3 square As and Aas matrices whose determinants are plotted in Figure 4(a) for the symmetric case and Figure 4(b) for the antisymmetric case. Two zeros are observed in this plot, marked as l=1 and l=2. Each zero corresponds to a propagating ridge mode. Values of kzl and the kxnl’s of the 6 constituent plate modes along with the magnitude of their coefficients anl are listed in Table 2. It can be seen that the main contribution to l=1 is from the 2 fundamental symmetric plate modes =1, and  n=2, both of which have complex kxn’s. Hence, this mode decays in depth and behaves as the fundamental propagating symmetric edge mode bound to the top surface of the ridge. The displacement field of this ridge mode, shown in Figure 5(a), resembles a surface wave. In fact, the propagation constant kz1 is almost equal to the propagation constant of a z-directed surface wave. For the case of l=2, the main contribution is from the fundamental antisymmetric plate mode n=3, which has a small but imaginary kx33. Hence, this can be called as the fundamental propagating antisymmetric ridge mode. Its displacement field, as shown in Figure 5(b), resembles a flexural plate mode decaying in depth. Comparing the values of kx11 and kx21 with kx33, it is seen that l=2 decays much slower compared to l=1. For a ridge waveguide standing on a substrate, in order to minimize reflections from the substrate, a higher aspect ratio (hR/wR in Figure 1(c)) should be used for the ridge carrying the antisymmetric mode. If it is carrying the symmetric mode, then the aspect ratio do not need to be as large. Particle velocity and stress fields within the cross-section for l=1 and l=2 are shown in Figure 6. In both cases, T2, T4, and T6 should be zero on the sides walls in order to satisfy boundary conditions of the plate, and T1, T5, and T6 should be zero on the top surface in order to satisfy boundary conditions of the ridge. In the present example, the values obtained for kymnl are found to be quite small compared to the width of the ridge. Thus, even though the boundary conditions were applied only at the center of the ridge, they appear to be satisfied across the entire range of y on the top surface to a good extent.Table 2
Constituent plate mode coefficients and propagation constants for the two cases identified in Figure4. l=1 is classified as the fundamental symmetric ridge mode and l=2 is classified as the fundamental antisymmetric ridge mode.
Mode coefficientsl=1(kzl/ksh=1.204)l=2(kzl/ksh=2.089)kxnl/ksh|anl|n=1Symmetric (L1)-0.5611+j0.56160.7071n=2Symmetric (SH0)0.5611+j0.56160.7071n=3Antisymmetric (F1)1.59970n=4Antisymmetric0+j2.35390n=5Antisymmetric0+j3.87510n=6Symmetric0+j8.29360.0017Determinant of the matrix of ridge boundary conditions separately applied to (a) 3 symmetric and (b) 3 antisymmetric plate modes tracked in Figure3. Two poles are identified with l=1 and l=2.
(a)(b)Displacement fields for the modes (a)l=1 and (b) l=2, identified in Figure 4. l=1 resembles a surface wave and l=2 resembles a flexural plate mode slowly decaying in depth.
(a)(b)Magnitudes of particle velocity and stress fields in the cross-section for the modes (a)l=1 and (b) l=2, identified in Figure 4. The x- and y-dimensions are in microns. In both cases, T2, T4, and T6 are zero on the sides walls, thus satisfying boundary conditions of the plate, and T1, T5, and T6 are zero on the top surface, thus satisfying boundary conditions of the ridge.
(a)
Magnitudes of particle velocity (m/s) and stress fields (Pa) in the cross-section for model=1(b)
Magnitudes of particle velocity (m/s) and stress fields (Pa) in the cross-section for model=2The ridge modes found here can also be obtained using the general formulation as a single expansion constituting all 6 plate modes with boundary conditions applied aty0=0 and p=0,1 (see Table 1). This results in a 6x6 square A matrix, whose determinant has zeros at the same values of kz as the ones found using the symmetric and antisymmetric formulations. In addition to these ridge solutions, some spurious zeros are observed that correspond to values of kz where the plate modes become degenerate, which happens when the curves of Figure 3 intersect. Such false zeros not representing a solution to the ridge can be rejected by verifying the boundary conditions from field plots similar to Figure 6.Earlier it was mentioned that, by increasing the number of plate modes included in the series, the accuracy of the ridge solution can be improved. This point has been illustrated for the current example through a simple convergence study, as presented in Figure7. It can be seen that both the symmetric and the antisymmetric solutions of the ridge indeed converge as the number of plate modes, N, is increased. The boundary conditions in each case are applied at various points in the range 0≤|y0|<wR/2 and p=0,1, depending on the number of plate modes chosen.A convergence study to illustrate that both the symmetric and the antisymmetric solutions of the ridge converge as the number of plate modes included in the series are increased.
(a)
Convergence study of the symmetric ridge mode(b)
Convergence study of the antisymmetric ridge mode


## 3.1. Pole Tracking

As mentioned earlier, the solutions of (15), kxn(kz) cannot be found analytically; hence, the zeros of the determinant of B or the poles of log(|det(B)|) must be found numerically. Performing a raster scan in the complex kx plane for each kzto track the poles is a computationally intensive task. To expedite the process, some heuristic methods can be employed to predict and locate the poles. A high-resolution raster scan is performed for a few initial values of kz to identify the precise location of poles. The subset of the plate modes used for synthesizing the ridge solution must be carefully selected, and these are then tracked over successive values of kz to ensure the continuity of the determinant. Using a polynomial fit, the approximate location of poles or kxn’s for each mode are estimated for successive kz’s. A more precise location upto a desired degree of accuracy is then arrived at by using a combination of gradient descent and a 2-dimensional version of the direct line search algorithms. A good prediction helps to make the gradient descent approach converge faster and reduces the search area for the direct search algorithm. For the gradient descent algorithm, the estimates are advanced at each step such that the drop in the value of the function remains constant. Usually the gradient descent algorithm is used to track minimax and tends to become slower when it approaches a zero. In the present formulation, since the algorithm is operating in the vicinity of a pole, it tends to converge faster as the estimate advances towards the location of the pole. In some cases, when kz is scanned, it is observed that two poles come close and then recede away in opposite orthogonal directions, rendering the predicted values grossly off the mark. Pathological cases such as these can cause the starting point to lie on a saddle point or a flat area. These are taken care of by limiting the jump in the estimate to within a certain radius from the current position or by resorting to the direct search algorithm. The use of prediction and precision search algorithms makes the overall pole-tracking process hundreds of times faster than a brute force raster scan. This also helps in tracking the behavior of individual plate modes and gives an estimate for their computed dispersion relations.Once the required plate modes are tracked overkz to get kxn(kz), the appropriate boundary conditions (19) for the ridge top surface are applied on them to construct the A matrix. The zeros of the determinant of A or the poles of log(|det(A)|) represent the propagation constants of the ridge modes. Locations of these poles are kzl with the constituent plate modes, having kx=kxnl(kzl). This entire process of pole-tracking and application of ridge boundary condition yields the ridge modes only at a given frequency ω or ridge thickness, and at a given ridge orientation in the anisotropic crystal. In order to get the full dispersion relations or behavior of the ridge modes over frequency or orientation, this entire process has to be repeated for each frequency point or ridge orientation. This process can again be time-consuming and may need significant manual intervention. To automate this process and make it much faster, the plate modes are directly tracked over ridge modes. That is, instead of tracking kxn(kz), kzl(ω), and kxnl(kz,ω), or kzl(c) and kxnl(kz,c), are tracked using a combination of gradient descent and direct line search.


## 3.2. Spurious Poles

Whenlog(|det(B)|) is plotted, some spurious poles are observed that do not correspond to any plate modes, that is, they do not satisfy the plate boundary conditions. These poles appear at locations where kym’s become redundant with the same velocity polarization vectors. Such degeneracies occur when two kym’s meet and split in orthogonal directions. In order to remove these spurious peaks and ensure smooth functioning of the automatic pole tracking algorithm, the 6 × 6 matrix B̃ in (23), is defined for the general case with m=1⋯6. Correspondingly, a 3×3 matrix B̃(s/as) in (24) and (25), is defined for the symmetric and antisymmetric cases, respectively, with m=1…3. The determinant of the main plate boundary condition matrix B or Bs or Bas is divided by the determinant of this pole nullifying matrix B̃ or B̃s or B̃as, respectively, to cancel out the spurious poles in the logarithm of their absolute values, as shown in (26). The pole tracking algorithms use d(kx,kz) while tracking the plate modes:(23)B̃:[B̃1m,B̃2m,B̃3m,B̃4m,B̃5m,B̃6m]T≜[uxm,uym,uzm,kymuxm,kymuym,kymuzm]T,(24)B̃s:[B̃1ms,B̃2ms,B̃3ms]T≜[uxm,kymuym,uzm]T,(25)B̃as:[B̃1mas,B̃2mas,B̃3mas]T≜[kymuxm,uym,kymuzm]T,(26)⁡d(kx,kz)=log⁡(|det⁡(B(s/as)(kx,kz))det⁡(B̃(s/as)(kx,kz))|).


## 3.3. Example Case

An example is presented to illustrate the procedure for finding ridge modes starting from plate modes. All the numerical computations have been done in MATLAB. Symbolic math has been used to compute the long expressions resulting for the general anisotropic cases. Propagating modes at 350 MHz are calculated for a 2-micron wide ridge deeply etched in a〈100〉 Silicon wafer with the z-axis along the 〈001〉 crystal axis. Silicon is a cubic crystal with a density ρ=2300kg/m3 and 3 independent stiffness constants: c11=165.7GPa, c12=63.9GPa, and c44=79.56GPa.The first step is to identify the plate modes for a 2 micron wide plate in Silicon with the same orientation by doing a raster scan in the complexkx plane for kz=0, as shown in Figure 2. Plate boundary conditions for the general case have been used in this calculation; hence, m=1⋯6. The axes in this figure and subsequent k-space figures are normalized to the propagation constant of shear waves in bulk Silicon defined as ksh=ωρ/c44. The plate modes are then tracked over kz to obtain the mode constellation shown in Figure 3. Individual plate modes have been labeled in this figure. The apparent convoluted evolutions of the curves are manifestations of the inherent modal dispersion relations of the plate modes in the complex plane. The n=1, n=2, and n=3 plate modes are purely propagating (real kx) for small values of kz/ksh and can be identified as the fundamental longitudinal L1, fundamental shear SH0, and fundamental flexural F1 modes of a plate, respectively. Since the crystal orientation is such that it has mirror symmetry in the y-z plane, the plate modes are seen to evolve symmetrically about the imaginary kx-axis. At higher values of kz, the n=1 and n=2 combine to give a hybrid shear-longitudinal plate mode. This is in contrast to the isotropic plate, where the SH0, and L1 modes can be separately distinguished for all values of kz.Figure 2
Plate modes of a 2-micron wide plate withx-axis along the 〈001〉 crystal axis and z-axis along the 〈001〉 crystal axis at 350 MHz for kz=0.Figure 3
Plate modes of Figure2 tracked over kz at 350 MHz. The n=1, n=2, and n=3 plate modes are identified as the L1, SH0, and F1 modes of a plate, respectively.In order to find the modes of the ridge, only 3 symmetric and 3 antisymmetric plate modes withkxn (n=1⋯6) are used in the expansion, and the top boundary conditions are applied at y0=0 and p=0,1 (see Table 1), separately, for the symmetric and antisymmetric cases. This results in 3 × 3 square As and Aas matrices whose determinants are plotted in Figure 4(a) for the symmetric case and Figure 4(b) for the antisymmetric case. Two zeros are observed in this plot, marked as l=1 and l=2. Each zero corresponds to a propagating ridge mode. Values of kzl and the kxnl’s of the 6 constituent plate modes along with the magnitude of their coefficients anl are listed in Table 2. It can be seen that the main contribution to l=1 is from the 2 fundamental symmetric plate modes =1, and  n=2, both of which have complex kxn’s. Hence, this mode decays in depth and behaves as the fundamental propagating symmetric edge mode bound to the top surface of the ridge. The displacement field of this ridge mode, shown in Figure 5(a), resembles a surface wave. In fact, the propagation constant kz1 is almost equal to the propagation constant of a z-directed surface wave. For the case of l=2, the main contribution is from the fundamental antisymmetric plate mode n=3, which has a small but imaginary kx33. Hence, this can be called as the fundamental propagating antisymmetric ridge mode. Its displacement field, as shown in Figure 5(b), resembles a flexural plate mode decaying in depth. Comparing the values of kx11 and kx21 with kx33, it is seen that l=2 decays much slower compared to l=1. For a ridge waveguide standing on a substrate, in order to minimize reflections from the substrate, a higher aspect ratio (hR/wR in Figure 1(c)) should be used for the ridge carrying the antisymmetric mode. If it is carrying the symmetric mode, then the aspect ratio do not need to be as large. Particle velocity and stress fields within the cross-section for l=1 and l=2 are shown in Figure 6. In both cases, T2, T4, and T6 should be zero on the sides walls in order to satisfy boundary conditions of the plate, and T1, T5, and T6 should be zero on the top surface in order to satisfy boundary conditions of the ridge. In the present example, the values obtained for kymnl are found to be quite small compared to the width of the ridge. Thus, even though the boundary conditions were applied only at the center of the ridge, they appear to be satisfied across the entire range of y on the top surface to a good extent.Table 2
Constituent plate mode coefficients and propagation constants for the two cases identified in Figure4. l=1 is classified as the fundamental symmetric ridge mode and l=2 is classified as the fundamental antisymmetric ridge mode.
Mode coefficientsl=1(kzl/ksh=1.204)l=2(kzl/ksh=2.089)kxnl/ksh|anl|n=1Symmetric (L1)-0.5611+j0.56160.7071n=2Symmetric (SH0)0.5611+j0.56160.7071n=3Antisymmetric (F1)1.59970n=4Antisymmetric0+j2.35390n=5Antisymmetric0+j3.87510n=6Symmetric0+j8.29360.0017Determinant of the matrix of ridge boundary conditions separately applied to (a) 3 symmetric and (b) 3 antisymmetric plate modes tracked in Figure3. Two poles are identified with l=1 and l=2.
(a)(b)Displacement fields for the modes (a)l=1 and (b) l=2, identified in Figure 4. l=1 resembles a surface wave and l=2 resembles a flexural plate mode slowly decaying in depth.
(a)(b)Magnitudes of particle velocity and stress fields in the cross-section for the modes (a)l=1 and (b) l=2, identified in Figure 4. The x- and y-dimensions are in microns. In both cases, T2, T4, and T6 are zero on the sides walls, thus satisfying boundary conditions of the plate, and T1, T5, and T6 are zero on the top surface, thus satisfying boundary conditions of the ridge.
(a)
Magnitudes of particle velocity (m/s) and stress fields (Pa) in the cross-section for model=1(b)
Magnitudes of particle velocity (m/s) and stress fields (Pa) in the cross-section for model=2The ridge modes found here can also be obtained using the general formulation as a single expansion constituting all 6 plate modes with boundary conditions applied aty0=0 and p=0,1 (see Table 1). This results in a 6x6 square A matrix, whose determinant has zeros at the same values of kz as the ones found using the symmetric and antisymmetric formulations. In addition to these ridge solutions, some spurious zeros are observed that correspond to values of kz where the plate modes become degenerate, which happens when the curves of Figure 3 intersect. Such false zeros not representing a solution to the ridge can be rejected by verifying the boundary conditions from field plots similar to Figure 6.Earlier it was mentioned that, by increasing the number of plate modes included in the series, the accuracy of the ridge solution can be improved. This point has been illustrated for the current example through a simple convergence study, as presented in Figure7. It can be seen that both the symmetric and the antisymmetric solutions of the ridge indeed converge as the number of plate modes, N, is increased. The boundary conditions in each case are applied at various points in the range 0≤|y0|<wR/2 and p=0,1, depending on the number of plate modes chosen.A convergence study to illustrate that both the symmetric and the antisymmetric solutions of the ridge converge as the number of plate modes included in the series are increased.
(a)
Convergence study of the symmetric ridge mode(b)
Convergence study of the antisymmetric ridge mode


## 4. Numerical Results and Validation

### 4.1. Verification with ANSYS

Finite element simulations were done in ANSYS to verify the calculated mode propagation velocities. In order to find the mode velocity in ANSYS, a modal analysis was performed for an anisotropic block of Silicon, as shown in Figure8. The block acts as a resonator imitating a section of the ridge waveguide. It has length lR along the z-direction and the same width wR and orientation as the ridge waveguide described in Section 3.3. The idea is to identify and study the resonant modes that resemble standing wave patterns resulting from propagating modes on the ridge waveguide. The vertical size of the block hRin the x-direction is chosen by trial and error such that the top surface bound modes decay to an acceptable level in depth. In all the simulations, the block has wR=2 μm. For the symmetric case, the size of the block was varied from lR=2.5 μm to 25 μm with  hR=3lR. For the antisymmetric case, the size of the block was varied from lR=2.5 μm to 7.5 μm with hR=10lR. Meshing was done using the SOLID185 element with a maximum element size of 0.5 μm. Thus, the number of elements generated in these simulations varied from 300 to 30000, depending on the size of the block. Boundary conditions imposed on the block are summarized in (27). Additional boundary conditions for the symmetric and antisymmetric cases are given in (28) and (29), respectively. Modes with half-wave resonance along z can be identified for the symmetric or antisymmetric cases and the resonant frequency fR for each mode is recorded. Plots of these modes obtained in ANSYS are shown in Figure 9. The propagation velocity of the constituent travelling wave modes can be estimated as, cR=2lRfR. This process is repeated for blocks of different lengths to obtain a dispersion curve over the frequency range of interest, as shown in Figure 10. Dispersion curves obtained in MATLAB by automatically tracking the l=1 and l=2 modes of Section 3.3 directly over ω, as described in Section 3.1, are also plotted. The ANSYS simulations and MATLAB computations match quite well within the frequency range of interest:
(27)vx=0,atx=-hR,vz=0,atz=-lR2,lR2,(28)symmetric⟶vy=0aty=0,(29)antisymmetric⟶vx=vz=0,aty=0,antisymmetric⟶vy=0,aty=0,z=0.Figure 8
Geometry and coordinate system of an anisotropic block of silicon simulated in ANSYS. The block imitates a resonating section of the ridge waveguide.Figure 9
Results of modal analysis in ANSYS showing the standing wave pattern of the symmetric ridge mode on left and the antisymmetric ridge mode on right. In both cases, the fields span half a wavelength along the length of the block. The color represents displacement in thex- and y-direction for the symmetric and antisymmetric case, respectively. Entire blocks are not shown because their heights are much larger than their lengths.Comparison of modal dispersion curves obtained for the (a) symmetric and (b) antisymmetric ridge modes shows good agreement between the ANSYS simulations and MATLAB computations. Note that they-axes in these plots have different scales.
(a)
Modal dispersion curve for the fundamental symmetric mode of a 2μm wide ridge in Si(b)
Modal dispersion curve for the fundamental antisymmetric mode of a 2μm wide ridge in Si


### 4.2. Example Cases

Results for a ridge etched in Silicon with various crystal orientations are presented in this section. Only the fundamental symmetric and antisymmetric modes of the ridge are considered. Dispersion curves for the ridge etched in a <100> Silicon wafer withx-axis along the <100> crystal axis and z-axis along the <001> crystal axis have already been presented in previous sections. Since the solutions only depend on the electrical dimensions of the ridge at a given frequency, the x-axis in the figures of this section has been generalized to the product of ridge width and frequency, having the units of velocity. Results for ridges aligned to the remaining symmetry planes of Silicon (x-axis along the wafer surface normal) are shown in Figure 11. The cases correspond to orientations with z-axis at 45 degrees from the <001> crystal axis on a <100> Silicon wafer, with z-axis at 0 and 90 degrees from the <001> crystal axis on a <110> Silicon wafer and with z-axis along the <112> crystal axis on a <111> Silicon wafer. Varying numbers of plate modes and an appropriate set of ridge boundary conditions, as mentioned in Section 2.3, were used in each case to construct the ridge solutions. The antisymmetric mode was found to be quite dispersive, whereas the symmetric mode does not exhibit a strong frequency dependence. This can be explained in view of the similarity of the symmetric mode to a surface wave. For surface waves, the propagation velocities are independent of frequency. Plots of Figure 11 also show that the variations of the velocity with respect to the angle of the ridge appear to be larger in the case of symmetric modes, where the ridge aligns with a plane of symmetry of the crystal. This dependence can be further investigated by using the general formulation and performing a continuous sweep of the ridge angle.Modal dispersion curves obtained for the fundamental (a) symmetric and (b) antisymmetric modes in a ridge waveguide oriented at various angles in a Silicon crystal. The velocity of the symmetric mode does not vary with frequency significantly, because this mode closely resembles the dispersionless surface wave.
(a)
Modal dispersion curves for the fundamental symmetric mode of ridge in Si oriented at various angles(b)
Modal dispersion curves for the fundamental antisymmetric mode of ridge in Si oriented at various anglesMode velocities as a function of ridge angle are presented in Figure12 for a ridge width and frequency product of 700. This corresponds to the behavior of a 2-micron wide ridge at 350 MHz. Note that the symmetric and antisymmetric are only meaningful when the ridge lies in a plane of mirror symmetry. For all intermediate angles, the ridge modes are formed by tracking the same plate modes used to construct the solutions at the symmetry planes, as a function of ridge angle. Therefore, these intermediate modes can only be classified as quasisymmetric and quasi-antisymmetric. A 〈100〉 Silicon wafer has 2 different mirror symmetry planes at 0 and 45 degrees from the <001> axis and rotation symmetry in 90 degrees, in the plane of the wafer. A 〈111〉 Silicon wafer has similar mirror symmetry planes at 0 and 60 degrees, but pointed in opposite directions, and a rotation symmetry in 120 degrees, in the plane of the wafer. Since the homogenous wave equation for uniform crystals gives the same mode velocity in the +z or –z direction, the plots are expected to be symmetric around 30 degrees for the 〈111〉 case. The plots of Figure 12 also confirm the stronger angular dependence of the quasisymmetric mode (~4% variation versus 1.5% for the quasi-antisymmetric case). The dependence on angle is smaller for the ridge waveguides in 〈111〉 wafer and is limited to <2%.Mode behavior for the fundamental (a) quasisymmetric (l=1) and (b) quasi-antisymmetric (l=2) ridge modes as a function of scan angle between planes of mirror symmetry for different Silicon wafers. Variations in velocity are quite subtle, suggesting that anisotropy may not be a significant issue while designing ridge waveguides in Silicon.
(a)
Modal propagation velocity as a function of scan angle for the fundamental quasi-symmetric mode of a ridge in Si withwR  ×  frequency = 700(b)
Modal propagation velocity as a function of scan angle for the fundamental quasi-antisymmetric mode of a ridge in Si withwR  ×  frequency = 700


## 4.1. Verification with ANSYS

Finite element simulations were done in ANSYS to verify the calculated mode propagation velocities. In order to find the mode velocity in ANSYS, a modal analysis was performed for an anisotropic block of Silicon, as shown in Figure8. The block acts as a resonator imitating a section of the ridge waveguide. It has length lR along the z-direction and the same width wR and orientation as the ridge waveguide described in Section 3.3. The idea is to identify and study the resonant modes that resemble standing wave patterns resulting from propagating modes on the ridge waveguide. The vertical size of the block hRin the x-direction is chosen by trial and error such that the top surface bound modes decay to an acceptable level in depth. In all the simulations, the block has wR=2 μm. For the symmetric case, the size of the block was varied from lR=2.5 μm to 25 μm with  hR=3lR. For the antisymmetric case, the size of the block was varied from lR=2.5 μm to 7.5 μm with hR=10lR. Meshing was done using the SOLID185 element with a maximum element size of 0.5 μm. Thus, the number of elements generated in these simulations varied from 300 to 30000, depending on the size of the block. Boundary conditions imposed on the block are summarized in (27). Additional boundary conditions for the symmetric and antisymmetric cases are given in (28) and (29), respectively. Modes with half-wave resonance along z can be identified for the symmetric or antisymmetric cases and the resonant frequency fR for each mode is recorded. Plots of these modes obtained in ANSYS are shown in Figure 9. The propagation velocity of the constituent travelling wave modes can be estimated as, cR=2lRfR. This process is repeated for blocks of different lengths to obtain a dispersion curve over the frequency range of interest, as shown in Figure 10. Dispersion curves obtained in MATLAB by automatically tracking the l=1 and l=2 modes of Section 3.3 directly over ω, as described in Section 3.1, are also plotted. The ANSYS simulations and MATLAB computations match quite well within the frequency range of interest:
(27)vx=0,atx=-hR,vz=0,atz=-lR2,lR2,(28)symmetric⟶vy=0aty=0,(29)antisymmetric⟶vx=vz=0,aty=0,antisymmetric⟶vy=0,aty=0,z=0.Figure 8
Geometry and coordinate system of an anisotropic block of silicon simulated in ANSYS. The block imitates a resonating section of the ridge waveguide.Figure 9
Results of modal analysis in ANSYS showing the standing wave pattern of the symmetric ridge mode on left and the antisymmetric ridge mode on right. In both cases, the fields span half a wavelength along the length of the block. The color represents displacement in thex- and y-direction for the symmetric and antisymmetric case, respectively. Entire blocks are not shown because their heights are much larger than their lengths.Comparison of modal dispersion curves obtained for the (a) symmetric and (b) antisymmetric ridge modes shows good agreement between the ANSYS simulations and MATLAB computations. Note that they-axes in these plots have different scales.
(a)
Modal dispersion curve for the fundamental symmetric mode of a 2μm wide ridge in Si(b)
Modal dispersion curve for the fundamental antisymmetric mode of a 2μm wide ridge in Si


## 4.2. Example Cases

Results for a ridge etched in Silicon with various crystal orientations are presented in this section. Only the fundamental symmetric and antisymmetric modes of the ridge are considered. Dispersion curves for the ridge etched in a <100> Silicon wafer withx-axis along the <100> crystal axis and z-axis along the <001> crystal axis have already been presented in previous sections. Since the solutions only depend on the electrical dimensions of the ridge at a given frequency, the x-axis in the figures of this section has been generalized to the product of ridge width and frequency, having the units of velocity. Results for ridges aligned to the remaining symmetry planes of Silicon (x-axis along the wafer surface normal) are shown in Figure 11. The cases correspond to orientations with z-axis at 45 degrees from the <001> crystal axis on a <100> Silicon wafer, with z-axis at 0 and 90 degrees from the <001> crystal axis on a <110> Silicon wafer and with z-axis along the <112> crystal axis on a <111> Silicon wafer. Varying numbers of plate modes and an appropriate set of ridge boundary conditions, as mentioned in Section 2.3, were used in each case to construct the ridge solutions. The antisymmetric mode was found to be quite dispersive, whereas the symmetric mode does not exhibit a strong frequency dependence. This can be explained in view of the similarity of the symmetric mode to a surface wave. For surface waves, the propagation velocities are independent of frequency. Plots of Figure 11 also show that the variations of the velocity with respect to the angle of the ridge appear to be larger in the case of symmetric modes, where the ridge aligns with a plane of symmetry of the crystal. This dependence can be further investigated by using the general formulation and performing a continuous sweep of the ridge angle.Modal dispersion curves obtained for the fundamental (a) symmetric and (b) antisymmetric modes in a ridge waveguide oriented at various angles in a Silicon crystal. The velocity of the symmetric mode does not vary with frequency significantly, because this mode closely resembles the dispersionless surface wave.
(a)
Modal dispersion curves for the fundamental symmetric mode of ridge in Si oriented at various angles(b)
Modal dispersion curves for the fundamental antisymmetric mode of ridge in Si oriented at various anglesMode velocities as a function of ridge angle are presented in Figure12 for a ridge width and frequency product of 700. This corresponds to the behavior of a 2-micron wide ridge at 350 MHz. Note that the symmetric and antisymmetric are only meaningful when the ridge lies in a plane of mirror symmetry. For all intermediate angles, the ridge modes are formed by tracking the same plate modes used to construct the solutions at the symmetry planes, as a function of ridge angle. Therefore, these intermediate modes can only be classified as quasisymmetric and quasi-antisymmetric. A 〈100〉 Silicon wafer has 2 different mirror symmetry planes at 0 and 45 degrees from the <001> axis and rotation symmetry in 90 degrees, in the plane of the wafer. A 〈111〉 Silicon wafer has similar mirror symmetry planes at 0 and 60 degrees, but pointed in opposite directions, and a rotation symmetry in 120 degrees, in the plane of the wafer. Since the homogenous wave equation for uniform crystals gives the same mode velocity in the +z or –z direction, the plots are expected to be symmetric around 30 degrees for the 〈111〉 case. The plots of Figure 12 also confirm the stronger angular dependence of the quasisymmetric mode (~4% variation versus 1.5% for the quasi-antisymmetric case). The dependence on angle is smaller for the ridge waveguides in 〈111〉 wafer and is limited to <2%.Mode behavior for the fundamental (a) quasisymmetric (l=1) and (b) quasi-antisymmetric (l=2) ridge modes as a function of scan angle between planes of mirror symmetry for different Silicon wafers. Variations in velocity are quite subtle, suggesting that anisotropy may not be a significant issue while designing ridge waveguides in Silicon.
(a)
Modal propagation velocity as a function of scan angle for the fundamental quasi-symmetric mode of a ridge in Si withwR  ×  frequency = 700(b)
Modal propagation velocity as a function of scan angle for the fundamental quasi-antisymmetric mode of a ridge in Si withwR  ×  frequency = 700


## 5. Conclusion

A semianalytical method for finding the approximate solution of the elastic equations in a ridge waveguide made of an anisotropic crystal was presented. The proposed formulation is general and leads to the propagating modes of the waveguide for the arbitrary orientation of the crystal and can be simplified to give symmetric and antisymmetric modes,s when the ridge is parallel with a plane of crystal symmetry. Numerical issues involving the tracking of the plate modes were discussed. The method has been used for the modal analysis of deep ridges etched in Silicon. Some of these results have been verified by finite element simulations in ANSYS. For the case of Silicon, the modal behavior appears to be quite similar along various crystal orientations, suggesting that the effect of the anisotropy may be neglected in a first-order analysis.



---

*Source: 685326-2012-04-17.xml*