Elastic Modes of an Anisotropic Ridge Waveguide

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.

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 𝑤_𝑅 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 𝑥 direction, and propagation is along 𝑧-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 |𝑦|≤𝑤_𝑅/2, 𝑥≤0 region.

(a)

(b)

(c)

Figure 1 

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

Section 2 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 superscript 𝑇.

2.1. Elastic Wave Equation

The homogenous time harmonic acoustic field equations are [20]∇⋅𝑇⃗=𝑗𝜌𝜔𝑣,⃗𝐜∇_s𝑣⃗=𝑗𝜔𝑇⃗, (1) where, 𝜔 is the frequency, 𝜌 is the material density, 𝑣⃗=[𝑣_𝑥,𝑣_𝑦,𝑣_𝑧]^𝑇=[𝑣₁,𝑣₂,𝑣₃]^𝑇 is the particle velocity field, 𝑇⃗=[𝑇_𝑥𝑥,𝑇_𝑦𝑦,𝑇_𝑧𝑧,𝑇_𝑦𝑧,𝑇_𝑥𝑧,𝑇_𝑥𝑦]^𝑇=[𝑇₁,𝑇₂,𝑇₃,𝑇₄,𝑇₅,𝑇₆]^𝑇is the abbreviated notation for the second-order stress tensor:𝑇𝐓=¹𝑇₆𝑇₅𝑇₆𝑇₂𝑇₄𝑇₅𝑇₄𝑇₃, (2) and 𝐜 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, 𝐜^𝐸, to the ridge coordinate axes using the bond transformation matrix, 𝐌_𝑏, given in [20]:𝐜=𝐌_𝑏𝐜^𝐸𝐌^𝑇_𝑏. (3) In a general triclinic crystal, which only has center point symmetry and no planes of mirror symmetry, 𝐜^𝐸has 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:∇⋅𝐜∇_s+𝜌𝜔²⃗𝑗𝑣⃗=0,𝑇=−𝐜∇_s𝑣⃗𝜔. (4)(5) For uniform plane waves of the form 𝑣⃗^𝑢𝑝𝑤=𝑒^{−𝑗(𝑘𝑧𝑧+𝑘𝑥𝑥+𝑘𝑦𝑦)}𝑢, 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):Γ𝑢=0. (6) Here, 𝑢=[𝑢_𝑥,𝑢_𝑦,𝑢_𝑧]^𝑇 is the velocity polarization, and 𝑘⃗=[𝑘_𝑥,𝑘_𝑦,𝑘_𝑧]^𝑇 is the propagation vector with 𝑘=|𝑘⃗|. The solutions are obtained by setting the characteristic determinant equal to zero. For a given material and frequency, the surfaces in 𝑘⃗-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, 𝑘_𝑧 is the propagation direction and has to be real, while 𝑘_𝑥 must be complex to allow for decay in the transverse direction. By setting 𝑘_𝑧∈ℝ and 𝑘_𝑥∈ℂ, a 6th degree polynomial in 𝑘_𝑦 is obtained, whose roots 𝑘_𝑦𝑚, in general, cannot be found analytically as stated by the Abel-Ruffini theorem [21] and must be found numerically as functions of 𝑘_𝑥 and 𝑘_𝑧 as 𝑘_𝑦𝑚=𝑘_𝑦𝑚(𝑘_𝑥,𝑘_𝑧). If 𝑘_𝑦𝑚 is not purely real, then it appears as a pair with its conjugate 𝑘^∗_𝑦𝑚. If the crystal has mirror symmetry in the x-z plane, and if 𝑘_𝑦𝑚 is not purely imaginary, then it also appears as a pair with −𝑘_𝑦𝑚. For every 𝑘_𝑦𝑚, the corresponding 𝑢_𝑚 is the eigenvector corresponding to the zero eigenvalue. In practice, where 𝑘_𝑦𝑚’s are found numerically, 𝑢_𝑚 may be found as the eigenvector corresponding to the smallest eigenvalue 𝜆_Γ as described in (8):detΓ𝑘_𝑥,𝑘_𝑦𝑚,𝑘_𝑧Γ𝑘=0,_𝑥,𝑘_𝑦𝑚,𝑘_𝑧𝑢_𝑚=min𝜆_Γ𝑢_𝑚,𝑚=1⋯6. (7)(8)

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):𝑣⃗^plate=⁶_𝑚=1𝑏_𝑚𝑣⃗_𝑚^𝑢𝑝𝑤. (9) If the plate is parallel to a plane of mirror symmetry of the crystal, then the roots of (7) appear in pairs of ±𝑘_𝑦𝑚 with polarizations [𝑢_𝑥,±𝑢_𝑦,𝑢_𝑧]^𝑇. 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):𝑣⃗^plate-s=𝑒^{−𝑗𝑘𝑧𝑧}𝑒^{−𝑗𝑘𝑥𝑥3}_𝑚=1𝑏^𝑠_𝑚cos𝑘_𝑦𝑚𝑦𝑢_𝑥𝑚−𝑗sin𝑘_𝑦𝑚𝑦𝑢_𝑦𝑚cos𝑘_𝑦𝑚𝑦𝑢_𝑧𝑚, (10)𝑣⃗^plate-as=𝑒^{−𝑗𝑘𝑧𝑧}𝑒^{−𝑗𝑘𝑥𝑥3}_𝑚=1𝑏_𝑚^𝑎𝑠−𝑗sin𝑘_𝑦𝑚𝑦𝑢_𝑥𝑚cos𝑘_𝑦𝑚𝑦𝑢_𝑦𝑚−𝑗sin𝑘_𝑦𝑚𝑦𝑢_𝑧𝑚. (11) 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 𝑛:𝐓^plate𝑛|_{atfreeboundary}=0. (12) The 2 free surfaces of the plate give rise to a matrix equation in 6 unknown coefficients represented by 𝑏̂=[𝑏₁,𝑏₂,𝑏₃,𝑏₄,𝑏₅,𝑏₆]^𝑇. Without loss of generality, the boundary conditions are applied at 𝑥=0, and 𝑧=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 𝑏̂^s/as=[𝑏₁,𝑏₂,𝑏₃]^𝑇 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 𝑠 and antisymmetric cases have superscript as:𝐁̂𝐓𝑏=^plate𝑦|_{𝑥=0,𝑦=𝑤𝑅/2,𝑧=0}−𝐓^plate𝑦|_{𝑥=0,𝑦=−𝑤𝑅/2,𝑧=0}=0, (13)𝐁^s/as𝑏̂^s/as=𝐓^plate-s/as𝑦|_{𝑥=0,𝑦=𝑤𝑅/2,𝑧=0}=0. (14) 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 𝑘_𝑥−𝑘_𝑧 plane and search for the zeros of the determinant (note that 𝑘_𝑦𝑚𝑛=𝑘_𝑦𝑚𝑛(𝑘_𝑥𝑛,𝑘_𝑧) according to (7)). For the infinite plate, the search is conducted for real values of 𝑘_𝑥 and 𝑘_𝑧, and results in closed contour in the 𝑘_𝑥−𝑘_𝑧 plane, or 𝑘_𝑥𝑛=𝑘_𝑥𝑛(𝑘_𝑧), 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 𝑘_𝑥 in the upper half of the complex plane, Im(𝑘_𝑥≥0), which corresponds to the edge-bound solutions. If 𝑘_𝑥𝑛 is not purely real, then it appears as a pair with its conjugate 𝑘^∗_𝑥𝑛. If the crystal has mirror symmetry in the y-z plane, and if 𝑘_𝑥𝑛 is not purely imaginary, then it also appears as a pair with −𝑘_𝑥𝑛. For every 𝑘_𝑥𝑛, the corresponding coefficient vector 𝑏̂_𝑛, theoretically representing the nullspace, is approximated as the eigenvector corresponding to the smallest eigenvalue 𝜆_𝐁 as described in (16):det𝐁𝑘_𝑥,𝑘_𝑧𝐁𝑘=0,_𝑥𝑛𝑏̂_𝑛=min𝜆_B𝑏̂_𝑛,∀𝑛. (15)(16) In Section 3, issues involved with locating and tracking the zeros of det(𝐁(𝑘_𝑥,𝑘_𝑧)) or the poles of log(|det(𝐁(𝑘_𝑥,𝑘_𝑧))|) are discussed.

2.3. Anisotropic Semi-Infinite Plate or Ridge

General solutions to the ridge can be constructed as a linear combination of 𝑁 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(𝑘_𝑥≥0) contribute to edge-bound solutions. If the crystal has mirror symmetry with respect to the y-z plane, and if 𝑘_𝑥𝑛 is not purely imaginary, then the modes should be chosen in pairs of [𝑘_𝑥𝑛,−𝑘^∗_𝑥𝑛] 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 𝑁. A simple convergence study to illustrate this has been presented in Section 3.3:𝑣⃗^ridge-(s/as)=^𝑁_𝑛=1𝑎_𝑛𝑣^plate-(s/as)_𝑛⃗. (17) Without loss of generality, the boundary condition is applied at 𝑧=0 on the top surface of the ridge to yield the matrix equation (18), which should be satisfied over the range |𝑦|≤𝑤_𝑅/2. Since the boundary conditions are a function of 𝑦, 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 𝐀 corresponding to the field quantities that must vanish on the top surface of the ridge and (2) finding the value of 𝑘_𝑧 that makes this matrix singular and the corresponding nullspace 𝑎. One way of expressing the boundary conditions is by expanding them in terms of appropriate domain functions over the range of 𝑦. 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 𝑦₀. For the case of free boundary conditions, setting the traction component and its higher order derivatives at a fixed point in 𝑦 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 𝑎=[𝑎₁,…,𝑎_𝑁]^𝑇 is the coefficient vector weighing each plate mode:𝐓^ridge𝑥|_{𝑥=0,|𝑦|≤𝑤𝑅/2,𝑧=0}=0, (18)𝑑^𝑝𝑑𝑦^𝑝𝐓^ridge𝑥∣_{𝑥=0,𝑦=𝑦0,𝑧=0}=𝐴⃗₁𝑝,𝑦₀|⋯|𝐴⃗_𝑁𝑝,𝑦₀𝑎=0;𝑝=0,1,…,𝑃. (19) For 0<|𝑦₀|<𝑤_𝑅/2, (19) gives 3 nontrivial equations for each value of 𝑝. For 𝑦=𝑤_𝑅/2 or −𝑤_𝑅/2, (19) only gives 2 nontrivial equations for 𝑝=0 because one traction component 𝑇_𝑥𝑦 is already zero by the plate boundary conditions. For 𝑦₀=0, and the symmetric case, 𝑇_𝑥𝑦, and it even, derivatives are zero at 𝑦₀, and (19) produces only 2 non-trivial equations for the even values of 𝑝. Similarly, the odd derivatives of 𝑇_𝑥𝑥 and 𝑇_𝑥𝑧 are zero, yielding only 1 non-trivial equation for every odd value of 𝑝. For 𝑦₀=0 and antisymmetric case, (19) gives 1 and 2 non-trivial equations for each even or odd value of 𝑝, respectively. The number of equations obtained for the 𝑝th derivative for different choices of 𝑦₀ 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:𝐀𝑎=0. (20) In order to have a square boundary condition matrix 𝐀^(s/as), the values of 𝑝 much be appropriately chosen such that the total number of equations is 𝑁.

Table 1 

Number of ridge boundary conditions obtained at each point on the 𝑦-axis for the 𝑝th derivative for the general, symmetric, and the antisymmetric cases.

After constructing a suitable 𝐀 matrix encapsulating the boundary conditions of the ridge, the next problem is to find its nullspace 𝑎. 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 ℂ^𝑁 of 𝑎. If the matrix is square, then the characteristic determinant in (21) is set to zero, and the values 𝑘_𝑧𝑙 corresponding to the propagating modes of the ridge are found numerically. The edge bound modes only correspond to those values of 𝑘_𝑧𝑙 for which none of the contributing 𝑘_𝑥𝑛𝑙’s with nonzero coefficients are purely real. For every 𝑘_𝑧𝑙, the respective nullspace 𝑎_𝑙 representing the coefficient vector is approximated as the eigenvector corresponding to the smallest eigenvalue 𝜆_𝐀 as described in (22). As functions of frequency, the dispersion relations of the ridge modes can be expressed as 𝑘_𝑧𝑙=𝑘_𝑧𝑙(𝜔), 𝑘_𝑥𝑛𝑙=𝑘_𝑥𝑛𝑙(𝑘_𝑧𝑙,𝜔) and, 𝑘_{𝑦𝑚𝑛𝑙}=𝑘_{𝑦𝑚𝑛𝑙}(𝑘_𝑥𝑛𝑙,𝑘_𝑧𝑙,𝜔):det𝐀𝑘_𝑧=0, (21)𝐀𝑘_𝑧𝑙𝑎_𝑙=min𝜆_A𝑎_𝑙,∀𝑙. (22) 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 𝑘_𝑧𝑙. If 𝑘^𝑁_𝑧𝑙 is the propagation constant obtained using 𝑁 plate modes, then lim_𝑁→∞𝑘^𝑁_𝑧𝑙=𝑘_𝑧𝑙.

3. Numerical Calculations

3.1. Pole Tracking

As mentioned earlier, the solutions of (15), 𝑘_𝑥𝑛(𝑘_𝑧) cannot be found analytically; hence, the zeros of the determinant of 𝐁 or the poles of log(|det(𝐁)|) must be found numerically. Performing a raster scan in the complex 𝑘_𝑥 plane for each 𝑘_𝑧to 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 𝑘_𝑧 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 𝑘_𝑧 to ensure the continuity of the determinant. Using a polynomial fit, the approximate location of poles or 𝑘_𝑥𝑛’s for each mode are estimated for successive 𝑘_𝑧’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 𝑘_𝑧 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 over 𝑘_𝑧 to get 𝑘_𝑥𝑛(𝑘_𝑧), the appropriate boundary conditions (19) for the ridge top surface are applied on them to construct the 𝐀 matrix. The zeros of the determinant of 𝐀 or the poles of log(|det(𝐀)|) represent the propagation constants of the ridge modes. Locations of these poles are 𝑘_𝑧𝑙 with the constituent plate modes, having 𝑘_𝑥=𝑘_𝑥𝑛𝑙(𝑘_𝑧𝑙). 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 𝑘_𝑥𝑛(𝑘_𝑧), 𝑘_𝑧𝑙(𝜔), and 𝑘_𝑥𝑛𝑙(𝑘_𝑧,𝜔), or 𝑘_𝑧𝑙(𝐜) and 𝑘_𝑥𝑛𝑙(𝑘_𝑧,𝐜), are tracked using a combination of gradient descent and direct line search.

3.2. Spurious Poles

When log(|det(𝐁)|) 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 𝑘_𝑦𝑚’s become redundant with the same velocity polarization vectors. Such degeneracies occur when two 𝑘_𝑦𝑚’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 𝐁 in (23), is defined for the general case with 𝑚=1⋯6. Correspondingly, a 3×3 matrix 𝐁^(s/as) in (24) and (25), is defined for the symmetric and antisymmetric cases, respectively, with 𝑚=1…3. The determinant of the main plate boundary condition matrix 𝐁 or 𝐁^s or 𝐁^as is divided by the determinant of this pole nullifying matrix 𝐁 or 𝐁^s or 𝐁^as, respectively, to cancel out the spurious poles in the logarithm of their absolute values, as shown in (26). The pole tracking algorithms use 𝑑(𝑘_𝑥,𝑘_𝑧) while tracking the plate modes:B∶𝐁_1𝑚,𝐁_2𝑚,𝐁_3𝑚,𝐁_4𝑚,𝐁_5𝑚,𝐁_6𝑚^𝑇≜𝑢_𝑥𝑚,𝑢_𝑦𝑚,𝑢_𝑧𝑚,𝑘_𝑦𝑚𝑢_𝑥𝑚,𝑘_𝑦𝑚𝑢_𝑦𝑚,𝑘_𝑦𝑚𝑢_𝑧𝑚^𝑇, (23)𝐁^s∶𝐁^s_1𝑚,𝐁^s_2𝑚,𝐁^s_3𝑚^𝑇≜𝑢_𝑥𝑚,𝑘_𝑦𝑚𝑢_𝑦𝑚,𝑢_𝑧𝑚^𝑇, (24)𝐁^as∶𝐁^as_1𝑚,𝐁^as_2𝑚,𝐁^as_3𝑚^𝑇≜𝑘_𝑦𝑚𝑢_𝑥𝑚,𝑢_𝑦𝑚,𝑘_𝑦𝑚𝑢_𝑧𝑚^𝑇, (25)𝑑𝑘_𝑥,𝑘_𝑧𝐁=logdet^(s/as)𝑘_𝑥,𝑘_𝑧det𝐁^(s/as)𝑘_𝑥,𝑘_𝑧. (26)

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 𝑧-axis along the ⟨001⟩ crystal axis. Silicon is a cubic crystal with a density 𝜌=2300kg/m³ and 3 independent stiffness constants: 𝑐₁₁=165.7GPa, 𝑐₁₂=63.9GPa, and 𝑐₄₄=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 complex 𝑘_𝑥 plane for 𝑘_𝑧=0, as shown in Figure 2. Plate boundary conditions for the general case have been used in this calculation; hence, 𝑚=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 𝑘_sh=𝜔𝜌/𝑐₄₄. The plate modes are then tracked over 𝑘_𝑧 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 𝑛=1, 𝑛=2, and 𝑛=3 plate modes are purely propagating (real 𝑘_𝑥) for small values of 𝑘_𝑧/𝑘_sh 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 𝑘_𝑥-axis. At higher values of 𝑘_𝑧, the 𝑛=1 and 𝑛=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 𝑘_𝑧.

Figure 2 

Plate modes of a 2-micron wide plate with 𝑥-axis along the ⟨001⟩ crystal axis and 𝑧-axis along the ⟨001⟩ crystal axis at 350 MHz for 𝑘𝑧=0.

Figure 3 

Plate modes of Figure 2 tracked over 𝑘𝑧 at 350 MHz. The 𝑛=1, 𝑛=2, and 𝑛=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 with 𝑘_𝑥𝑛 (𝑛=1⋯6) are used in the expansion, and the top boundary conditions are applied at 𝑦₀=0 and 𝑝=0,1 (see Table 1), separately, for the symmetric and antisymmetric cases. This results in 3 × 3 square 𝐀^s and 𝐀^as 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 𝑙=1 and 𝑙=2. Each zero corresponds to a propagating ridge mode. Values of 𝑘_𝑧𝑙 and the 𝑘_𝑥𝑛𝑙’s of the 6 constituent plate modes along with the magnitude of their coefficients 𝑎_𝑛𝑙 are listed in Table 2. It can be seen that the main contribution to 𝑙=1 is from the 2 fundamental symmetric plate modes =1, and 𝑛=2, both of which have complex 𝑘_𝑥𝑛’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 𝑘_𝑧1 is almost equal to the propagation constant of a z-directed surface wave. For the case of 𝑙=2, the main contribution is from the fundamental antisymmetric plate mode 𝑛=3, which has a small but imaginary 𝑘_𝑥33. 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 𝑘_𝑥11 and 𝑘_𝑥21 with 𝑘_𝑥33, it is seen that 𝑙=2 decays much slower compared to 𝑙=1. For a ridge waveguide standing on a substrate, in order to minimize reflections from the substrate, a higher aspect ratio (ℎ_𝑅/𝑤_𝑅 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 𝑙=1 and 𝑙=2 are shown in Figure 6. In both cases, 𝑇₂, 𝑇₄, and 𝑇₆ should be zero on the sides walls in order to satisfy boundary conditions of the plate, and 𝑇₁, 𝑇₅, and 𝑇₆ should be zero on the top surface in order to satisfy boundary conditions of the ridge. In the present example, the values obtained for 𝑘_{𝑦𝑚𝑛𝑙} 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 𝑦 on the top surface to a good extent.

Table 2 

Constituent plate mode coefficients and propagation constants for the two cases identified in Figure 4. 𝑙=1 is classified as the fundamental symmetric ridge mode and 𝑙=2 is classified as the fundamental antisymmetric ridge mode.

(a)

(b)

Figure 4 

Determinant of the matrix of ridge boundary conditions separately applied to (a) 3 symmetric and (b) 3 antisymmetric plate modes tracked in Figure 3. Two poles are identified with 𝑙=1 and 𝑙=2.

(a)

(b)

Figure 5 

Displacement fields for the modes (a) 𝑙=1 and (b) 𝑙=2, identified in Figure 4. 𝑙=1 resembles a surface wave and 𝑙=2 resembles a flexural plate mode slowly decaying in depth.

(a) Magnitudes of particle velocity (m/s) and stress fields (Pa) in the cross-section for mode 𝑙=1

(b) Magnitudes of particle velocity (m/s) and stress fields (Pa) in the cross-section for mode 𝑙=2

Figure 6 

Magnitudes of particle velocity and stress fields in the cross-section for the modes (a) 𝑙=1 and (b) 𝑙=2, identified in Figure 4. The 𝑥- and 𝑦-dimensions are in microns. In both cases, 𝑇2, 𝑇4, and 𝑇6 are zero on the sides walls, thus satisfying boundary conditions of the plate, and 𝑇1, 𝑇5, and 𝑇6 are zero on the top surface, thus satisfying boundary conditions of the ridge.

The 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 at 𝑦₀=0 and 𝑝=0,1 (see Table 1). This results in a 6x6 square 𝐀 matrix, whose determinant has zeros at the same values of 𝑘_𝑧 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 𝑘_𝑧 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 Figure 7. 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≤|𝑦₀|<𝑤_𝑅/2 and 𝑝=0,1, depending on the number of plate modes chosen.

(a) Convergence study of the symmetric ridge mode

(b) Convergence study of the antisymmetric ridge mode

Figure 7 

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.

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 Figure 8. The block acts as a resonator imitating a section of the ridge waveguide. It has length 𝑙_𝑅 along the z-direction and the same width 𝑤_𝑅 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 ℎ_𝑅in 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 𝑤_𝑅=2 μm. For the symmetric case, the size of the block was varied from 𝑙_𝑅=2.5 μm to 25 μm with ℎ_𝑅=3𝑙_𝑅. For the antisymmetric case, the size of the block was varied from 𝑙_𝑅=2.5 μm to 7.5 μm with ℎ_𝑅=10𝑙_𝑅. 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 𝑧 can be identified for the symmetric or antisymmetric cases and the resonant frequency 𝑓_𝑅 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, 𝑐_𝑅=2𝑙_𝑅𝑓_𝑅. 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 𝑙=1 and 𝑙=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: 𝑣_𝑥=0,at𝑥=−ℎ_𝑅,𝑣_𝑧𝑙=0,at𝑧=−^𝑅2,𝑙_𝑅2,symmetric⟶𝑣_𝑦=0at𝑦=0,antisymmetric⟶𝑣_𝑥=𝑣_𝑧=0,at𝑦=0,antisymmetric⟶𝑣_𝑦=0,at𝑦=0,𝑧=0. (27)(28)(29)

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 the x- and y-direction for the symmetric and antisymmetric case, respectively. Entire blocks are not shown because their heights are much larger than their lengths.

(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

Figure 10 

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 the y-axes in these plots have different scales.

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 with 𝑥-axis along the <100> crystal axis and 𝑧-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 𝑥-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 𝑧-axis at 45 degrees from the <001> crystal axis on a <100> Silicon wafer, with 𝑧-axis at 0 and 90 degrees from the <001> crystal axis on a <110> Silicon wafer and with 𝑧-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.

(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 angles

Figure 11 

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.

Mode velocities as a function of ridge angle are presented in Figure 12 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%.

(a) Modal propagation velocity as a function of scan angle for the fundamental quasi-symmetric mode of a ridge in Si with 𝑤𝑅  ×  frequency = 700

(b) Modal propagation velocity as a function of scan angle for the fundamental quasi-antisymmetric mode of a ridge in Si with 𝑤𝑅  ×  frequency = 700

Figure 12 

Mode behavior for the fundamental (a) quasisymmetric (𝑙=1) and (b) quasi-antisymmetric (𝑙=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.

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.

Appendix

Elastic Wave Equation
𝜕∇⋅=𝜕𝜕𝑥000𝜕𝜕𝑧0𝜕𝜕𝑦0𝜕𝜕𝑦0𝜕𝜕𝑧𝜕𝜕𝑥00𝜕𝜕𝑧𝜕𝜕𝑦0𝜕𝑥, (A.1)∇_s=(∇⋅)^𝑇, (A.2)Γ=𝑘²Γ_ch−𝜌𝜔²𝐈, (A.3)Γ_ch=Γ_𝑘𝐜Γ^𝑇_𝑘𝑘², (A.4)Γ_𝑘=𝑘_𝑥000𝑘_𝑦000𝑘_𝑧0𝑘_𝑦𝑘_𝑧𝑘_𝑧0𝑘_𝑥𝑘_𝑦𝑘_𝑥0, (A.5)Γ₁₁=𝑘_𝑥𝑐₁₁𝑘_𝑥+𝑐₁₆𝑘_𝑦+𝑐₁₅𝑘_𝑧+𝑘_𝑦𝑐₁₆𝑘_𝑥+𝑐₆₆𝑘_𝑦+𝑐₅₆𝑘_𝑧+𝑘_𝑧𝑐₁₅𝑘_𝑥+𝑐₅₆𝑘_𝑦+𝑐₅₅𝑘_𝑧−𝜌𝜔²,Γ₁₂=𝑘_𝑥𝑐₁₆𝑘_𝑥+𝑐₆₆𝑘_𝑦+𝑐₅₆𝑘_𝑧+𝑘_𝑦𝑐₁₂𝑘_𝑥+𝑐₂₆𝑘_𝑦+𝑐₂₅𝑘_𝑧+𝑘_𝑧𝑐₁₄𝑘_𝑥+𝑐₄₆𝑘_𝑦+𝑐₄₅𝑘_𝑧,Γ₁₃=𝑘_𝑥𝑐₁₅𝑘_𝑥+𝑐₅₆𝑘_𝑦+𝑐₅₅𝑘_𝑧+𝑘_𝑦𝑐₁₄𝑘_𝑥+𝑐₄₆𝑘_𝑦+𝑐₄₅𝑘_𝑧+𝑘_𝑧𝑐₁₃𝑘_𝑥+𝑐₃₆𝑘_𝑦+𝑐₃₅𝑘_𝑧,Γ₂₁=𝑘_𝑥𝑐₁₆𝑘_𝑥+𝑐₁₂𝑘_𝑦+𝑐₁₄𝑘_𝑧+𝑘_𝑦𝑐₆₆𝑘_𝑥+𝑐₂₆𝑘_𝑦+𝑐₄₆𝑘_𝑧+𝑘_𝑧𝑐₅₆𝑘_𝑥+𝑐₂₅𝑘_𝑦+𝑐₄₅𝑘_𝑧,Γ₂₂=𝑘_𝑥𝑐₆₆𝑘_𝑥+𝑐₂₆𝑘_𝑦+𝑐₄₆𝑘_𝑧+𝑘_𝑦𝑐₂₆𝑘_𝑥+𝑐₂₂𝑘_𝑦+𝑐₂₄𝑘_𝑧+𝑘_𝑧𝑐₄₆𝑘_𝑥+𝑐₂₄𝑘_𝑦+𝑐₄₄𝑘_𝑧−𝜌𝜔²,Γ₂₃=𝑘_𝑥𝑐₅₆𝑘_𝑥+𝑐₂₅𝑘_𝑦+𝑐₄₅𝑘_𝑧+𝑘_𝑦𝑐₄₆𝑘_𝑥+𝑐₂₄𝑘_𝑦+𝑐₄₄𝑘_𝑧+𝑘_𝑧𝑐₃₆𝑘_𝑥+𝑐₂₃𝑘_𝑦+𝑐₃₄𝑘_𝑧,Γ₃₁=𝑘_𝑥𝑐₁₅𝑘_𝑥+𝑐₁₄𝑘_𝑦+𝑐₁₃𝑘_𝑧+𝑘_𝑦𝑐₅₆𝑘_𝑥+𝑐₄₆𝑘_𝑦+𝑐₃₆𝑘_𝑧+𝑘_𝑧𝑐₅₅𝑘_𝑥+𝑐₄₅𝑘_𝑦+𝑐₃₅𝑘_𝑧,Γ₃₂=𝑘_𝑥𝑐₅₆𝑘_𝑥+𝑐₄₆𝑘_𝑦+𝑐₃₆𝑘_𝑧+𝑘_𝑦𝑐₂₅𝑘_𝑥+𝑐₂₄𝑘_𝑦+𝑐₂₃𝑘_𝑧+𝑘_𝑧𝑐₄₅𝑘_𝑥+𝑐₄₄𝑘_𝑦+𝑐₃₄𝑘_𝑧,Γ₃₃=𝑘_𝑥𝑐₅₅𝑘_𝑥+𝑐₄₅𝑘_𝑦+𝑐₃₅𝑘_𝑧+𝑘_𝑦𝑐₄₅𝑘_𝑥+𝑐₄₄𝑘_𝑦+𝑐₃₄𝑘_𝑧+𝑘_𝑧𝑐₃₅𝑘_𝑥+𝑐₃₄𝑘_𝑦+𝑐₃₃𝑘_𝑧−𝜌𝜔². (A.6)

Plate Boundary Conditions (𝑚=1…6 for the General Case and 𝑚=1…3 for the Symmetric/Antisymmetric Cases)
𝐵⃗_𝑚=^{⃗𝐷1𝑚+⃗𝐷2𝑚𝑒𝑦𝑚𝑤𝑅−(𝑗𝑘/2)⃗𝐷1𝑚+⃗𝐷2𝑚𝑒𝑦𝑚𝑤𝑅+(𝑗𝑘/2)}𝜔,𝐵⃗^s_𝑚=𝑗𝐷⃗¹_𝑚sin𝑘_𝑦𝑚𝑤_𝑅/2−𝐷⃗²_𝑚cos𝑘_𝑦𝑚𝑤_𝑅/2𝜔𝐵⃗_as𝑚,=𝑗𝐷⃗²_𝑚sin𝑘_𝑦𝑚𝑤_𝑅/2−𝐷⃗¹_𝑚cos𝑘_𝑦𝑚𝑤_𝑅/2𝜔𝐷⃗_1𝑚,=𝑐₄₆𝑘_𝑦𝑚𝑢_𝑧𝑚+𝑘_𝑧𝑢_𝑦𝑚+𝑐₆₆𝑘_𝑥𝑢_𝑦𝑚+𝑘_𝑦𝑚𝑢_𝑥𝑚𝑐₂₄𝑘_𝑦𝑚𝑢_𝑧𝑚+𝑘_𝑧𝑢_𝑦𝑚+𝑐₂₆𝑘_𝑥𝑢_𝑧𝑚+𝑘_𝑦𝑚𝑢_𝑥𝑚𝑐₄₄𝑘_𝑦𝑚𝑢_𝑧𝑚+𝑘_𝑧𝑢_𝑦𝑚+𝑐₄₆𝑘_𝑥𝑢_𝑦𝑚+𝑘_𝑦𝑚𝑢_𝑥𝑚𝐷⃗_2𝑚,=𝑐₁₆𝑘_𝑥𝑢_𝑥𝑚+𝑐₂₆𝑘_𝑦𝑚𝑢_𝑦𝑚+𝑐₃₆𝑘_𝑧𝑢_𝑧𝑚+𝑐₅₆𝑘_𝑥𝑢_𝑧𝑚+𝑘_𝑧𝑢_𝑥𝑚𝑐₁₂𝑘_𝑥𝑢_𝑥𝑚+𝑐₂₂𝑘_𝑦𝑚𝑢_𝑦𝑚+𝑐₂₃𝑘_𝑧𝑢_𝑧𝑚+𝑐₂₅𝑘_𝑥𝑢_𝑧𝑚+𝑘_𝑧𝑢_𝑥𝑚𝑐₁₄𝑘_𝑥𝑢_𝑥𝑚+𝑐₂₄𝑘_𝑦𝑚𝑢_𝑦𝑚+𝑐₃₄𝑘_𝑧𝑢_𝑧𝑚+𝑐₄₅𝑘_𝑥𝑢_𝑧𝑚+𝑘_𝑧𝑢_𝑥𝑚. (A.7)(A.8)(A.9)(A.10)(A.11)

Ridge Boundary Conditions (𝑛=1⋯𝑁)
𝐴⃗_𝑛=∑𝑝,𝑦^6𝑚=1𝑏_𝑚𝑛−𝑗𝑘_𝑦𝑚𝑛^𝑝𝐶⃗¹_𝑚𝑛+𝐶⃗²_𝑚𝑛𝑒^{−𝑗𝑘𝑦𝑚𝑛𝑦}𝜔𝐴⃗_s𝑛,=∑𝑝,𝑦^3𝑚=1𝑏_𝑚𝑛𝑘_𝑦𝑚^𝑝𝑗𝐶⃗¹_𝑚𝑛sin𝑝𝜋/2+𝑘_𝑦𝑚𝑛𝑦−𝐶⃗²_𝑚𝑛cos𝑝𝜋/2+𝑘_𝑦𝑚𝑛𝑦𝜔𝐴⃗_as𝑛,=∑𝑝,𝑦^3𝑚=1𝑏_𝑚𝑛𝑘_𝑦𝑚^𝑝𝑗𝐶⃗²_𝑚𝑛sin𝑝𝜋/2+𝑘_𝑦𝑚𝑛𝑦−𝐶⃗¹_𝑚𝑛cos𝑝𝜋/2+𝑘_𝑦𝑚𝑛𝑦𝜔,𝐶⃗¹_𝑚𝑛=𝑐₁₄𝑘_𝑦𝑚𝑛𝑢_𝑧𝑚𝑛+𝑘_𝑧𝑢_𝑦𝑚𝑛+𝑐₁₆𝑘_𝑥𝑛𝑢_𝑦𝑚𝑛+𝑘_𝑦𝑚𝑛𝑢_𝑥𝑚𝑛𝑐₄₆𝑘_𝑦𝑚𝑛𝑢_𝑧𝑚𝑛+𝑘_𝑧𝑢_𝑦𝑚𝑛+𝑐₆₆𝑘_𝑥𝑛𝑢_𝑦𝑚𝑛+𝑘_𝑦𝑚𝑛𝑢_𝑥𝑚𝑛𝑐₄₅𝑘_𝑦𝑚𝑛𝑢_𝑧𝑚𝑛+𝑘_𝑧𝑢_𝑦𝑚𝑛+𝑐₅₆𝑘_𝑥𝑛𝑢_𝑦𝑚𝑛+𝑘_𝑦𝑚𝑛𝑢_𝑥𝑚𝑛𝐶⃗_2𝑚𝑛,=𝑐₁₁𝑘_𝑥𝑛𝑢_𝑥𝑚𝑛+𝑐₁₂𝑘_𝑦𝑚𝑛𝑢_𝑦𝑚𝑛+𝑐₁₃𝑘_𝑧𝑢_𝑧𝑚𝑛+𝑐₁₅𝑘_𝑥𝑛𝑢_𝑧𝑚𝑛+𝑘_𝑧𝑢_𝑥𝑚𝑛𝑐₁₆𝑘_𝑥𝑛𝑢_𝑥𝑚𝑛+𝑐₂₆𝑘_𝑦𝑚𝑛𝑢_𝑦𝑚𝑛+𝑐₃₆𝑘_𝑧𝑢_𝑧𝑚𝑛+𝑐₅₆𝑘_𝑥𝑛𝑢_𝑧𝑚𝑛+𝑘_𝑧𝑢_𝑥𝑚𝑛𝑐₁₅𝑘_𝑥𝑛𝑢_𝑥𝑚𝑛+𝑐₂₅𝑘_𝑦𝑚𝑛𝑢_𝑦𝑚𝑛+𝑐₃₅𝑘_𝑧𝑢_𝑧𝑚𝑛+𝑐₅₅𝑘_𝑥𝑛𝑢_𝑧𝑚𝑛+𝑘_𝑧𝑢_𝑥𝑚𝑛. (A.12)(A.13)(A.14)(A.15)(A.16)