ARRAW: anti-resonant reflecting acoustic waveguides

We first consider the case of pure shear waves, with solely out-of-plane displacement fields (${u}_{x}^{\left(j\right)}={u}_{z}^{\left(j\right)}=0$). These are uncoupled from longitudinal modes throughout the structure (see the derivations in appendix A or [24]). The guidance condition (vanishing amplitude of incoming waves in the outermost layer) yields a rather complex transcendental equation relating Ω, β and ${k}_{\mathrm{s}/\mathrm{p}}^{\left(j\right)}$ (not shown here). However, it can be simplified by considering separately the modes which are symmetric and anti-symmetric with respect to the $\hat{\mathbf{y}}$-$\hat{\mathbf{z}}$ symmetry plane of the structure at x = 0. The resulting transcendental equations for the symmetric and anti-symmetric modes can be expressed as

$$\begin{equation}\left(-1+{\text{e}}^{\text{i}c{k}_{\mathrm{s}}^{\left(1\right)}}{r}_{\mathrm{s}}\right)-{\text{e}}^{2\text{i}d{k}_{\mathrm{s}}^{\left(2\right)}}{r}_{\mathrm{s}}\left(-{r}_{\mathrm{s}}+{\text{e}}^{\text{i}c{k}_{\mathrm{s}}^{\left(1\right)}}\right)=0,\end{equation} \tag{ 3 }$$

$$\begin{equation}\left(1+{\text{e}}^{\text{i}c{k}_{\mathrm{s}}^{\left(1\right)}}{r}_{\mathrm{s}}\right)-{\text{e}}^{2\text{i}d{k}_{\mathrm{s}}^{\left(2\right)}}{r}_{\mathrm{s}}\left({r}_{\mathrm{s}}+{\text{e}}^{\text{i}c{k}_{\mathrm{s}}^{\left(1\right)}}\right)=0,\end{equation} \tag{ 4 }$$

respectively, where as marked in figure 1, c and d are the thicknesses of the core and cladding layers, and

$$\begin{equation}{r}_{\mathrm{s}}=\frac{{k}_{\mathrm{s}}^{\left(1\right)}{\mu }^{\left(1\right)}-{k}_{\mathrm{s}}^{\left(2\right)}{\mu }^{\left(2\right)}}{{k}_{\mathrm{s}}^{\left(1\right)}{\mu }^{\left(1\right)}+{k}_{\mathrm{s}}^{\left(2\right)}{\mu }^{\left(2\right)}},\end{equation} \tag{ 5 }$$

is the acoustic reflection coefficient for a pure S wave propagating in medium '1', reflecting off the interface with medium '2'.

2.1.1. Optics-like anti-resonance condition

We note that equations (3) and (4) also describe the symmetric and antisymmetric s-polarized optical waveguiding modes, if we replace r_s with the optical reflection coefficient for non-magnetic materials r^opt = (k⁽¹⁾ − k⁽²⁾)/(k⁽¹⁾ + k⁽²⁾). This is thanks to the purely transverse nature of these modes, and the close mapping between acoustic and optical boundary conditions (with the continuity of tangential components of the electric field, and of normal components of the magnetic field, mirroring the continuity of acoustic displacement and normal stress components, respectively).

Furthermore, to further explore the analogy to optics for ARRAWs, we can simplify the conditions given in equations (3) and (4) by eliminating the dependence on the thickness of the core layer (c) from the transcendental equations. To this end, we choose the core radius c to correspond to the lowest order symmetric modes supported by the core: [14]

$$\begin{equation}{k}_{\mathrm{s}}^{\left(1\right)}c=\left(2n+1\right)\pi ,\end{equation} \tag{ 6 }$$

for n = 0, 1, 2, ..., and arrive at the simplified symmetric ARRAW condition from (3)

$$\begin{equation}{\text{e}}^{2\text{i}{k}_{\mathrm{s}}^{\left(2\right)}d}={r}_{\mathrm{s}}.\end{equation} \tag{ 7 }$$

Similarly, the lowest order antisymmetric modes (4) are found for a core width satisfying

$$\begin{equation}{k}_{\mathrm{s}}^{\left(1\right)}c=2n\pi ,\end{equation} \tag{ 8 }$$

which simplifies (4) to the ARRAW condition given in (7), identical to the symmetric ARRAW modes.

To simplify this condition even further, we can consider waves propagating in the core at a glancing incidence to the cladding interface, where the reflection coefficient r_s ≈ −1. We then retrieve from (7) the approximate relation:

$$\begin{equation}{k}_{\mathrm{s}}^{\left(2\right)}d=\left(2\;m+1\right)\frac{\pi }{2},\end{equation} \tag{ 9 }$$

for m = 0, 1, 2, ..., which can be used to identify anti-resonant behavior in the dispersion of acoustic waveguides. This condition is analogous to that found for ARROWs, which reads k⁽²⁾d = (2m + 1)π/2.

2.1.2. Exact dispersion relation

The complete, exact dispersion relation of the symmetric, pure shear modes of the planar waveguide with arbitrary core radius c is found by solving (3) numerically, and includes families of both conventionally guided and leaky modes. To differentiate between the two, and to establish a parallel between acoustic and optical frameworks, we arbitrarily choose the largest shear velocity of the system ${v}_{\mathrm{s}}^{\left(1\right)}$ as a reference, and define an acoustic effective mode index

$$\begin{equation}{n}_{\text{eff}}=\frac{\text{Re}\left(\beta \right)}{{k}_{\mathrm{s},0}^{\left(1\right)}},\end{equation} \tag{ 10 }$$

where ${k}_{\mathrm{s},0}^{\left(1\right)}={\Omega}/{v}_{\mathrm{s}}^{\left(1\right)}$. Conventionally guided modes, characterized by a real longitudinal wavenumber β, and the acoustic field propagating predominantly in the slow cladding layer thus correspond to $1{< }{n}_{\text{eff}}{< }{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{s}}^{\left(2\right)}$. Conversely, leaky modes found for n_eff < 1 have complex longitudinal wavenumber.

These features are demonstrated in figure 2, where we analyze the symmetric, pure shear modes of a Si(core)/SiO₂(inner cladding)/Si(outer cladding) planar waveguide. For this setup, the ratio of shear velocities ${v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{s}}^{\left(2\right)}\approx 1.42$. The core width is set to c = 1 μm, and we vary the cladding width d from approximately 0 to d = 0.5 μm—a range of sizes which, for acoustic modes propagating at frequency Ω/2π = 15 GHz, spans ${k}_{\mathrm{s},0}^{\left(2\right)}d\approx 0$ to 4π. In figure 2(a) we plot the effective mode index n_eff, and in (b) the normalized attenuation parameter $\text{Im}\left(\beta \right)/{k}_{\mathrm{s},0}^{\left(1\right)}$.

Zoom In Zoom Out Reset image size

In figure 2(a) we clearly identify the families of modes conventionally guiding the acoustic waves in the inner cladding layer, characterized by $1{< }{n}_{\text{eff}}{< }{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{s}}^{\left(2\right)}$. These modes exhibit purely real transverse wavenumbers in the cladding (${k}_{\mathrm{s}}^{\left(2\right)}$), and purely imaginary transverse wavenumbers in the core (${k}_{\mathrm{s}}^{\left(1\right)}$). This indicates oscillatory behavior of the fields in the cladding, along the transverse direction $\hat{\mathbf{x}}$, and exponential localization of the fields to the interfaces in the core. We clearly identify these features in the two upper panels of figure 2(c), which show the displacement field profile of two modes depicted as A and B in figure 2(a) and (b), differentiated by the number of nodes of the displacement field in the cladding (with ${k}_{\mathrm{s}}^{\left(2\right)}d\approx 3\pi /2$ in A and 5π/2 in B).

Modes below the shear core sound line (n_eff < 1) are characterized by complex wavenumbers β, and ${k}_{\mathrm{s}}^{\left(j\right)}$'s, meaning that the field in the core does not simply exponentially decay with x away from the core/cladding interface, but exhibits oscillations governed by $\text{Re}\left({k}_{\mathrm{s}}^{\left(1\right)}\right)$. Simultaneously, in a manner consistent with optical leaky modes, the imaginary part of ${k}_{\mathrm{s}}^{\left(1\right)}$ becomes negative, and the outgoing fields increase exponentially with x outside the structure. Just below the n_eff = 1 line, the glancing modes form flat-dispersion sections where acoustic waves are confined predominantly to the core (see panels C and E in figure 2(c); modes C, D and E were selected to match the local minima of normalized loss $\text{Im}\left(\beta \right)/{k}_{\mathrm{s}}^{\left(1\right)}$. This behavior is described approximately by the horizontal gray dotted lines depicting resonances of the core (n = 0 in (6)), which cross with the colored dashed lines indicating the cladding anti-resonance (m = 0, 1, 2 in (9)) near the minimum of the loss function. This agreement breaks down for mode D, which lies far below the n_eff = 1 line, and does not meet the glancing incidence criterion. We also identify a prominent anti-crossing behavior near the crossings between the resonances of the core (horizontal dotted gray lines) and resonances of the cladding (${k}_{\mathrm{s}}^{\left(2\right)}d=m\pi$, anti-crossings marked with hollow circles). At these points the cladding transmission reaches a local maximum, suppressing the formation of a waveguiding mode.

This simple analysis allows us to identify the ARRAW modes as the leaky acoustic modes localized to the fast medium, and found at the local minima of loss, right below the sound-line corresponding to the velocity of waves in the fastest medium (core)—here illustrated as modes C and E. As we show below, this definition naturally extends to polarizations and media supporting the propagation of P waves.

In-plane polarization of the waveguiding modes (${u}_{y}^{\left(i\right)}=0$) necessarily couples the in-plane shear (S) and longitudinal (P) waves. Consequently, the transcendental equation for the modes is more complex (see derivation in appendix A), even if we focus on a particular symmetry of both the components. Nevertheless, it submits to numerical solution. As in the out-of-plane polarization case, we focus on symmetric modes for simplicity. Furthermore, to observe anti-resonant behavior of longitudinal waves, characterized by substantially longer wavelengths (since ${v}_{\mathrm{p}}^{\left(2\right)}/{v}_{\mathrm{s}}^{\left(2\right)}\approx 1.6$), we consider claddings of larger thickness, comparable to the longitudinal wavelength in the cladding medium.

In the following discussion, it is useful to draw from the analytical framework presented in appendix A, and treat S (shear) and P (longitudinal) components of the acoustic wave (denoted as u_s and u_p) as if they were independent, and characterize regimes in which these components behave as evanescent (exponentially decaying away from interfaces in either direction), conventionally guided (in the inner cladding) or leaky waves. While, as we show below, the coupling between the S and P waves necessarily blurs the characteristics of these regimes, this approach is instructive in developing intuition about the in-plane polarization acoustic guidance.

The dispersion relations and loss of the symmetric, in-plane modes are shown in figures 3(a) and (b), respectively. The field distributions (separated into S and P components described as the norms of u_s and u_p fields defined in appendix A) of modes marked as A–F are shown in (c), chosen to best represent 4 distinct regimes of acoustic guidance. The first, denoted in figure 3(a) as conventional S and evanescent P ($1{< }{n}_{\text{eff}}{< }{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{s}}^{\left(2\right)}$), is characterized by real β, ${k}_{\mathrm{s}}^{\left(2\right)}$, and imaginary ${k}_{\mathrm{s}}^{\left(1\right)}$ (indicating conventional S guidance in the cladding), as well as purely imaginary transverse P wave ${k}_{\mathrm{p}}^{\left(i\right)}$'s (indicating surface states with fields exponentially decaying away from the interface in the core and the outer layer). Two examples of such modes are shown in panels A and B in figure 3(c).

Zoom In Zoom Out Reset image size

For n_eff < 1, β becomes complex, and the shear components of the fields diverge exponentially outside of the structure as in the out-of-plane case (see panels C–F), as expected for the leaky S modes. In particular, as n_eff is reduced below ${v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(2\right)}\approx 0.89$, we first find the regime (${n}_{\text{eff}}{ >}{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(1\right)}\approx 0.6$) in which P components become conventionally guided in the inner cladding (see modes D and E). Modes found in this regime exhibit very particular characteristics, as they mix the conventional-like localization of the P waves with the leaky-like exponential increase of displacement fields in the outer cladding layer (see blue lines in D and E) due to the negative imaginary component of ${k}_{\mathrm{p}}^{\left(1\right)}$ (resulting from complex β, or equivalently, coupling to the S waves). Furthermore, for particular geometric parameters the normalized loss of these modes (figure 3(b)) appears to dip toward 0—a behavior which we identify with the onset of the simultaneous anti-resonant guidance of S waves (note the localization of S waves represented by orange lines to the core in D and E) and conventional P guidance.

Finally, further reducing n_eff below ${v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(1\right)}$ brings us to the leaky S and P regime, in which the displacement fields of both the P and S contributions oscillate in the core and in the inner cladding layer. We can therefore expect to find modes for which both the P and S waves build up an anti-resonance in the inner cladding, by approximately meeting the ARRAW condition given for S waves in (9). These conditions will not be fulfilled exactly, since the transverse wavenumbers ${k}_{\mathrm{s}}^{\left(2\right)}$ and ${k}_{\mathrm{p}}^{\left(2\right)}$ are imaginary, and the two components are coupled. Nevertheless, we can find ARRAW modes, such as the one shown in panel F, which meet our previous definition: they are characterized by a local minimum of loss, lie immediately below the ${n}_{\text{eff}}={v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(1\right)}$ sound line, and exhibit a strong localization of both the P and S components in the core.

We should also note that the combination of Si/SiO₂ materials forbids the formation of modes with both S and P components which are conventionally guided, since the regions of effective mode velocities $\left({v}_{\mathrm{s}}^{\left(2\right)},{v}_{\mathrm{s}}^{\left(1\right)}\right)$ and $\left({v}_{\mathrm{p}}^{\left(2\right)},{v}_{\mathrm{p}}^{\left(1\right)}\right)$ do not overlap. Such regions could be found for other combinations of materials, e.g. Si/As₂S₃ or SiO₂/As₂S₃.

In figure 4(a) we present the dispersion of torsional modes, and identify two families of modes: conventional S, with $1{< }{n}_{\text{eff}}{< }{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{s}}^{\left(2\right)}$, and leaky S, with n_eff < 1. The azimuthal (and only) component of the displacement field is shown, for a collection of modes, in figure 4(c). For A and B, the conventionally guided S modes are localized to the cladding layer, and decay exponentially in the outermost layer due to the purely imaginary transverse wavenumber ${k}_{\mathrm{s}}^{\left(3\right)}=\text{i}{\kappa }_{\mathrm{s}}^{\left(3\right)}$ (with ${\kappa }_{\mathrm{s}}^{\left(3\right)}{ >}0$). For the leaky modes C–E, we find that the effective mode index n_eff of all the modes increases with the cladding thickness d, until the dispersion crosses into the conventional S guidance regime. In particular, by tracing the evolution of the branch with modes C and D, we see that as the dispersion approaches the ${v}_{\mathrm{s}}^{\left(1\right)}$ sound line, oscillations in the cladding layer become more pronounced, and the anti-resonant response becomes significantly stronger, leading to substantial quenching of losses, as observed previously in optical ARROW systems.

Zoom In Zoom Out Reset image size

Dilatational modes share the fundamental dispersion characteristics of in-plane modes in planar waveguides, with the four regimes marked in figure 5(a) and illustrated in figure 5(c), mixing evanescent, conventional, and leaky characteristics of the S (orange lines) and P (blue lines) waves of the displacement field. Parameters of the waveguides considered in this calculation are identical to those discussed in section 3.1.

Zoom In Zoom Out Reset image size

In particular, for ${v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(1\right)}{< }{n}_{\text{eff}}{< }{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(2\right)}$ we again find the peculiar modes (C and D) characterized by very small loss, for which the P waves have the dual characteristic of conventional-like localization to the inner cladding, and leaky-like exponential growth in the outer cladding. For these modes, S waves are clearly localized through the anti-resonance to the core of the waveguide.

Furthermore, as for the in-plane modes of the planar waveguide, we expect to find the ARRAW modes right below the fastest sound velocity ${v}_{\mathrm{p}}^{\left(1\right)}$ line, in the leaky S and P regime. While these modes are not clearly defined by a dip in loss, we can nevertheless identify ARRAW behavior in the field distributions of the mode. For example, mode E carries the characteristics of both the S and P components localized to the core, and oscillatory behavior of the components in the cladding layer.

The above-identified ARRAW modes would clearly not constitute very good acoustic waveguiding channels, as their normalized loss barely reaches 10⁻² (see point F in figure 3 and point E in figure 5). This is mostly due to the fact that the core layers of both the planar and cylindrical systems investigated to this point are too narrow to fit multiple wavelengths of the longer, P acoustic waves (${k}_{\mathrm{p}}^{\left(1\right)}a,{k}_{\mathrm{p}}^{\left(1\right)}c\ll \pi$). While their geometric parameters were chosen to provide insights into the different regimes of operation of layered acoustic waveguides, we can now consider systems with larger core thicknesses that will provide much better examples of ARRAW guidance, characterized by lower losses and stronger localization of the field to the core, for application in nonlinear Brillouin scattering.

The mechanical quality factor Q_m is optimized by tuning the geometric parameters of the structure (core radius a and cladding thickness d), and mechanical frequency Ω. Since the underlying physics of acoustic systems is linear, the modes are invariant under simultaneous rescaling of the geometric parameters and wavelength (or inverse frequency), and consequently, we can consider two of these three parameters as independent. Here, we fix Ω/2π = 15 GHz, and optimize a and d. For clarity, in figure 6(a) we show the dispersion relation and mechanical quality factors of the highest-Q_m ARRAW modes lying right below the last sound line (${n}_{\text{eff}}{\leqslant}{v}_{\mathrm{s}}^{\left(1\right)}/{v}_{\mathrm{p}}^{\left(1\right)}$), as function of the cladding thickness d. For each of the ARRAW modes denoted with solid blue lines, we identify a resonant dependence of the Q_m on d. These features clearly point to the resonances in the reflection within the cladding layer, and the ARRAW nature of the guidance.

Zoom In Zoom Out Reset image size

Since we aim to optimize the waveguiding properties of the structure, we consider larger core radii a = 2 μm, which offer mechanical quality factors over 1000 (structures with different core radii are analyzed in table 1). As the maximum Q_m increases with the core radius (see results in table 1), it is tempting to favor the larger core geometries. However, we should note that our calculations do not account for the viscous losses in the material, which typically limit the Q_m to the order of ∼1000, and thus suppress the advantage of eliminating the acoustic radiative dissipation channels. Furthermore, it was shown in previous reports that the PE contribution to SBS gain exhibits an approximate a⁻² scaling with the transverse dimension of the waveguide [21].

The particular choice of the optical mode which can undergo BSBS coupling via a selected ARRAW mode is predominantly limited by the phase-matching condition (11), which dictates the magnitude of the longitudinal wavenumber k₁ of the optical mode. Therefore, the only parameters which we can optimize here are the type of optical mode (TE, TM or hybrid), and its effective mode index ${n}_{\text{eff}}^{\text{opt}}$. We first focus on the lowest ${n}_{\text{eff}}^{\text{opt}}$ (and high-order) modes supported by the core, and later (figure 6(d)) explore the dependence on the order of the optical mode. In particular, in figure 6(b) we analyze the BSBS gain coupling between these three types of modes (TM₄, TE₄ and EH₁₇/HE₁₇) as a function of cladding thickness—a parameter which should have little effect on the optical guidance, but governs the ARRAW behavior. We find that the gain closely follows the dependence of the quality factor Q_m, suggesting that neither the overlap integral ${\mathcal{Q}}_{1}^{\left(\text{PE}\right)}$ nor the normalization factors ${\mathcal{E}}_{b}$ and ${\mathcal{P}}^{\left(i\right)}$ change significantly with d.

The orange dot in the plots in figure 6(a) and (b)denotes the parameters of optical and acoustic modes analyzed in detail in figure 6(c). The radial profile of the acoustic field, shown in the top panel of figure 6(a), indicates significant localization to the core—as expected for large-Q_m ARRAW modes. Together with the core-localized optical mode (middle panel), this results in a strong localization of the ${\tilde {\mathcal{Q}}}_{1}^{\left(\text{PE}\right)}$(r) overlap integral kernel shown in the bottom panel.

To further enhance the Brillouin gain, we can consider changing the geometric parameters of the waveguide (e.g. core radius a), operating mechanical frequency Ω, or explore coupling to different orders of the optical modes. We provide a comparison of selected BSBS ARRAW systems in table 1, and find a few guiding principles for designing such systems.

• As reported by Rakich et al [21], the smaller cross section waveguides yield larger BSBS gain—however, this principle trades off against the sharp decrease in mechanical quality factor for small core radii; simultaneous increase of mechanical frequencies (e.g. toward 30 GHz frequency) should allow us to retain high Q_m's due to the linear nature of the acoustic physics, but the overall gain would likely become suppressed by the increased non-radiative acoustic losses at higher frequencies,
• Dependence on the order of optical mode (or effective optical mode index) is not monotonic, and in fact is the smallest for the most homogeneous field of the lowest order mode; this dependence is shown, for a number of hybrid modes from HE₁₇ to HE₁₂, in figure 6(d).

It is instructive to compare the nonlinear performance of the investigated ARRAW to other BSBS waveguiding systems. The mechanical quality factor Q_m of our structure can easily reach 10³, which is a result comparable to that found in suspended waveguides, in which the radiative dissipation of acoustic waves is suppressed, and Q_m becomes limited by the intrinsic viscous losses in silicon [8, 28]. Furthermore, the maximum Brillouin gain found in our system (nearly 1000 (W m)⁻¹) is of a similar order to that reported in GHz BSBS systems, including those based on sub-micron silicon slot waveguides where the optoacoustic interaction is dominated by radiation pressure [29], or relying on materials with different optoacoustic properties, such as chalcogenides [12].

Our simple designs can be further modified to better suit integrated platforms by considering finite outer-cladding layers or multiple anti-resonant layers. The principles of ARRAW guidance can also be combined with other mechanisms, for example by using the anti-resonant reflection to suppress the acoustic dissipation from exposed cores of rib waveguides into the substrate. Alternatively, the entire designs could be reversed and implement optical anti-resonant guidance and conventional acoustic TIR.

We start with the simpler case of out-of-plane polarization, in which modes are composed of shear waves only, and displacement fields have u_y components only, meaning that the boundary conditions demand the continuity of u_y components of the displacement field and continuity of the T_xy element of the stress tensor.

We use the ansatz

$$\begin{equation}{\mathbf{u}}^{\left(j\right)}\left(x,z,t\right)=\hat{\mathbf{y}}\left({u}_{+}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}+{u}_{-}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\right){\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}+\text{c.c.},\end{equation} \tag{ A.1 }$$

where ${u}_{+}^{\left(j\right)}$ and ${u}_{-}^{\left(j\right)}$ are undetermined coefficients for the right- and left-propagating fields in the jth layer respectively. We can thus represent the boundary conditions at x = x_j between the jth and (j + 1)th layer as

$$\begin{equation}\underset{{M}^{\left(j\right)}\left({x}_{j}\right)}{\underbrace{\left[\begin{bmatrix}{ccc}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}{x}_{j}}& \hfill \hfill & {\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}{x}_{j}}\\ \hfill {k}_{\mathrm{s}}^{\left(j\right)}{\mu }^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}{x}_{j}}\hfill & \hfill \hfill & \hfill -{k}_{\mathrm{s}}^{\left(j\right)}{\mu }^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}{x}_{j}}\hfill \end{bmatrix}\right]}}\underset{{\tilde {\mathbf{u}}}^{\left(j\right)}}{\underbrace{\left[\begin{bmatrix}{c}\hfill {u}_{+}^{\left(j\right)}\hfill \\ \hfill {u}_{-}^{\left(j\right)}\hfill \end{bmatrix}\right]}}=\underset{{M}^{\left(j+1\right)}\left({x}_{j}\right)}{\underbrace{\left[\begin{bmatrix}{ccc}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j+1\right)}{x}_{j}}& \hfill \hfill & {\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j+1\right)}{x}_{j}}\\ \hfill {k}_{\mathrm{s}}^{\left(j+1\right)}{\mu }^{\left(j+1\right)}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j+1\right)}{x}_{j}}\hfill & \hfill \hfill & \hfill -{k}_{\mathrm{s}}^{\left(j+1\right)}{\mu }^{\left(j+1\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j+1\right)}{x}_{j}}\hfill \end{bmatrix}\right]}}\underset{{\tilde {\mathbf{u}}}^{\left(j+1\right)}}{\underbrace{\left[\begin{bmatrix}{c}\hfill {u}_{+}^{\left(j+1\right)}\hfill \\ \hfill {u}_{-}^{\left(j+1\right)}\hfill \end{bmatrix}\right]}}.\end{equation} \tag{ A.2 }$$

We now multiply both sides by ${\left[{M}^{\left(j+1\right)}\left({x}_{j}\right)\right]}^{-1}$ to obtain

$$\begin{equation}{\left[{M}^{\left(j+1\right)}\left({x}_{j}\right)\right]}^{-1}{M}^{\left(j\right)}\left({x}_{j}\right){\tilde {\mathbf{u}}}^{\left(j\right)}\doteq {S}^{\left(j\right)}{\tilde {\mathbf{u}}}^{\left(j\right)}={\tilde {\mathbf{u}}}^{\left(j+1\right)}.\end{equation} \tag{ A.3 }$$

This formalism relates the fields in the last medium j = 5 with the first medium j = 1 (${S}^{\left(\text{all}\right)}{\tilde {\mathbf{u}}}^{\left(1\right)}={\tilde {\mathbf{u}}}^{\left(5\right)}$), through the matrix

$$\begin{equation}{S}^{\left(\text{all}\right)}={S}^{\left(4\right)}{S}^{\left(3\right)}{S}^{\left(2\right)}{S}^{\left(1\right)}.\end{equation} \tag{ A.4 }$$

To find waveguiding modes, we require that fields in the first and the last mediums should be outgoing only, i.e.

$$\begin{equation}{\tilde {\mathbf{u}}}^{\left(1\right)}=\left[\begin{bmatrix}{c}0\\ \hfill {u}_{-}^{\left(1\right)}\hfill \end{bmatrix}\right],\quad {\tilde {\mathbf{u}}}^{\left(5\right)}=\left[\begin{bmatrix}{c}\hfill {u}_{+}^{\left(5\right)}\hfill \\ 0\end{bmatrix}\right].\end{equation} \tag{ A.5 }$$

This condition only holds when ${S}_{22}^{\left(\text{all}\right)}=0$, which defines a 4th order polynomial in the quantities ${k}_{\mathrm{s}}^{\left(i\right)}$. We can simplify it by considering the symmetric and anti-symmetric modes of the structure separately. Expressing the field u⁽³⁾ through the transfer matrix, imposing the above symmetry conditions, and accounting for the vanishing of the incident wave in medium 1, we arrive at the resonance condition:

$$\begin{equation}{\left({W}^{{\pm}}{S}^{\left(2\right)}{S}^{\left(1\right)}\right)}_{22}=0,\end{equation} \tag{ A.6 }$$

where W^± imposes the symmetry/antisymmetry of u⁽³⁾ at x = 0, and is defined by

$$\begin{equation}{W}^{{\pm}}=\left[\begin{bmatrix}{ccc}\hfill 1\hfill & \hfill \hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill \hfill & \hfill {\pm}1\hfill \end{bmatrix}\right].\end{equation} \tag{ A.7 }$$

The condition (A.6) is equivalent to equations (3) and (4).

For in-plane polarization, the displacement fields lie in the $\hat{\mathbf{x}}$-$\hat{\mathbf{y}}$ plane, and they have both longitudinal and transverse contributions. Importantly, these two contributions exhibit different transverse wavevectors. The displacement field components thus have two right- and left-propagating field coefficients each (${u}_{\mathrm{p},{\pm}}^{\left(j\right)}$ and ${u}_{\mathrm{s},{\pm}}^{\left(j\right)}$):

$$\begin{equation}{\mathbf{u}}_{\mathrm{p}}^{\left(j\right)}\left(x,z,t\right)=\left[\left(\hat{\mathbf{x}}{k}_{\mathrm{p}}^{\left(j\right)}+\hat{\mathbf{z}}\beta \right){u}_{\mathrm{p}+}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x}\right.+\;\left.\left(-\hat{\mathbf{x}}{k}_{\mathrm{p}}^{\left(j\right)}+\hat{\mathbf{z}}\beta \right){u}_{\mathrm{p}-}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x}\right]{\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}+\text{c.c.},\end{equation} \tag{ A.8 }$$

$$\begin{equation}{\mathbf{u}}_{\mathrm{s}}^{\left(j\right)}\left(x,z,t\right)=\left[\left(-\hat{\mathbf{x}}\beta +\hat{\mathbf{z}}{k}_{\mathrm{s}}^{\left(j\right)}\right){u}_{\mathrm{s}+}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\right.+\;\left.\left(\hat{\mathbf{x}}\beta +\hat{\mathbf{z}}{k}_{\mathrm{s}}^{\left(j\right)}\right){u}_{\mathrm{s}-}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\right]{\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}+\text{c.c.},\end{equation} \tag{ A.9 }$$

$$\begin{equation}{\mathbf{u}}^{\left(j\right)}\left(x,z,t\right)={\mathbf{u}}_{\mathrm{p}}^{\left(j\right)}\left(x,z,t\right)+{\mathbf{u}}_{\mathrm{s}}^{\left(j\right)}\left(x,z,t\right).\end{equation} \tag{ A.10 }$$

The relevant elements of the stress tensor which need to be continuous at each interface are ${T}_{xx}^{\left(j\right)}$, and ${T}_{xz}^{\left(j\right)}$.

The matrices M^(j)(x) and vectors ${\tilde {\mathbf{u}}}^{\left(j\right)}$ are thus given by:

$$\begin{equation}{\tilde {\mathbf{u}}}^{\left(j\right)}={\left[\begin{bmatrix}{cccc}\hfill {u}_{\mathrm{p}+}^{\left(j\right)}\hfill & \hfill {u}_{\mathrm{p}-}^{\left(j\right)}\hfill & \hfill {u}_{\mathrm{s}+}^{\left(j\right)}\hfill & \hfill {u}_{\mathrm{s}-}^{\left(j\right)}\hfill \end{bmatrix}\right]}^{\mathrm{T}},\end{equation} \tag{ A.11 }$$

$$\begin{equation}{M}^{\left(j\right)}\left(x\right)=\left[\begin{bmatrix}{ccccccc}\hfill {k}_{\mathrm{p}}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x}\hfill & \hfill \hfill & \hfill -{k}_{\mathrm{p}}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x}\hfill & \hfill \hfill & \hfill -\beta {\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\hfill & \hfill \hfill & \hfill \beta {\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\hfill \\ \hfill \beta {\text{e}}^{\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x}\hfill & \hfill \hfill & \hfill \beta {\text{e}}^{-\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x}\hfill & \hfill \hfill & \hfill {k}_{\mathrm{s}}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\hfill & \hfill \hfill & \hfill {k}_{\mathrm{s}}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}\hfill \\ \hfill {m}_{31}\hfill & \hfill \hfill & \hfill {m}_{32}\hfill & \hfill \hfill & \hfill {m}_{33}\hfill & \hfill \hfill & \hfill {m}_{34}\hfill \\ \hfill {m}_{41}\hfill & \hfill \hfill & \hfill {m}_{42}\hfill & \hfill \hfill & \hfill {m}_{43}\hfill & \hfill \hfill & \hfill {m}_{44}\hfill \end{bmatrix}\right],\end{equation} \tag{ A.12 }$$

with

$$\begin{equation}{m}_{31}=\left[\left({\lambda }^{\left(j\right)}+2{\mu }^{\left(j\right)}\right){\left({k}_{\mathrm{p}}^{\left(j\right)}\right)}^{2}+{\lambda }^{\left(j\right)}{\beta }^{2}\right]{\text{e}}^{\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x},\end{equation} \tag{ A.13 }$$

$$\begin{equation}{m}_{32}=\left[\left({\lambda }^{\left(j\right)}+2{\mu }^{\left(j\right)}\right){\left({k}_{\mathrm{p}}^{\left(j\right)}\right)}^{2}+{\lambda }^{\left(j\right)}{\beta }^{2}\right]{\text{e}}^{-\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x},\end{equation} \tag{ A.14 }$$

$$\begin{equation}{m}_{33}=-2{\mu }^{\left(j\right)}\beta {k}_{\mathrm{s}}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x},\end{equation} \tag{ A.15 }$$

$$\begin{equation}{m}_{34}=-2{\mu }^{\left(j\right)}\beta {k}_{\mathrm{s}}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x},\end{equation} \tag{ A.16 }$$

$$\begin{equation}{m}_{41}=2{\mu }^{\left(j\right)}\beta {k}_{\mathrm{p}}^{\left(j\right)}{\text{e}}^{\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x},\end{equation} \tag{ A.17 }$$

$$\begin{equation}{m}_{42}=-2{\mu }^{\left(j\right)}\beta {k}_{\mathrm{p}}^{\left(j\right)}{\text{e}}^{-\text{i}{k}_{\mathrm{p}}^{\left(j\right)}x},\end{equation} \tag{ A.18 }$$

$$\begin{equation}{m}_{43}={\mu }^{\left(j\right)}\left[{\left({k}_{\mathrm{s}}^{\left(j\right)}\right)}^{2}-{\beta }^{2}\right]{\text{e}}^{\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x},\end{equation} \tag{ A.19 }$$

$$\begin{equation}{m}_{44}=-{\mu }^{\left(j\right)}\left[{\left({k}_{\mathrm{s}}^{\left(j\right)}\right)}^{2}-{\beta }^{2}\right]{\text{e}}^{-\text{i}{k}_{\mathrm{s}}^{\left(j\right)}x}.\end{equation} \tag{ A.20 }$$

As in the previous section, this formulation of matrices M^(j) allows us to construct matrices S^(j) as in (A.4), which will readily mix the S and P components of the fields in layers j and j + 1. This behavior signals coupling between the S and P waves through scattering at interfaces.

To find waveguiding modes, we again require that fields in the first and the last layers should be outgoing only, i.e.

$$\begin{equation}{\tilde {\mathbf{u}}}^{\left(5\right)}=\left[\begin{matrix}{c}{u}_{\mathrm{p}+}^{\left(5\right)}\\ 0\\ {u}_{\mathrm{s}+}^{\left(5\right)}\\ 0\end{matrix}\right],\quad {\tilde {\mathbf{u}}}^{\left(1\right)}=\left[\begin{matrix}{c}0\\ {u}_{\mathrm{p}-}^{\left(1\right)}\\ 0\\ {u}_{\mathrm{s}-}^{\left(1\right)}\end{matrix}\right].\end{equation} \tag{ A.21 }$$

Thus, we formulate the transfer problem as ${\tilde {\mathbf{u}}}^{\left(5\right)}={S}^{\left(\text{all}\right)}{\tilde {\mathbf{u}}}^{\left(1\right)}$, and look for the solutions to

$$\begin{equation}\mathrm{det}\left[\begin{bmatrix}{cc}\hfill {\left({S}^{\left(\text{all}\right)}\right)}_{22}\hfill & \hfill {\left({S}^{\left(\text{all}\right)}\right)}_{24}\hfill \\ \hfill {\left({S}^{\left(\text{all}\right)}\right)}_{42}\hfill & \hfill {\left({S}^{\left(\text{all}\right)}\right)}_{44}\hfill \end{bmatrix}\right]=0,\end{equation} \tag{ A.22 }$$

for discrete values of β and fixed Ω. As discussed in the main text, this transcendental equation will have both real and complex solutions, corresponding to conventional and ARRAW modes, respectively. Treating the real and imaginary parts of β as independent variables and applying one of the many standard multi-dimensional root-finding algorithms is computationally intensive and time-consuming. However, by treating β as a complex variable, we are able to apply methods of complex analysis; (A.22) is analytic everywhere except along a branch cut on the real axis. In order to calculate the dispersion relations presented in the main text, we applied the global complex roots and pole finding (GRPF) algorithm developed in [30].

The problem in question is that of a cylindrical rod of isotropic material and infinite length, as shown in figure 1(b). The problem is solved using cylindrical coordinates (r, θ, z) and assuming a time-harmonic solution of the form

$$\begin{equation}\mathbf{u}\left(r,\theta ,z,t\right)=\mathbf{U}\left(r,\theta \right){\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}+\text{c.c.}\end{equation} \tag{ B.1 }$$

The relevant elements of the stress tensor in cylindrical coordinates are identified as [24]

$$\begin{equation}{T}_{rr}=\lambda \left({\partial }_{r}{u}_{r}+\frac{1}{r}{\partial }_{\theta }{u}_{\theta }+\frac{{u}_{r}}{r}+{\partial }_{z}{u}_{z}\right)+2\mu {\partial }_{r}{u}_{r},\end{equation} \tag{ B.2 }$$

$$\begin{equation}{T}_{r\theta }=\mu \left(\frac{1}{r}{\partial }_{\theta }{u}_{r}-\frac{{u}_{\theta }}{r}+{\partial }_{r}{u}_{\theta }\right),\end{equation} \tag{ B.3 }$$

$$\begin{equation}{T}_{rz}=\mu \left({\partial }_{z}{u}_{r}+{\partial }_{r}{u}_{z}\right).\end{equation} \tag{ B.4 }$$

A convenient expansion for the modes of the cylinder can be derived by expressing the displacement fields u through two scalar potentials Φ and Ψ as

$$\begin{equation}\mathbf{u}=A\nabla {\Phi}+B\nabla {\times}\left(\hat{\mathbf{z}}{\Psi}\right)+C\nabla {\times}\nabla {\times}\left(\hat{\mathbf{z}}{\Psi}\right),\end{equation} \tag{ B.5 }$$

where A, B and C are free coefficients [24]. These potentials are solutions to the scalar wave equations

$$\begin{equation}\left[{\nabla }_{\perp }^{2}+\left(\frac{{\partial }^{2}}{\partial {z}^{2}}-\frac{1}{{v}_{\mathrm{p}}^{2}}\frac{{\partial }^{2}}{\partial {t}^{2}}\right)\right]{\Phi}\left(r,\theta ,z,t\right)=\left[{\nabla }_{\perp }^{2}+\left(\frac{{{\Omega}}^{2}}{{v}_{\mathrm{p}}^{2}}-{\beta }^{2}\right)\right]\phi \left(r,\theta \right){\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}=0,\end{equation} \tag{ B.6 }$$

$$\begin{equation}\left[{\nabla }_{\perp }^{2}+\left(\frac{{\partial }^{2}}{\partial {z}^{2}}-\frac{1}{{v}_{\mathrm{s}}^{2}}\frac{{\partial }^{2}}{\partial {t}^{2}}\right)\right]{\Psi}\left(r,\theta ,z,t\right)=\left[{\nabla }_{\perp }^{2}+\left(\frac{{{\Omega}}^{2}}{{v}_{\mathrm{s}}^{2}}-{\beta }^{2}\right)\right]\psi \left(r,\theta \right){\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}=0,\end{equation} \tag{ B.7 }$$

where we have made a similar ansatz for the scalar fields as for u (see (B.1)). The resulting equations for the transverse scalar fields ϕ and ψ can be solved by supposing an azimuthal dependence of e^imθ, which yields the radial dependence of the potentials as solutions to the Bessel equation χ_m(k_ir) (we will comment on the particular choice of χ), where ${k}_{i}^{2}={\left({\Omega}/{v}_{i}\right)}^{2}-{\beta }^{2}$ for i = p (ϕ) or s (ψ).

We thus arrive at the general form of the displacement field:

$$\begin{equation}{U}_{r}=A{k}_{\mathrm{p}}{\chi }_{m}^{\prime }\left({k}_{\mathrm{p}}r\right)+B\text{i}\beta {k}_{\mathrm{s}}{\chi }_{m}^{\prime }\left({k}_{\mathrm{s}}r\right)+C\text{i}m\frac{{\chi }_{m}\left({k}_{\mathrm{s}}r\right)}{r},\end{equation} \tag{ B.8 }$$

$$\begin{equation}{U}_{\theta }=A\text{i}m\frac{{\chi }_{m}\left({k}_{\mathrm{p}}r\right)}{r}-Bm\beta \frac{{\chi }_{m}\left({k}_{\mathrm{s}}r\right)}{r}-C{k}_{\mathrm{s}}{\chi }_{m}^{\prime }\left({k}_{\mathrm{s}}r\right),\end{equation} \tag{ B.9 }$$

$$\begin{equation}{U}_{z}=A\text{i}\beta {\chi }_{m}\left({k}_{\mathrm{p}}r\right)+B{k}_{\mathrm{s}}^{2}{\chi }_{m}\left({k}_{\mathrm{s}}r\right).\end{equation} \tag{ B.10 }$$

Therefore, we can observe that in the axially symmetric case, when m = 0, the displacement field decouples into torsional modes with only a u_θ component (A = B = 0), and dilatational modes (C = 0) with u_r and u_z components. We will focus on this case throughout.

Just as for the planar waveguide structure, we can decompose the displacement field in each layer into a basis of outgoing and incoming waves. This is determined by the choice of the functions χ from the family of Bessel and Hankel functions. To represent the outgoing wave, we choose χ in the form of Hankel functions of the first kind ${H}_{0}^{\left(1\right)}$; we will suppress the (1) superscript in all the following equations for clarity. To describe the incoming waves, we could use the Hankel functions of the second kind. However, since H⁽²⁾ diverge at r = 0, we instead consider the Bessel functions of the first kind J⁽¹⁾—while this component will not strickly describe incoming waves only, our expansion will remain complete. The displacement field is thus given by

$$\begin{equation}{\mathbf{u}}^{\left(j\right)}\left(r,z,t\right)=\hat{\boldsymbol{\theta }}\left({u}_{+}^{\left(j\right)}{H}_{0}^{\prime }\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)+{u}_{-}^{\left(j\right)}{J}_{0}^{\prime }\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\right){\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}.\end{equation} \tag{ B.11 }$$

Here the derivatives of special Bessel and Hankel functions (denoted by ') are calculated with respect to the entire argument of the respective function. The only relevant non-zero component of the stress tensor is T_rθ. Thus, the matrices M^(j)(r) and coefficient vectors ${\tilde {\mathbf{u}}}^{\left(j\right)}$ defined in appendix A are given by

$$\begin{equation}{\tilde {\mathbf{u}}}^{\left(j\right)}={\left[\begin{bmatrix}{cc}\hfill {u}_{+}^{\left(j\right)}\hfill & \hfill {u}_{-}^{\left(j\right)}\hfill \end{bmatrix}\right]}^{\mathrm{T}},\end{equation} \tag{ B.12 }$$

$$\begin{align}\hfill {M}^{\left(j\right)}\left(r\right)& =\left[\begin{bmatrix}{ccc}{H}_{0}^{\prime }\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)& \hfill \hfill & {J}_{0}^{\prime }\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\\ \hfill {\mu }^{\left(j\right)}{k}_{\mathrm{s}}^{\left(j\right)}\left({2H}_{0}^{{\prime\prime}}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)+{H}_{0}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\right)\hfill & \hfill \hfill & \hfill {\mu }^{\left(j\right)}{k}_{\mathrm{s}}^{\left(j\right)}\left({2J}_{0}^{{\prime\prime}}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)+{J}_{0}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\right)\hfill \end{bmatrix}\right]\hfill \\ \hfill & =\left[\begin{bmatrix}{ccc}\hfill -{H}_{1}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\hfill & \hfill \hfill & \hfill -{J}_{1}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\hfill \\ \hfill {\mu }^{\left(j\right)}{k}_{\mathrm{s}}^{\left(j\right)}{H}_{2}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\hfill & \hfill \hfill & \hfill {\mu }^{\left(j\right)}{k}_{\mathrm{s}}^{\left(j\right)}{J}_{2}\left({k}_{\mathrm{s}}^{\left(j\right)}r\right)\hfill \end{bmatrix}\right].\hfill \end{align} \tag{ B.13 }$$

The waveguiding condition is formulated in exactly the same manner as for the planar structure, by connecting the vectors of coefficients in the outermost ${\tilde {\mathbf{u}}}^{\left(3\right)}$ and innermost layer ${\tilde {\mathbf{u}}}^{\left(1\right)}$ through the M matrices (see A.3), and requiring that the innermost field is non-singular at r = 0 (no outgoing wave in layer (1); ${u}_{+}^{\left(1\right)}=0$) and the outermost field has no incoming wave components (since the Bessel function can be expressed as a sum of incoming and outgoing waves, we put ${u}_{-}^{\left(3\right)}=0$).

Following the same arguments as in the previous subsection, we again write down the displacement field as a sum of '+' and '−' components, with each one including contributions from both the P and S waves:

$$\begin{equation}{\mathbf{u}}_{\mathrm{p}}^{\left(j\right)}\left(r,z,t\right)=\left[\left(\hat{\mathbf{r}}{k}_{\mathrm{p}}^{\left(j\right)}{H}_{0\mathrm{p}}^{\prime }+\hat{\mathbf{z}}\text{i}\beta {H}_{0\mathrm{p}}\right){u}_{\mathrm{p}+}^{\left(j\right)}\right.\left.+\left(\hat{\mathbf{r}}{k}_{\mathrm{p}}^{\left(j\right)}{J}_{0\mathrm{p}}^{\prime }+\hat{\mathbf{z}}\text{i}\beta {J}_{0\mathrm{p}}\right){u}_{\mathrm{p}-}^{\left(j\right)}\right]{\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}+\text{c.c.},\end{equation} \tag{ B.14 }$$

$$\begin{equation}{\mathbf{u}}_{\mathrm{s}}^{\left(j\right)}\left(x,z,t\right)=\left[\left(\hat{\mathbf{r}}\text{i}\beta {H}_{0\mathrm{s}}^{\prime }+\hat{\mathbf{z}}{k}_{\mathrm{s}}^{\left(j\right)}{H}_{0\mathrm{s}}\right){u}_{\mathrm{s}+}^{\left(j\right)}\right.\left.+\left(\hat{\mathbf{r}}\text{i}\beta {J}_{0\mathrm{s}}^{\prime }+\hat{\mathbf{z}}{k}_{\mathrm{s}}^{\left(j\right)}{J}_{0\mathrm{s}}\right){u}_{\mathrm{s}-}^{\left(j\right)}\right]{\text{e}}^{\text{i}\left(\beta z-{\Omega}t\right)}+\text{c.c.},\end{equation} \tag{ B.15 }$$

$$\begin{equation}{\mathbf{u}}^{\left(j\right)}\left(x,z,t\right)={\mathbf{u}}_{\mathrm{p}}^{\left(j\right)}\left(x,z,t\right)+{\mathbf{u}}_{\mathrm{s}}^{\left(j\right)}\left(x,z,t\right),\end{equation} \tag{ B.16 }$$

where we have introduced the abbreviations ${H}_{0i}={H}_{0}\left({k}_{i}^{\left(j\right)}r\right)$ and ${J}_{0i}={J}_{0}\left({k}_{i}^{\left(j\right)}r\right)$, for i = s, p. The relevant components of the stress tensor are T_rr and T_rz. Therefore, we obtain the coefficients vectors

$$\begin{equation}{\tilde {\mathbf{u}}}^{\left(j\right)}={\left[\begin{bmatrix}{cccc}\hfill {u}_{\mathrm{p}+}^{\left(j\right)}\hfill & \hfill {u}_{\mathrm{p}-}^{\left(j\right)}\hfill & \hfill {u}_{\mathrm{s}+}^{\left(j\right)}\hfill & \hfill {u}_{\mathrm{s}-}^{\left(j\right)}\hfill \end{bmatrix}\right]}^{\mathrm{T}},\end{equation} \tag{ B.17 }$$

and matrices M for each layer j

$$\begin{equation}{M}^{\left(j\right)}\left(r\right)=\left[\begin{bmatrix}{ccccccc}\hfill {k}_{\mathrm{p}}^{\left(j\right)}{H}_{0\mathrm{p}}^{\prime }\hfill & \hfill \hfill & \hfill {k}_{\mathrm{p}}^{\left(j\right)}{J}_{0\mathrm{p}}^{\prime }\hfill & \hfill \hfill & \hfill \text{i}\beta {H}_{0\mathrm{s}}^{\prime }\hfill & \hfill \hfill & \hfill \text{i}\beta {J}_{0\mathrm{s}}^{\prime }\hfill \\ \hfill \text{i}\beta {H}_{0\mathrm{p}}\hfill & \hfill \hfill & \hfill \text{i}\beta {J}_{0\mathrm{p}}\hfill & \hfill \hfill & \hfill {k}_{\mathrm{s}}^{\left(j\right)}{H}_{0\mathrm{s}}\hfill & \hfill \hfill & \hfill {k}_{\mathrm{s}}^{\left(j\right)}{J}_{0\mathrm{s}}\hfill \\ \hfill {m}_{31}\hfill & \hfill \hfill & \hfill {m}_{32}\hfill & \hfill \hfill & \hfill {m}_{33}\hfill & \hfill \hfill & \hfill {m}_{34}\hfill \\ \hfill {m}_{41}\hfill & \hfill \hfill & \hfill {m}_{42}\hfill & \hfill \hfill & \hfill {m}_{43}\hfill & \hfill \hfill & \hfill {m}_{44}\hfill \end{bmatrix}\right],\end{equation} \tag{ B.18 }$$

where

$$\begin{equation}{m}_{31}=2{\mu }^{\left(j\right)}{\left({k}_{\mathrm{p}}^{\left(j\right)}\right)}^{2}{H}_{0\mathrm{p}}^{{\prime\prime}}-{\lambda }^{\left(j\right)}\left({\left({k}_{\mathrm{p}}^{\left(j\right)}\right)}^{2}+{\beta }^{2}\right){H}_{0\mathrm{p}},\end{equation} \tag{ B.19 }$$

$$\begin{equation}{m}_{32}=2{\mu }^{\left(j\right)}{\left({k}_{\mathrm{p}}^{\left(j\right)}\right)}^{2}{J}_{0\mathrm{p}}^{{\prime\prime}}-{\lambda }^{\left(j\right)}\left({\left({k}_{\mathrm{p}}^{\left(j\right)}\right)}^{2}+{\beta }^{2}\right){J}_{0\mathrm{p}},\end{equation} \tag{ B.20 }$$

$$\begin{equation}{m}_{33}=2\text{i}\beta {\mu }^{\left(j\right)}{\left({k}_{\mathrm{s}}^{\left(j\right)}\right)}^{2}{H}_{0\mathrm{s}}^{{\prime\prime}},\end{equation} \tag{ B.21 }$$

$$\begin{equation}{m}_{34}=2\text{i}\beta {\mu }^{\left(j\right)}{\left({k}_{\mathrm{s}}^{\left(j\right)}\right)}^{2}{J}_{0\mathrm{s}}^{{\prime\prime}},\end{equation} \tag{ B.22 }$$

$$\begin{equation}{m}_{41}=2\text{i}\beta {\mu }^{\left(j\right)}{k}_{\mathrm{p}}^{\left(j\right)}{H}_{0\mathrm{p}}^{\prime },\end{equation} \tag{ B.23 }$$

$$\begin{equation}{m}_{42}=2\text{i}\beta {\mu }^{\left(j\right)}{k}_{\mathrm{p}}^{\left(j\right)}{J}_{0\mathrm{p}}^{\prime },\end{equation} \tag{ B.24 }$$

$$\begin{equation}{m}_{43}={\mu }^{\left(j\right)}{k}_{\mathrm{s}}^{\left(j\right)}\left({\left({k}_{\mathrm{s}}^{\left(j\right)}\right)}^{2}-{\beta }^{2}\right){H}_{0\mathrm{s}}^{\prime },\end{equation} \tag{ B.25 }$$

$$\begin{equation}{m}_{44}={\mu }^{\left(j\right)}{k}_{\mathrm{s}}^{\left(j\right)}\left({\left({k}_{\mathrm{s}}^{\left(j\right)}\right)}^{2}-{\beta }^{2}\right){J}_{0\mathrm{s}}^{\prime }.\end{equation} \tag{ B.26 }$$

Again, we formulate the waveguiding condition by relating vectors ${\tilde {\mathbf{u}}}^{\left(1\right)}$ and ${\tilde {\mathbf{u}}}^{\left(3\right)}$ via M matrices, and require that the innermost (outermost) layer has no outgoing (incoming) components. The resulting transcendental equation can also be solved numerically using the GRPF algorithm.