Nonlinear Phononics and Coupled Mode Theory

Abstract

Understanding phononic wave propagation in nonlinear and coupled media is inevitable in designing devices to harness phonons. In photonics, the nonlinear Schrodinger equation (NSE) and coupled mode theory (CMT) are extensively used to manipulate the propagation of electromagnetic (EM) waves. Since phonons as the quasi-particles related to mechanical vibrations are similar to their photonic counterparts (photons), inspired by EM, the NSE and CMT are developed for elastic or phononic waves. In the first section, from the novel point of view, the nonlinearity is dealt with and consequently applied to the equation of motion to derive the NSE. In the following, a hard limiter is introduced as an application of nonlinear wave propagation. The coupled-mode equations are obtained for the elastic slab waveguides in the second section. These equations are applied to explain the wave propagation in two parallel and nonparallel waveguides.

1. Introduction

Phonons, the quantas of mechanical vibration, play a significant role not only in phononic systems like phononic crystals but also in optoelectronic and optomechanical devices where electron-phonon and photon-phonon interactions have an undeniable impact on the operation of devices. Furthermore, they can influence the operation of photonic devices and nano and micro-electromechanical systems (N/MEMS). Raman and Brillouin scattering are the most famous phonon-involved phenomena. Cavity optomechanics, superconducting circuits, cell manipulation, quantum acoustics, and communication devices are some hot topic issues utilizing phononic waves [1]. The electron-phonon interaction affects the mobility of carriers in semiconductors, the electrical resistance of metals, carrier density in graphene, and optical absorption of the indirect semiconductors [2, 3]. Moreover, the phonon-photon interaction can be manipulated to change the group velocity used in the integrated quantum optomechanical memories and signal processing [4, 5]. Due to these numerous applications and effects, harnessing phonons is of interest to engineers and physicists. The nonlinear media and coupled mode theory are significant in controlling wave propagation regardless of its type, which can be electronic, photonic, or phononic. This book chapter is generally classified into nonlinear phononics and coupled mode theory.

2. Nonlinear phononics

Nonlinearity has broad applications in electronics and photonics. Designing Schmitt triggers in electronics, full-optical analog to digital converters, and optical microscopy are some examples of the use of nonlinearity [6, 7, 8, 9, 10]. Generally, a nonlinear medium is referred to as a medium in which the basic parameters depend on the amplitude of the transmitted signal. The nonlinear media provide the possibility of controlling wave propagation via other waves. In electromagnetics and photonics, the nonlinear Schrödinger equation (NSE) was introduced and has been extensively utilized to investigate nonlinear phenomena such as second harmonic generation, four-wave mixing, and modulation instability [11]. In this section, The NSE is introduced for phononic or elastic waves in an elastic medium. Furthermore, this equation is applied to describe the operation of a hard limiter in which the amplitude of the output wave is confined to the desired value.

2.1 Elastic waves

As the electromagnetic wave equation describes the optical wave propagation, the phononic waves in solids are governed by the equation of motion given by Eq. (1) in which $\rho$, $u_i$ and ${\sigma }_{\mathit{ij}}$ are mass density, displacement, and stress tensor [12].

$$\rho {\ddot{u}}_i=\frac{\partial {\sigma }_{\mathit{ij}}}{\partial x_j}$$E1

The stress tensor is related to the strain tensor via Eq. (2), where $c_{\mathit{\text{ijkl}}}$ is the stiffness tensor and $e_{\mathit{kl}}$ is the strain tensor defined by Eq. (3) [12].

$${\sigma }_{\mathit{ij}}=c_{\mathit{\text{ijkl}}}{e}_{\mathit{kl}}$$E2

$$e_{\mathit{ij}}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$E3

Using Voigt notation, the stiffness tensor for cubic elasticity is expressed by a 6 × 6 matrix in which $\lambda$ and $\mu$ are Lamb parameters [12].

$$c_{\mathit{ij}}=\left(\begin{array}{cccccc}\lambda +2\mu & \lambda & 0& 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ 0& 0& 0& \lambda & 0& 0\\ 0& 0& 0& 0& \lambda & 0\\ 0& 0& 0& 0& 0& \lambda \end{array}\right)$$E4

The wave equation for the displacement is obtained using Eqs. (1) to (4) and given by [13]:

$$\mu {\nabla }^2\overrightarrow{u}+\left(\lambda +\mu \right)\nabla \left(\nabla .\overrightarrow{u}\right)=\rho \frac{{\partial }^2\overrightarrow{u}}{\partial t^2}$$E5

2.2 Mass density as the fundamental parameter of nonlinearity

Lennard-Jones potential, shown in Figure 1, is a simple mathematical model for the interaction between neutral atoms and molecules. According to the mass-spring model, an atom located in the equilibrium point with the minimum potential energy is considered a particle with a certain mass connected to the fixed particle at the origin through a spring. Around this equilibrium point, the potential can be approximated quadratically; hence, the force, which is the gradient of potential energy, named by Hooke’s law, is linear. Increasing the force on the particle causes increasing displacement; hence, more force changes the particle’s position more than before. Therefore, the quadratic approximation will no longer be valid. It leads to the nonlinear relationship between displacement and force. The greater the displacement from the equilibrium distance, the more the compression or expansion of the spring. Consequently, the total mass in unit length is increased or decreased due to the flux of the displacement field [14].

Figure 1.

Lennard-Jones potential energy.

2.2.1 Mathematical modeling of nonlinearity

The deformation is defined as Eq. (6), where ΔV, V_0, and u are volume change, initial volume, and displacement, respectively [14].

$$\frac{\mathit{\Delta V}}{V_0}=\nabla .\overrightarrow{u}$$E6

Using the conservation of mass and Eq. (6), it can be obtained:

$${\rho }_0{V}_0=\mathit{\rho V}=\rho V_0\left(1+\nabla .\overrightarrow{u}\right)$$E7

where $V=V_0+\mathit{\Delta V}$ and $\rho$ is the perturbed mass density. Therefore $\rho$ is expressed as:

$$\rho =\frac{{\rho }_0}{1+\nabla .\overrightarrow{u}}$$E8

Since the deformation is small compared to the initial volume, the perturbed mass density can be described as:

$$\rho ={\rho }_0\left(1-\nabla .\overrightarrow{u}+{\left(\nabla .\overrightarrow{u}\right)}^2+\cdots \right)$$E9

The displacement can be represented by scalar ($\phi$) and vector ($\psi$) potential as the following equation [13, 15, 16].

$$u=\nabla \phi +\nabla \times \overrightarrow{\psi }$$E10

Applying Eq. (10) to Eq. (5) leads to the wave equations for scalar and vector potential [13, 16].

$${\nabla }^2\phi =\left(\frac{{\rho }_0}{\lambda +2\mu }\right)\frac{{\partial }^2\phi }{\partial t^2}$$E11

$${\nabla }^2\overrightarrow{\psi }=\left(\frac{{\rho }_0}{\mu }\right)\frac{{\partial }^2\overrightarrow{\psi }}{\partial t^2}$$E12

Eq. (10) indicates that $\nabla .\overrightarrow{u}={\nabla }^2\phi$ and hence, Eq. (9) can be expressed by scalar potential as:

$$\rho ={\rho }_0\left(1-{\nabla }^2\phi +{\left({\nabla }^2\phi \right)}^2+\cdots \right)$$E13

2.3 Nonlinear schrödinger equation (NSE)

The wave equation for scalar potential is expressed as Eq. (13), in which $n_0=\sqrt{\frac{{\rho }_0}{\lambda +2\mu }}$ is considered the refractive index similar to EM waves [14].

$${\nabla }^2\phi =-n_0^2{\omega }^2\phi$$E14

The scalar potential is assumed to be $\phi \left(r,t\right)=\left(1/2\right)\left\{F\left(x,y\right)A(z,t)exp\left[i\left({\beta }_{0}z-\mathit{\omega t}\right)\right]+c.c\right\}$ in the time domain and $\phi \left(r,\omega -{\omega }_0\right)=\left(1/2\right)\left\{F\left(x,y\right)A(z,\omega -{\omega }_0)expi\left({\beta }_{0}z\right)+c.c\right\}$ in the frequency domain. A and F are slowly varying envelope and mode shapes, respectively. ${\beta }_0$ is the wave number and ${\omega }_0$ is the central frequency of the input pulse. Inserting this assumed response in the wave equation results in [14]:

$$\frac{1}{A}\left(\frac{{\partial }^{2}A}{\partial z^2}+2j{\beta }_0\frac{\partial A}{\partial z}-{\beta }_0^{2}A\right)+\frac{1}{F}\left(\frac{{\partial }^{2}F}{\partial x^2}+\frac{{\partial }^{2}F}{\partial y^2}\right)+\frac{{\rho }_0}{\lambda +2\mu }{\omega }^2=0$$E15

Two eigenvalue equations for A and F (Eqs. 16 and 17) are obtained using the separation of variable method. In Eq. (15), the second derivation of amplitude is neglected because of its slow variation.

$$\frac{{\partial }^{2}F}{\partial x^2}+\frac{{\partial }^{2}F}{\partial y^2}+\left(\frac{{\rho }_0}{\lambda +2\mu }{\omega }^2-{\tilde{\beta }}^2\right)F=0$$E16

$$2j{\beta }_0\frac{\partial A}{\partial z}+\left({\tilde{\beta }}^2-{\beta }_0^2\right)A=0$$E17

The perturbed mass density ($\rho$) equals ${\rho }_0\left(1-{\nabla }^2\phi \right)$ in the first order perturbation. This mass density causes the refractive index to be perturbed:

$$n=n_0+\mathit{\Delta n}\triangleq \sqrt{\frac{{\rho }_0+\mathit{\Delta \rho }}{\lambda +2\mu }}$$E18

As $\mathit{\Delta \rho }\ll {\rho }_0$ the perturbation in the refractive index is obtained by Eq. (19), which changes the eigenvalue ($\tilde{\beta }$) (Eq. 20). The value of this alteration is calculated by Eq. (21) [11].

$$\mathit{\Delta n}=-\frac{1}{2}{\left(\frac{{\rho }_0}{\lambda +2\mu }\right)}^{\frac{3}{2}}{\omega }_0^2\left|\phi \right|$$E19

$$\tilde{\beta }\left(\omega \right)=\beta \left(\omega \right)+\mathit{\Delta \beta }$$E20

$$\mathit{\Delta \beta }=-\frac{1}{2}{\left(\frac{{\rho }_0}{\lambda +2\mu }\right)}^{\frac{3}{2}}{\omega }_0^3\frac{1}{A_{\mathit{eff}}}\left|A\right|$$E21

$A_{\mathit{eff}}$ is defined as Eq. (22). It is the effective area where most wave energy is confined.

$$A_{\mathit{eff}}=\frac{{\left(\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{\left|F\right|}^2\mathit{\text{dxdy}}\right)}^2}{\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}{\left|F\right|}^3\mathit{\text{dxdy}}}$$E22

The eigenvalue equation for the amplitude is written as Eq. (23) using Eq. (21).

$$\frac{\partial A}{\partial z}=j\left(\beta \left(\omega \right)+\mathit{\Delta \beta }-{\beta }_0\right)A$$E23

According to the dependence of the wave number on the angular frequency, it can be expanded around the central frequency of the input pulse in the form of Eq. (24), where $1/{\beta }_1$ is the group velocity.

$$\beta \left(\omega \right)={\beta }_0+\left(\omega -{\omega }_0\right){\beta }_1+\frac{1}{2}{\left(\omega -{\omega }_0\right)}^2{\beta }_2+\frac{1}{6}{\left(\omega -{\omega }_0\right)}^3{\beta }_3+\dots$$E24

Inserting the above equation in Eq. (23) and applying the inverse Fourier transform results in the nonlinear Schrödinger equation (NSE) [14]:

$$\frac{\partial A}{\partial z}=-\left({\beta }_1\frac{\partial A}{\partial t}+i{\beta }_2\frac{{\partial }^{2}A}{\partial t^2}+\dots \right)-i{{\omega }_0}^3\left({\left(\frac{{\rho }_0}{\lambda +2\mu }\right)}^{\frac{3}{2}}\left(\frac{1}{A_{\mathit{eff}}}\right)\right)\left|A\right|A+\dots$$E25

The NSE includes both the dispersion and nonlinear effects. ${\beta }_1$, ${\beta }_2$, and higher orders $\beta$ indicate the first, second, and higher order dispersions. The term containing the effective area represents nonlinearity. The NSE clearly expresses how dispersion and nonlinearity affect wave amplitude in a specific position during wave propagation.

2.4 Hard limiter

NSE can be utilized to describe nonlinear phenomena such as wave mixing and phase modulation. Wave mixing is used to design phononic limiters. A limiter confines the output to a particular value when the input exceeds a threshold value. As the mass density is dynamically changed due to the wave propagation, the wave number is considered $\beta \left(\omega \right)=\omega \widehat{\rho }\left(\omega \right)$ with $\widehat{\rho }\left(\omega \right)=\sqrt{\rho \left(\omega \right)/\left(\lambda +2\mu \right)}$. $\widehat{\rho }\left(\omega \right)$ is introduced as the elastic wave refractive index using duality between elastic and electromagnetic waves [17].

In the 1-dimensional structure depicted in Figure 2, Eq. (5) is expressed as:

Figure 2.

The 1D periodic structure includes two types of layers with lengths of l₁ and l₂. ρ₀ and ρ₀⁽¹⁾ are the linear and first nonlinear parts of mass density, respectively [14].

$$\frac{{\partial }^{2}u\left(z\right)}{\partial z^2}=-{\left[\omega \widehat{\rho }\left(\omega \right)\right]}^{2}u\left(z\right)$$E26

The nonlinear mass density can be written as:

$$\rho ={\rho }_0\left(1-\frac{u}{c_0}+{\left(\frac{u}{c_0}\right)}^2+\dots \right)$$E27

where c₀ is the lattice constant in the direction of displacement. The refractive index in Eq. (26) for the structure of Figure 2 can be expressed as [18]:

$$\begin{array}{lll}{\left[\widehat{\rho }\left(z\right)\right]}^2& =& {\widehat{\rho }}_{\mathit{av}}^2+2{\widehat{\rho }}_{\mathit{av}}\left({\widehat{\rho }}_1-{\widehat{\rho }}_2\right)f\left(z\right)+i2{\widehat{\rho }}_{\mathit{av}}\left[{\kappa }_0+\left({\kappa }_1-{\kappa }_2\right)f\left(z\right)\right]\\ & +& 2{\widehat{\rho }}_{\mathit{av}}\left[{\widehat{\rho }}_{\mathit{NL}}+\left({\widehat{\rho }}_{1\mathit{NL}}-{\widehat{\rho }}_{2\mathit{NL}}\right)f\left(z\right)\right]\left|u\right|\end{array}$$E28

where ${\widehat{\rho }}_1=\sqrt{{\rho }_{01}/\left({\lambda }_1+2{\mu }_1\right)}$ and ${\widehat{\rho }}_2=\sqrt{{\rho }_{02}/\left({\lambda }_2+2{\mu }_2\right)}$ are linear refractive indices and ${\widehat{\rho }}_{\mathit{av}}$ is the average refractive index given by:

$${\widehat{\rho }}_{\mathit{av}}=\left(1-\frac{l_2}{L}\right){\widehat{\rho }}_1+\left(\frac{l_2}{L}\right){\widehat{\rho }}_2$$E29

$f\left(z\right)$ is the periodic function defined by Eq. (30) in one period.

$$f\left(z\right)=\left\{\begin{array}{l}\frac{l_2}{L},0<z<l_1\\ \frac{l_2}{L}-1,l_1<z<L\end{array}\right.$$E30

κ₁ and κ₂ are the imaginary parts of refractive indices, and κ₀ is their average. Similar to the linear parts, Eq. (31) defines the average nonlinear refractive index ${\widehat{\rho }}_{\mathit{NL}}$. ${\widehat{\rho }}_{1\mathit{NL}}$ and ${\widehat{\rho }}_{2\mathit{NL}}$ are the nonlinear indices obtained through the Taylor expansion of the refractive index. These indices are expressed via Eq. (33) considering the nonlinear relation for mass density given by Eq. (32), [14].

$${\widehat{\rho }}_{\mathit{Nl}}=\left(1-\frac{l_2}{L}\right){\widehat{\rho }}_{1\mathit{NL}}+\left(\frac{l_2}{L}\right){\widehat{\rho }}_{2\mathit{NL}}$$E31

$$\rho ={\rho }_0+{\rho }_0^{\left(1\right)}u$$E32

$${\widehat{\rho }}_{1\mathit{NL}}=\frac{{\rho }_{01}^{\left(1\right)}}{2\sqrt{{\rho }_{01}\left({\lambda }_1+2{\mu }_1\right)}}$$E33

$${\widehat{\rho }}_{2\mathit{NL}}=\frac{{\rho }_{02}^{\left(1\right)}}{2\sqrt{{\rho }_{02}\left({\lambda }_2+2{\mu }_2\right)}}$$E34

The displacement field is the sum of forward and backward waves given by: $u\left(z\right)=A\left(z\right)e^{\mathit{ikz}}+B\left(z\right)e^{-\mathit{ikz}}$, in which $k=\omega {\widehat{\rho }}_{\mathit{av}}$. Inserting this equation and Eq. (28) in Eq. (26) and using Eq. (35), which is the Fourier transform of f(z), results in Eq. (36) in which the slowly varying approximation is applied [18].

$$f\left(z\right)=-\sum _{m\ne 0}{e}^{\mathit{im\pi }\raisebox{1ex}{$l_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\mathit{m\pi }\raisebox{1ex}{$l_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)z}$$E35

$$\begin{array}{ll}-i\left(\frac{\mathit{dA}\left(z\right)}{\mathit{dz}}{e}^{\mathit{ikz}}-\frac{\mathit{dB}\left(z\right)}{\mathit{dz}}{e}^{-\mathit{ikz}}\right)& =-\omega \left({\widehat{\rho }}_1-{\widehat{\rho }}_2\right)A\left(z\right)\sum _{m\ne 0}{e}^{\raisebox{1ex}{$\mathit{im\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\raisebox{1ex}{$\mathit{m\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.+k\right)z}\\ & -\mathit{i\omega }\left({\kappa }_1-{\kappa }_2\right)A\left(z\right)\sum _{m\ne 0}{e}^{\raisebox{1ex}{$\mathit{im\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\raisebox{1ex}{$\mathit{m\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.+k\right)z}\\ & -\omega \left({\widehat{\rho }}_{1\mathit{NL}}-{\widehat{\rho }}_{2\mathit{NL}}\right)\left|u\left(z\right)\right|A\left(z\right)\sum _{m\ne 0}{e}^{\raisebox{1ex}{$\mathit{im\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\raisebox{1ex}{$\mathit{m\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.+k\right)z}\\ & -\omega \left({\widehat{\rho }}_1-{\widehat{\rho }}_2\right)B\left(z\right)\sum _{m\ne 0}{e}^{\raisebox{1ex}{$\mathit{im\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\raisebox{1ex}{$\mathit{m\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.-k\right)z}\\ & -\mathit{i\omega }\left({\kappa }_1-{\kappa }_2\right)B\left(z\right)e^{\raisebox{1ex}{$\mathit{im\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\raisebox{1ex}{$\mathit{m\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.-k\right)z}\\ & -\omega \left({\widehat{\rho }}_{1\mathit{NL}}-{\widehat{\rho }}_{2\mathit{NL}}\right)\left|u\left(z\right)\right|B\left(z\right)\sum _{m\ne 0}{e}^{\raisebox{1ex}{$\mathit{im\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.}\frac{sin\left(\raisebox{1ex}{$\mathit{m\pi }{l}_2$}\!\left/ \!\raisebox{-1ex}{$L$}\right.\right)}{\mathit{m\pi }}{e}^{i\left(\raisebox{1ex}{$2\mathit{m\pi }$}\!\left/ \!\raisebox{-1ex}{$L$}\right.-k\right)z}\\ & +\omega \left(i{\kappa }_0+{\widehat{\rho }}_{\mathit{NL}}\left|u\right|\right)\left(A\left(z\right)e^{\mathit{ikz}}+B\left(z\right)e^{-\mathit{ikz}}\right)\end{array}$$E36

Assuming the wavelength around the first stop band ($k\approx \raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$L$}\right.$), the phase-matching condition leads to the coupled equations for the amplitudes of forward and backward waves given by Eq. (36), [8, 18].

$$\begin{array}{lll}i\frac{\mathit{dA}\left(z\right)}{\mathit{dz}}& =& \omega \left[\left({\widehat{\rho }}_1-{\widehat{\rho }}_2\right)+\left({\widehat{\rho }}_{1\mathit{NL}}-{\widehat{\rho }}_{2\mathit{NL}}\right)\left|u\left(z\right)\right|\right]e^{\frac{\mathit{i\pi }{l}_2}{L}}\frac{sin\left(\frac{\pi l_2}{L}\right)}{\pi }{e}^{i\left(\frac{2\pi }{L}-2k\right)z}B\left(z\right)\\ & & -{\widehat{\rho }}_{\mathit{NL}}\left|u\left(z\right)\right|A\left(z\right)\end{array}$$E37

$$\begin{array}{lll}i\frac{\mathit{dB}\left(z\right)}{\mathit{dz}}& =& \omega \left[\left({\widehat{\rho }}_2-{\widehat{\rho }}_1\right)+\left({\widehat{\rho }}_{2\mathit{NL}}-{\widehat{\rho }}_{1\mathit{NL}}\right)\left|u\left(z\right)\right|\right]e^{\frac{-\mathit{i\pi }{l}_2}{L}}\frac{sin\left(\frac{\pi l_2}{L}\right)}{\pi }{e}^{-i\left(\frac{2\pi }{L}-2k\right)z}A\left(z\right)\\ & & +{\widehat{\rho }}_{\mathit{NL}}\left|u\left(z\right)\right|B\left(z\right)\end{array}$$E38

These equations will be solved for the lengths given by Eqs. (39) and (40).

$$l_1=\frac{\lambda }{2\left({\widehat{\rho }}_1-{\widehat{\rho }}_2\frac{{\widehat{\rho }}_{1\mathit{NL}}}{{\widehat{\rho }}_{2\mathit{NL}}}\right)}$$E39

$$l_2=\frac{\lambda }{2\left({\widehat{\rho }}_2-{\widehat{\rho }}_1\frac{{\widehat{\rho }}_{2\mathit{NL}}}{{\widehat{\rho }}_{1\mathit{NL}}}\right)}$$E40

The similarity of the linear characteristic of two layers in the structure of Figure 2 means that the propagating medium is homogeneous in the absence of nonlinearity; hence, there is no reflected wave. However, the different nonlinear features of the two layers cause the medium to be inhomogeneous during wave propagation. The intensity of the reflected wave corresponds to the value of the nonlinear part of the refractive index.

The linear and nonlinear parts of the refractive indices are assumed to be ${\widehat{\rho }}_1={\widehat{\rho }}_2=1.5\times {10}^{-4}$, ${\widehat{\rho }}_{1\mathit{NL}}={10}^5$, and ${\widehat{\rho }}_{2\mathit{NL}}=-{10}^5$ to investigate the operation of the limiter. The transmission for different number of layers at the input intensity of 1.44 and different intensities at N = 50 are depicted in Figure 3. The thickness of layers is set to hinder the transmission of 10 MHz. The transmission is decreased due to the increasing number of layers and intensity because of boosting the nonlinear effect [14].

Figure 3.

${\widehat{\rho }}_{\mathit{1}}={\widehat{\rho }}_{\mathit{2}}=\mathit{1.5}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{1}\mathit{NL}}={\mathit{10}}^{\mathit{5}}$, and${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-{\mathit{10}}^{\mathit{5}}$. (a) for different number of layers at input intensity of 1.44. (b) for different intensities at N = 50 [14].

Characteristic curves of the limiter are shown in Figure 4. Figure 4a illustrates that more layers result in a better limiting process. For a nonlinear index difference (NID) of 2 × 10⁵, 100 layers are required to confine the output to 1. Furthermore, as seen in Figure 4b, the more the NID, the better the limiting. This figure indicates that NID = 3 × 10⁵ is the best for 100 layers to limit the output to 1. Moreover, Figure 4c represents that NID can be adjusted so that the output is confined to the desirable value. In this figure, NIDs are selected so that the output is limited to 0.64, 0.81, and 1.

Figure 4.

The characteristic curves for the phononic hard limiter. (a) Different number of layers when the linear and nonlinear parts of the refractive index are ${\widehat{\rho }}_{\mathit{1}}=\mathit{1.5}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{2}}=\mathit{1.7}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{1}\times {\mathit{10}}^{\mathit{5}}$, and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{1}\times {\mathit{10}}^{\mathit{5}}$. (b) N = 100. Dotted line: ${\widehat{\rho }}_{\mathit{1}}=\mathit{1.5}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{2}}=\mathit{1.6}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{0.5}\times {\mathit{10}}^{\mathit{5}}$, and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{0.5}\times {\mathit{10}}^{\mathit{5}}$. Dashed line: ${\widehat{\rho }}_{\mathit{1}}=\mathit{1.5}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{2}}=\mathit{1.7}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{1}\times {\mathit{10}}^{\mathit{5}}$, and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{1}\times {\mathit{10}}^{\mathit{5}}$. Solid line: ${\widehat{\rho }}_{\mathit{1}}=\mathit{1.5}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{2}}=\mathit{1.8}\times {\mathit{10}}^{-\mathit{4}}$, ${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{1.5}\times {\mathit{10}}^{\mathit{5}}$, and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{1.5}\times {\mathit{10}}^{\mathit{5}}$. (c) N = 100, ${\widehat{\rho }}_{\mathit{1}}=\mathit{1.4}\times {\mathit{10}}^{-\mathit{4}}$, and ${\widehat{\rho }}_{\mathit{2}}=\mathit{1.6}\times {\mathit{10}}^{-\mathit{4}}$. Dotted line: The output is limited to 1 with ${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{1}\times {\mathit{10}}^{\mathit{5}}$ and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{1}\times {\mathit{10}}^{\mathit{5}}$. Dashed line: The limited level is 0.81 with${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{1.1}\times {\mathit{10}}^{\mathit{5}}$ and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{1.1}\times {\mathit{10}}^{\mathit{5}}$. Solid line: The limited level is 0.64 with${\widehat{\rho }}_{\mathit{1}\mathit{NL}}=\mathit{1.25}\times {\mathit{10}}^{\mathit{5}}$ and ${\widehat{\rho }}_{\mathit{2}\mathit{NL}}=-\mathit{1.25}\times {\mathit{10}}^{\mathit{5}}$ [14].

2.5 Future challenge

Metamaterials are required to realize the negative refractive index. These artificial materials were introduced in 1968 for electromagnetic waves [19]. More efforts were made to implement these materials afterward [20, 21, 22, 23, 24, 25, 26, 27]. Furthermore, the nonlinear properties of metamaterials are also of interest [28, 29, 30]. Phononic metamaterials can also be designed so that their effective Young’s modulus or mass density is negative [31, 32, 33]. Realization of the phononic hard limiter requires the nonlinear metamaterial. Research in this area can be an exciting subject and help design a variety of devices in phononics. Furthermore, a phononic limiter can be utilized to manipulate phonon-photon interaction in optomechanics and electron-phonon interaction in optoelectronic or electronic devices.

3. Phononic wave coupling in slab waveguides

Understanding the wave coupling is crucial to design devices based on wave propagation. The Coupled mode theory (CMT) discusses the energy exchange between different modes or coupled waveguides [34, 35, 36]. CMT in the base of Machzender interferometer, microsphere resonators, couplers, switches, and filters in photonics [37]. Herein, the CMT is applied to the elastic slab waveguide to obtain the coupled mode equations (CME) for shear and Rayleigh-Lame modes inspired by electromagnetic wave propagation.

3.1 Coupled mode equations

3.1.1 Shear waves

The wave equation for the shear waves (SH) is expressed by Eq. (37) for a waveguide orientated in the direction of z (Figure 5), where ${\epsilon }^0=\raisebox{1ex}{$\rho $}\!\left/ \!\raisebox{-1ex}{$\mu $}\right.$. This equation is for an unperturbed system containing a single waveguide. The solution for the shear wave equation is given by Eq. (42), in which a_m, U_m(x), and β_m are the amplitude, mode shape, and wave number, respectively, and m is the mode number.

Figure 5.

The elastic slab waveguide. a and b are the width of the core and cladding regions [38].

$$\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}\right)u_y+{\epsilon }^0{\omega }^2{u}_y=0$$E41

$$u_y=\sum _m{a}_{\mathit{ym}}{U}_{\mathit{ym}}\left(x\right)e^{-j{\beta }_{m}z}$$E42

The wave equation for the mode shape is obtained via inserting Eq. (42) in Eq. (41) and given by:

$$\frac{{\partial }^2{U}_{\mathit{ym}}\left(x\right)}{\partial x^2}+\left({\epsilon }^0{\omega }^2-{\beta }_m^2\right)U_{\mathit{ym}}\left(x\right)=0$$E43

For the waveguides that interchange energy, assuming they are not very close together, the perturbed wave equation is written as Eq. (44). ${\epsilon }^1$ represents the perturbation.

$$\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}\right)u_y+\left({\epsilon }^0+{\epsilon }^1\right){\omega }^2{u}_y=0$$E44

The perturbation causes the amplitude not to be constant in the propagation direction; hence, the displacement should be written as Eq. (45), where $a_{\mathit{ym}}\left(z\right)$ is the slowly varying amplitude.

$$u_y=\sum _m{a}_{\mathit{ym}}\left(z\right)U_{\mathit{ym}}\left(x\right)e^{-j{\beta }_{m}z}$$E45

Inserting Eq. (45) in Eq. (44) results in Eq. (46).

$$\begin{array}{l}\sum _m{a}_{\mathit{ym}}\left(z\right)\left(\frac{{\partial }^2{U}_{\mathit{ym}}\left(x\right)}{\partial x^2}+\left({\epsilon }^0{\omega }^2-{\beta }_m^2\right)U_{\mathit{ym}}\left(x\right)\right)e^{-j{\beta }_{m}z}+\\ \sum _m\left(\frac{{\partial }^2{a}_{\mathit{ym}}\left(z\right)}{\partial z^2}-j2{\beta }_m\frac{\partial a_{\mathit{ym}}\left(z\right)}{\partial z}\right)U_{\mathit{ym}}\left(x\right)e^{-j{\beta }_{m}z}\\ =-{\epsilon }^1{\omega }^2\sum _m{a}_{\mathit{ym}}\left(z\right)U_{\mathit{ym}}\left(x\right)e^{-j{\beta }_{m}z}\end{array}$$E46

Eq. (43) specifies that the first term in Eq. (46) is zero. Furthermore, the second derivative of amplitude is neglected because of its slow variation. The orthogonality of modes leads to Eq. (47) with the coupling coefficient obtained by Eq. (48) [38].

$$\frac{\partial a_{\mathit{yn}}\left(z\right)}{\partial z}=-j\frac{{\omega }^2}{2{\beta }_n}\sum _m{\kappa }_{\left(\mathit{SH}\right)n,m}{a}_{\mathit{ym}}\left(z\right)e^{-j\left({\beta }_m-{\beta }_n\right)z}$$E47

$${\kappa }_{\left(\mathit{SH}\right)n,m}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{yn}}^{\ast }\left(x\right){\epsilon }^1{U}_{\mathit{ym}}\left(x\right)\mathit{dx}}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{yn}}^{\ast }\left(x\right)U_{\mathit{yn}}\left(x\right)\mathit{dx}}$$E48

3.1.2 Rayleigh-lame waves

Rayleigh-Lame (RL) modes have displacement in the x and z directions for the waveguide depicted in Figure 5. Therefore, the wave equation for the RL modes is expressed by [15]:

$$\mu {\nabla }^2{u}_i+\left(\lambda +\mu \right)\frac{\partial }{\partial x_i}\sum _{i=1,2}\frac{\partial u_i}{\partial x_i}+{\mathit{\rho \omega }}^2{u}_i=0$$E49

where i = 1 and 2 indicate the direction of x and z, respectively. Consequently, the wave equations for displacement components are given by:

$${\eta }_s^0\left(\frac{{\partial }^2{u}_x}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial z^2}\right)+{\eta }_l^0\left(\frac{{\partial }^2{u}_x}{\partial x^2}+\frac{{\partial }^2{u}_z}{\partial x\partial z}\right)+{\omega }^2{u}_x=0$$E50

$${\eta }_s^0\left(\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2{u}_z}{\partial z^2}\right)+{\eta }_l^0\left(\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial x\partial z}\right)+{\omega }^2{u}_z=0$$E51

where ${\eta }_s^0=\mu /\rho$ and ${\eta }_l^0=\left(\lambda +\mu \right)/\rho$ are unperturbed parameters of the waveguide. Coupling results in the perturbation in ${\eta }_s^0$ and ${\eta }_l^0$; hence, the above equations change to Eqs. (52) and (53), where ${\eta }_s^1$ and ${\eta }_l^1$ are the perturbation terms.

$$\left({\eta }_s^0+{\eta }_s^1\right)\left(\frac{{\partial }^2{u}_x}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial z^2}\right)+\left({\eta }_l^0+{\eta }_l^1\right)\left(\frac{{\partial }^2{u}_x}{\partial x^2}+\frac{{\partial }^2{u}_z}{\partial x\partial z}\right)+{\omega }^2{u}_x=0$$E52

$$\left({\eta }_s^0+{\eta }_s^1\right)\left(\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2{u}_z}{\partial z^2}\right)+\left({\eta }_l^0+{\eta }_l^1\right)\left(\frac{{\partial }^2{u}_z}{\partial x^2}+\frac{{\partial }^2{u}_x}{\partial x\partial z}\right)+{\omega }^2{u}_z=0$$E53

Similar to the shear waves, the displacement components (u_x and u_z) are assumed to be:

$$u_x=\sum _m{a}_{\mathit{xm}}\left(z\right)U_{\mathit{xm}}\left(x\right)e^{-j{\beta }_{m}z}$$E54

$$u_x=\sum _m{a}_{\mathit{xm}}\left(z\right)U_{\mathit{xm}}\left(x\right)e^{-j{\beta }_{m}z}$$E55

CMEs for RL modes are expressed by Eqs. (56) and (57) using the solutions mentioned above, the orthogonality of modes, and the slow variation of amplitude. Eqs. (58) and (59) give the coupling coefficients in CMEs [38].

$$2j{\beta }_n\frac{\partial a_{\mathit{xn}}}{\partial z}+{\beta }_n^2{a}_{\mathit{xn}}=\sum _m\left[{\kappa }_{\left(\mathit{xx}\right)n,m}{a}_{\mathit{xm}}+{\kappa }_{\left(\mathit{xz}\right)n,m}^{\prime }\left(\frac{\partial a_{\mathit{zm}}}{\partial z}-j{\beta }_m{a}_{\mathit{zm}}\right)\right]e^{-j\left({\beta }_m-{\beta }_n\right)z}$$E56

$$2j{\beta }_n\frac{\partial a_{\mathit{zn}}}{\partial z}+{\beta }_n^2{a}_{\mathit{zn}}=\sum _m\left[{\kappa }_{\left(\mathit{zz}\right)n,m}{a}_{\mathit{zm}}+{\kappa }_{\left(\mathit{zx}\right)n,m}^{\prime }\left(\frac{\partial a_{\mathit{xm}}}{\partial z}-j{\beta }_m{a}_{\mathit{xm}}\right)\right]e^{-j\left({\beta }_m-{\beta }_n\right)z}$$E57

$${\kappa }_{\left(\mathit{xx}\right)n,m}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{xn}}^{\ast }\left(x\right)\left(1+\frac{{\eta }_l^1\left(x\right)}{{\eta }_s^1\left(x\right)}\right)\frac{{\partial }^2{U}_{\mathit{xm}}\left(x\right)}{\partial x^2}\mathit{dx}}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{xn}}^{\ast }\left(x\right)U_{\mathit{xn}}\left(x\right)\mathit{dx}}$$E58

$${\kappa }_{\left(\mathit{xz}\right)n,m}^{\prime }=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{xn}}^{\ast }\left(x\right)\left(\frac{{\eta }_l^1\left(x\right)}{{\eta }_s^1\left(x\right)}\right)\frac{\partial U_{\mathit{zm}}\left(x\right)}{\partial x}\mathit{dx}}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{xn}}^{\ast }\left(x\right)U_{\mathit{xn}}\left(x\right)\mathit{dx}}$$E59

$${\kappa }_{\left(\mathit{zz}\right)n,m}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{zn}}^{\ast }\left(x\right)\left(1+\frac{{\eta }_l^1\left(x\right)}{{\eta }_s^1\left(x\right)}\right)\frac{{\partial }^2{U}_{\mathit{zm}}\left(x\right)}{\partial x^2}\mathit{dx}}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{zn}}^{\ast }\left(x\right)U_{\mathit{zn}}\left(x\right)\mathit{dx}}$$E60

$${\kappa }_{\left(\mathit{zx}\right)n,m}^{\prime }=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{zn}}^{\ast }\left(x\right)\left(\frac{{\eta }_l^1\left(x\right)}{{\eta }_s^1\left(x\right)}\right)\frac{\partial U_{\mathit{xm}}\left(x\right)}{\partial x}\mathit{dx}}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{\mathit{zn}}^{\ast }\left(x\right)U_{\mathit{zn}}\left(x\right)\mathit{dx}}$$E61

The modal analysis should be performed to obtain the mode shapes for calculating the coupling coefficients. Eqs. (56) and (57) can be simplified by assuming ${\widehat{a}}_{\mathit{xn}}\triangleq a_{\mathit{xn}}{e}^{\raisebox{1ex}{$-j{\beta }_{n}z$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$, ${\widehat{a}}_{\mathit{zn}}\triangleq a_{\mathit{zn}}{e}^{\raisebox{1ex}{$-j{\beta }_{n}z$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$, ${\widehat{a}}_{\mathit{xm}}\triangleq a_{\mathit{xm}}{e}^{\raisebox{1ex}{$-j{\beta }_{m}z$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$, and ${\widehat{a}}_{\mathit{zm}}\triangleq a_{\mathit{zm}}{e}^{\raisebox{1ex}{$-j{\beta }_{m}z$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ given by:

$$\frac{\partial {\widehat{a}}_{\mathit{xn}}}{\partial z}=-\frac{j}{2{\beta }_n}\sum _m\left[{\kappa }_{\left(\mathit{xx}\right)n,m}{\widehat{a}}_{\mathit{xm}}+{\kappa }_{\left(\mathit{xz}\right)n,m}^{\prime }\left(\frac{\partial {\widehat{a}}_{\mathit{zm}}}{\partial z}-j\frac{{\beta }_m}{2}{\widehat{a}}_{\mathit{zm}}\right)\right]e^{-\frac{j\left({\beta }_m-{\beta }_n\right)z}{2}}$$E62

$$\frac{\partial {\widehat{a}}_{\mathit{zn}}}{\partial z}=-\frac{j}{2{\beta }_n}\sum _m\left[{\kappa }_{\left(\mathit{zz}\right)n,m}{\widehat{a}}_{\mathit{zm}}+{\kappa }_{\left(\mathit{zx}\right)n,m}^{\prime }\left(\frac{\partial {\widehat{a}}_{\mathit{xm}}}{\partial z}-j\frac{{\beta }_m}{2}{\widehat{a}}_{\mathit{xm}}\right)\right]e^{-\frac{j\left({\beta }_m-{\beta }_n\right)z}{2}}$$E63

3.2 Two coupled waveguides

3.2.1 SH modes

The couple mode equations for two single-mode waveguides, illustrated in Figure 6, are given by Eqs. (64) and (65), using Eq. (47). The indices 1 and 2 refer to the two coupled waveguides.

Figure 6.

Two coupled waveguides. c represents the distance between waveguides. a and b are the widths of the core and cladding regions [38].

$$\frac{\partial a_{y1}\left(z\right)}{\partial z}=-j\frac{{\omega }^2}{2{\beta }_1}{\kappa }_{\left(\mathit{SH}\right)1,2}{a}_{y2}\left(z\right)e^{-j\left({\beta }_2-{\beta }_1\right)z}$$E64

$$\frac{\partial a_{y2}\left(z\right)}{\partial z}=-j\frac{{\omega }^2}{2{\beta }_2}{\kappa }_{\left(\mathit{SH}\right)2,1}{a}_{y1}\left(z\right)e^{-j\left({\beta }_1-{\beta }_2\right)z}$$E65

The CMEs are solved for two coupled waveguides, shown in Figure 6, considering both core width (a) and distance between waveguides (c) to be 200 μm and the cladding width (b) to be 4 mm to evaluate the coupling process. The result is illustrated in Figure 7. The materials of core and cladding regions are presumed copper and silicon, respectively. The coupling gives rise to the energy transfer between waveguides periodically. The value of elastic energy in each waveguide depends on the length of the coupled system [38].

Figure 7.

The amplitudes of waves in two coupled single-mode waveguides for the coupling of SH modes [38].

3.2.2 RL modes

The CMEs of Rayleigh-Lame modes for the coupled waveguides in Figure 6 are given by Eqs. (66–69) using Eqs. (62) and (63).

$$\frac{\partial {\widehat{a}}_{x1}}{\partial z}=-\frac{j}{2{\beta }_1}\left[{\kappa }_{\left(\mathit{xx}\right)1,2}{\widehat{a}}_{x2}+{\kappa }_{\left(\mathit{xz}\right)1,2}^{\prime }\left(\frac{\partial {\widehat{a}}_{z2}}{\partial z}-j\frac{{\beta }_2}{2}{\widehat{a}}_{z2}\right)\right]e^{-\frac{j\left({\beta }_2-{\beta }_1\right)z}{2}}$$E66

$$\frac{\partial {\widehat{a}}_{z1}}{\partial z}=-\frac{j}{2{\beta }_1}\left[{\kappa }_{\left(\mathit{zz}\right)1,2}{\widehat{a}}_{z2}+{\kappa }_{\left(\mathit{zx}\right)1,2}^{\prime }\left(\frac{\partial {\widehat{a}}_{x2}}{\partial z}-j\frac{{\beta }_2}{2}{\widehat{a}}_{x2}\right)\right]e^{-\frac{j\left({\beta }_2-{\beta }_1\right)z}{2}}$$E67

$$\frac{\partial {\widehat{a}}_{x2}}{\partial z}=-\frac{j}{2{\beta }_2}\left[{\kappa }_{\left(\mathit{xx}\right)2,1}{\widehat{a}}_{x1}+{\kappa }_{\left(\mathit{xz}\right)2,1}^{\prime }\left(\frac{\partial {\widehat{a}}_{z1}}{\partial z}-j\frac{{\beta }_1}{2}{\widehat{a}}_{z1}\right)\right]e^{-\frac{j\left({\beta }_1-{\beta }_2\right)z}{2}}$$E68

$$\frac{\partial {\widehat{a}}_{z2}}{\partial z}=-\frac{j}{2{\beta }_1}\left[{\kappa }_{\left(\mathit{zz}\right)2,1}{\widehat{a}}_{z1}+{\kappa }_{\left(\mathit{zx}\right)2,1}^{\prime }\left(\frac{\partial {\widehat{a}}_{x1}}{\partial z}-j\frac{{\beta }_1}{2}{\widehat{a}}_{x1}\right)\right]e^{-\frac{j\left({\beta }_1-{\beta }_2\right)z}{2}}$$E69

The displacement components in one waveguide are coupled to both displacement components of another one. The CMEs for the waveguide are solved with the dimensions mentioned in the previous section. The modal analysis should be first performed to calculate the coupling coefficients. Figure 8 illustrates the coupling results when the input displacement is considered in the x direction. Figure 8a and b depict the amplitudes in the x and z directions, respectively. On the contrary, in Figure 9, the input displacement is totally in the z-direction. Figure 9a and b are for displacement amplitude in the x and z directions, respectively.

Figure 8.

The coupling process for the first RL mode when ${\widehat{a}}_{x\mathit{1}}\left(\mathit{0}\right)=\mathit{1}$, ${\widehat{a}}_{z\mathit{1}}\left(\mathit{0}\right)=\mathit{0}$, ${\widehat{a}}_{x\mathit{2}}\left(\mathit{0}\right)=\mathit{0}$, and ${\widehat{a}}_{z\mathit{2}}\left(\mathit{0}\right)=\mathit{0}$. (a) and (b) are for amplitudes of displacements in the x and z directions, respectively [38].

Figure 9.

The coupling process for the first RL mode when ${\widehat{a}}_{x\mathit{1}}\left(\mathit{0}\right)=\mathit{0}$, ${\widehat{a}}_{z\mathit{1}}\left(\mathit{0}\right)=\mathit{1}$, ${\widehat{a}}_{x\mathit{2}}\left(\mathit{0}\right)=\mathit{0}$, and ${\widehat{a}}_{z\mathit{2}}\left(\mathit{0}\right)=\mathit{0}$. (a) and (b) are for amplitudes of displacements in the x and z directions, respectively [38].

3.2.3 SH wave coupling between non-parallel waveguides

The SH modes of a slab waveguide are given by Eq. (70), where t and C are the thickness of the core and normalization constant, respectively. The indices co and cl indicate the core and cladding regions. p and h are defined as Eqs. (71) and (72).

$$U_y\left(x\right)=\left\{\begin{array}{l}{\mathit{Ce}}^{-\mathit{px}},0\le x<+\infty \\ C\left[cos\left(\mathit{hx}\right)-\frac{{\mu }_{\mathit{cl}}p}{{\mu }_{\mathit{co}}h}sin\left(\mathit{hx}\right)\right],-t\le x\le 0\\ C\left[cos\left(\mathit{ht}\right)+\frac{{\mu }_{\mathit{cl}}p}{{\mu }_{\mathit{co}}h}sin\left(\mathit{ht}\right)\right]e^{p\left(x+t\right)},-\infty \le x<-t\end{array}\right.$$E70

$$p={\left({\beta }^2-n_{\mathit{cl}}^2{\omega }^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$E71

$$h={\left(n_{\mathit{co}}^2{\omega }^2-{\beta }^2\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$E72

$n_{\mathit{co}}$ and $n_{\mathit{cl}}$ are the refractive indices of the core and cladding, equaling the square root of ${\epsilon }^0$. The elastic power flux is expressed as Eq. (73), in which $\sigma$ is the stress tensor [39].

$$\iint -\left(\sigma .\frac{\partial \overrightarrow{u}}{\partial t}\right).\widehat{n}\mathit{dS}$$E73

The normalization constant is calculated by Eq. (74), assuming the power flux to be 1.

$$\left|C\right|=\frac{{\left({\mu }_{\mathit{co}}h\right)}^2+{\left({\mu }_{\mathit{cl}}p\right)}^2}{{\mu }_{\mathit{co}}\sqrt{2\mathit{\beta \omega }{\left({\mu }_{\mathit{cl}}p\right)}^3}}$$E74

The above-introduced modes are applied to the structure containing two non-parallel coupled waveguides illustrated in Figure 10 to investigate the coupling process. The displacement for SH waves can be written as Eq. (75) when the waveguides are not very close. The wave equation for the total displacement is given by Eq. (76), where $\mathrm{\Delta }{n}_{\mathit{co}1}^2\left(x_1\right)=n_{\mathit{co}1}^2\left(x_1\right)-n_{\mathit{cl}}^2$ and $\mathrm{\Delta }{n}_{\mathit{co}2}^2\left(x_2\right)=n_{\mathit{co}2}^2\left(x_2\right)-n_{\mathit{cl}}^2$ [40].

Figure 10.

Two non-parallel coupled waveguides [38].

$$u_y=a_{y1}\left(z_1\right)U_{y1}\left(x_1\right)e^{-\mathit{j\beta }{z}_1}+a_{y2}\left(z_2\right)U_{y2}\left(x_2\right)e^{-\mathit{j\beta }{z}_2}$$E75

$$\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial z^2}\right)u_y+{\omega }^2\left[n_{\mathit{cl}}^2+\mathrm{\Delta }{n}_{\mathit{co}1}^2\left(x_1\right)+\mathrm{\Delta }{n}_{\mathit{co}2}^2\left(x_2\right)\right]u_y=0$$E76

Eqs. (79–82) gives the CMEs for the structure depicted in Figure 10, considering the orthogonality of modes and slowly varying amplitude approximation. The coupling coefficients are calculated through Eqs. (79–82) [38].

$$\frac{\partial a_{y1}\left(z_1\right)}{\partial z_1}=-j\frac{{\omega }^2}{2\beta }\left[{\kappa }_{12}{a}_{y2}\left(z_2\right)e^{\mathit{j\beta }\left(1-cos\alpha \right)}+{\kappa }_{11}{a}_{y1}\left(z_1\right)\right]$$E77

$$\frac{\partial a_{y2}\left(z_2\right)}{\partial z_2}=-j\frac{{\omega }^2}{2\beta }\left[{\kappa }_{21}{a}_{y1}\left(z_1\right)e^{\mathit{j\beta }\left(1-cos\alpha \right)}+{\kappa }_{22}{a}_{y2}\left(z_2\right)\right]$$E78

$${\kappa }_{11}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y1}^{\ast }\left(x_1\right)\Delta n_{\mathit{co}2}^2\left(x_1,z_1\right)U_{y1}\left(x_1\right){\mathit{dx}}_1}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y1}^{\ast }\left(x_1\right)U_{y1}\left(x_1\right){\mathit{dx}}_1}$$E79

$${\kappa }_{12}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y1}^{\ast }\left(x_1\right)\Delta n_{\mathit{co}1}^2\left(x_1\right)U_{y2}\left(x_1,z_1\right)e^{-\mathit{j\beta }sin\alpha x_1}{\mathit{dx}}_1}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y1}^{\ast }\left(x_1\right)U_{y1}\left(x_1\right){\mathit{dx}}_1}$$E80

$${\kappa }_{22}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y2}^{\ast }\left(x_2\right)\Delta n_{\mathit{co}1}^2\left(x_2,z_2\right)U_{y2}\left(x_2\right){\mathit{dx}}_2}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y2}^{\ast }\left(x_2\right)U_{y2}\left(x_2\right){\mathit{dx}}_2}$$E81

$${\kappa }_{21}=\frac{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y2}^{\ast }\left(x_2\right)\Delta n_{\mathit{co}2}^2\left(x_2\right)U_{y1}\left(x_2,z_2\right)e^{-\mathit{j\beta }sin\alpha x_2}{\mathit{dx}}_2}{\underset{-\infty }{\overset{+\infty }{\int }}{U}_{y2}^{\ast }\left(x_2\right)U_{y2}\left(x_2\right){\mathit{dx}}_2}$$E82

For small angles between waveguides, the coupling coefficients are obtained by [38]:

$${\kappa }_{11}=\frac{\left(h^2+p^2\right)\left(e^{2\mathit{pt}}-1\right)}{{\omega }^2\left(2+\frac{2{\mu }_{\mathit{cl}}{p}^2}{{\mu }_{\mathit{co}}{h}^2}+\left(1+\frac{{\mu }_{\mathit{cl}}^2{p}^2}{{\mu }_{\mathit{co}}^2{h}^2}\right)\mathit{pt}\right)}{e}^{-2{\mathit{pz}}_{2}sin\alpha }={\kappa }_0{e}^{-2{\mathit{pz}}_{1}tan\alpha }$$E83

$${\kappa }_{12}=\frac{2p^2\left(\left(1+\frac{{\mu }_{\mathit{cl}}}{{\mu }_{\mathit{co}}}\right)\left(e^{\mathit{pt}}+1\right)-2\right)}{{\omega }^2\left(2+\frac{2{\mu }_{\mathit{cl}}{p}^2}{{\mu }_{\mathit{co}}{h}^2}+\left(1+\frac{{\mu }_{\mathit{cl}}^2{p}^2}{{\mu }_{\mathit{co}}^2{h}^2}\right)\mathit{pt}\right)}{e}^{-{\mathit{pz}}_{1}sin\alpha }={\kappa }_1{e}^{-{\mathit{pz}}_{1}sin\alpha }$$E84

$${\kappa }_{22}=\frac{\left(h^2+p^2\right)\left(e^{2\mathit{pt}}-1\right)}{{\omega }^2\left(2+\frac{2{\mu }_{\mathit{cl}}{p}^2}{{\mu }_{\mathit{co}}{h}^2}+\left(1+\frac{{\mu }_{\mathit{cl}}^2{p}^2}{{\mu }_{\mathit{co}}^2{h}^2}\right)\mathit{pt}\right)}{e}^{-2{\mathit{pz}}_{1}sin\alpha }={\kappa }_0{e}^{-2{\mathit{pz}}_{2}tan\alpha }$$E85

$${\kappa }_{21}=\frac{2p^2\left(\left(1+\frac{{\mu }_{\mathit{cl}}}{{\mu }_{\mathit{co}}}\right)\left(e^{\mathit{pt}}+1\right)-2\right)}{{\omega }^2\left(2+\frac{2{\mu }_{\mathit{cl}}{p}^2}{{\mu }_{\mathit{co}}{h}^2}+\left(1+\frac{{\mu }_{\mathit{cl}}^2{p}^2}{{\mu }_{\mathit{co}}^2{h}^2}\right)\mathit{pt}\right)}{e}^{-{\mathit{pz}}_{2}sin\alpha }={\kappa }_1{e}^{-{\mathit{pz}}_{2}sin\alpha }$$E86

The CMEs for the waveguides, whose core and cladding are similar to section 3.2.1, are solved, and the results are illustrated in Figure 11 for two different situations: the waveguides getting away from each other and vice versa. In the case that two waveguides move away from each other after a certain length due to the increase in the distance between the two waveguides, the interaction between them tends to zero, and the coupling stops. Consequently, the amplitude remains constant. Conversely, when two waveguides gradually get closer to each other, the intensity of coupling increases, and as a result, the elastic power is exchanged faster between two waveguides.

Figure 11.

(a) Two waveguides are moving away from each other. (b) Two waveguides are getting closer to each other [38].

3.3 Future challenges

Modal analysis of waveguides with inhomogeneity, such as impurities in materials or rough boundaries between the core and cladding regions, which cause mode coupling, can help to reach a deep understanding of the wave propagation in waveguides. Furthermore, investigation of the coupling in a system including more than two waveguides may be required to know how we can guide the desired value of elastic energy toward a specific point in optoelectronic or photonic devices. Moreover, exploring the quantum mechanical behavior of phonons, named by quantum acoustics, is a significant issue because of the growing demand for quantum computing; hence, designing devices that can generate, transport, and detect a single phonon is required.

4. Conclusion

Since phonons have an essential role in physics, designing suitable devices that allow physicists and engineers to control and manipulate them and their interaction with electrons and photons has been of interest to researchers. As nonlinearity is widely used to design desirable devices in electromagnetics, inspired by photonics and the similarity between photons and phonons, the nonlinear Schrödinger equation (NSE), which describes the nonlinear wave propagation in nonlinear and dispersive media, was developed for phononic waves. A phononic wave hard limiter, which confines the amplitude of the output wave to a desirable value, was introduced using NSE. The coupled mode theory was also developed for elastic slab waveguides, and the coupled mode equations are obtained for shear and Rayleigh-Lame modes. Finally, the shear wave coupling was investigated analytically for two non-parallel slab waveguides.