Efficient Analysis of Confined Guided Modes in Phoxonic Crystal Slabs

TL;DR

This paper introduces an improved plane wave expansion method for analyzing confined guided modes in phoxonic crystal slabs, overcoming convergence issues of standard methods and enabling efficient, accurate mode calculations for practical photonic and phononic devices.

Contribution

The authors propose a modified PWE method that uses eigenmodes of an unstructured slab for better convergence and accuracy in analyzing phoxonic crystal slabs.

Findings

The new method is faster and more accurate than conventional PWE.

It effectively estimates confined propagating modes in phoxonic slabs.

Application to a non-reciprocal photonic device demonstrates practical utility.

Abstract

Today's standard fabrication processes are just capable of manufacturing slab of photonic and phononic crystals, so an efficient method for analysis of these crystals is indispensable. Plane wave expansion (PWE) as a widely used method in studying photonic and phononic (phoxonic) crystals in full three dimensions is not suitable for slab analysis in its standard form, because of convergence and stability issues. Here, we propose a modification to this method which overcomes these limitations. This improved method can be utilized for calculation of both photonic and phononic modes in phoxonic slabs. While in the standard three-dimensional PWE, Fourier series are used to estimate the field dependence across the normal component of the slab, we expand the fields across the third dimension using eigenmodes of a plain unstructured slab. Despite its approximate nature, this approach is observed to be both much faster and more accurate than the conventional PWE method and can give a very accurate estimation of confined propagating modes. As an application example of the proposed method we investigate a non-reciprocal photonic device. This device is a phoxonic slab waveguide which changes modes of optical waves by elastic waves and can be used as an optical insulator or mode converter.

Figures

Tables

Table 1. TABLE I: Coupling Strength Between Optical Waves Through Elastic Waves due to Photoelastic Effect

	${\xi }_{S_1,S_3}$	${\xi }_{S_2,S_3}$
$T_1$	$0.9287+ı0.0875$	$0$
$T_2$	$-0.0046+ı0.7713$	$0$
$T_3$	$0$	$0.1011-ı0.2290$

Equations

$$(\lambda +\mu )\nabla \nabla \cdot \mathbf{u}+\mu {\nabla }^2\mathbf{u}=\rho \ddot{\mathbf{u}},$$

$$\sum _{j=1}^3\left[(\lambda +\mu )\frac{{\partial }^2{u}_j}{\partial x_j\partial x_i}+\mu \frac{{\partial }^2{u}_i}{\partial x_j^2}\right]=\rho \ddot{u_i},(i=1,2,3).$$

$$u_{y}(x,z)=\left\{\begin{array}[]{lr}A_{1}\cos(kz)e^{-\imath(qx-\Omega t)}&\text{Symmetric}\\ A_{2}\sin(kz)e^{-\imath(qx-\Omega t)}&\text{Anti-symmetric}\end{array}\right.,$$

$$k=n\pi /\text{t},n=0,1,2,\dots .$$

$$u_x=\left[A_2ıq\cos (k_{1}z)-B_1{k}_2\cos (k_{2}z)\right]e^{-ı(qx-\Omega t)},$$

$$u_z=\left[-A_2{k}_1\sin (k_{1}z)+B_1ıq\sin (k_{2}z)\right]e^{-ı(qx-\Omega t)},$$

$$\frac{\tan (k_2\text{t}/2)}{\tan (k_1\text{t}/2)}=-\frac{4k_1{k}_2{q}^2}{{\left(q^2-k_2^2\right)}^2},$$

$$\frac{A_2}{B_1}=-ı\frac{2q{k}_2\cos (k_2\text{t}/2)}{\left(q^2-k_2^2\right)\cos (k_1\text{t}/2)},$$

$$u_x=\left[A_1ıq\sin (k_{1}z)+B_2{k}_2\sin (k_{2}z)\right]e^{-ı(qx-\Omega t)},$$

$$u_z=\left[A_1{k}_1\cos (k_{1}z)+B_2ıq\cos (k_{2}z)\right]e^{-ı(qx-\Omega t)},$$

$$\frac{\tan (k_2\text{t}/2)}{\tan (k_1\text{t}/2)}=-\frac{{\left(q^2-k_2^2\right)}^2}{4k_1{k}_2{q}^2},$$

$$\frac{A_1}{B_2}=ı\frac{2q{k}_2\sin (k_2\text{t}/2)}{\left(q^2-k_2^2\right)\sin (k_1\text{t}/2)}.$$

$$\frac{\tan (k_2\text{t}/2)}{\tan (k_1\text{t}/2)}$$

$$k_1^2$$

$$k_2^2$$

$$\mathbf{u}(\mathbf{r})$$

$$=\sum _{\mathbf{G}}\left[u_{{\mathbf{G}}_x}(z)\hat{x}+u_{{\mathbf{G}}_y}(z)\hat{y}+u_{{\mathbf{G}}_z}(z)\hat{z}\right]e^{-ı(\mathbf{q}+\mathbf{G})\cdot {\mathbf{r}}_{xy}},$$

$$u_{{\mathbf{G}}_x}(z)$$

$$u_{{\mathbf{G}}_y}(z)$$

$$u_{{\mathbf{G}}_z}(z)$$

$$u_{{\mathbf{G}}_x}(z)$$

$$u_{{\mathbf{G}}_y}(z)$$

$$u_{{\mathbf{G}}_z}(z)$$

$$\sum _{j=1}^3$$

$$=\rho \ddot{u_i},(i=1,2,3).$$

$$\mathbb{L}|\mathbf{u}(\mathbf{r})⟩={\Omega }^2\mathbb{R}|\mathbf{u}(\mathbf{r})⟩,$$

$${\Omega }^2=\frac{⟨\mathbf{u}(\mathbf{r})|\mathbb{L}|\mathbf{u}(\mathbf{r})⟩}{⟨\mathbf{u}(\mathbf{r})|\mathbb{R}|\mathbf{u}(\mathbf{r})⟩},$$

$$\mathbb{M}{\left(\begin{array}{c}{A}_{{\mathbf{G}}_x}\\ B_{{\mathbf{G}}_x}\\ A_{{\mathbf{G}}_y}\\ B_{{\mathbf{G}}_y}\\ A_{{\mathbf{G}}_z}\\ B_{{\mathbf{G}}_z}\end{array}\right)}_{\mathbf{G}}+{\Omega }^2\mathbb{N}{\left(\begin{array}{c}{A}_{{\mathbf{G}}_x}\\ B_{{\mathbf{G}}_x}\\ A_{{\mathbf{G}}_y}\\ B_{{\mathbf{G}}_y}\\ A_{{\mathbf{G}}_z}\\ B_{{\mathbf{G}}_z}\end{array}\right)}_{\mathbf{G}}=0.$$

$$\mathbb{K}|\mathbf{H}(\mathbf{r})⟩=\frac{{\omega }^2}{c^2}|\mathbf{H}(\mathbf{r})⟩,$$

$$\mathbf{H}(\mathbf{r})$$

$$=\sum _{\mathbf{G}}\left[H_{{\mathbf{G}}_{\Vert }}(z){\mathbf{e}}_{{\mathbf{G}}_{\Vert }}+H_{{\mathbf{G}}_{\perp }}(z){\mathbf{e}}_{{\mathbf{G}}_{\perp }}+H_{{\mathbf{G}}_z}(z)\hat{z}\right]$$

$$\cdot e^{-ı(\mathbf{\kappa }+\mathbf{G})\cdot {\mathbf{r}}_{xy}},$$

$$H_{{\mathbf{G}}_{\Vert }}(z)$$

$$H_{{\mathbf{G}}_{\perp }}(z)$$

$$H_{{\mathbf{G}}_z}(z)$$

$$f_{{\mathbf{G}}_{\Vert }}(z)$$

$$f_{{\mathbf{G}}_{\perp }}(z)$$

$$f_{{\mathbf{G}}_z}(z)$$

$${\alpha }_{\mathbf{G}}$$

Full text

I Introduction

There exist many different types of waves in the nature, among them elastic and electromagnetic waves being the most famous and known ones. Because of the availability of sound waves as the scalar polarization of elastic pressure waves, their behavior was well known to physicists many years before the formulation of Maxwell’s equations. With the advent of electromagnetism and special relativity, the focus of study in wave phenomena shifted primarily to the realm of optics and electromagnetics and consequently photonic crystals (PtC) were introduced [1, 2] before their phononic counterparts [3]. They are called crystals since their constituents are periodic. Periodicity of permittivity in photonic crystals, and density and stiffness in phononic crystals (PnC) may be held responsible for their characteristics. The main property of these crystals is their ability to prevent some continuous frequency ranges of waves from propagating inside them. These frequency ranges are called forbidden photonic/phononic frequency gaps, which are caused by destructive interference of waves inside the crystal.

Because of this interesting property, photonic and phononic crystals have found many applications in different fields of science and technology [4], among them are optomechanics [5, 6, 7, 8, 9], sensors [10], acoustic and optical insulators, waveguides, filters and high $Q$ cavities [11, 12, 13, 14, 15, 16]. With the help of simultaneous phononic and photonic band gaps, one can create complete optical insulation [17]. Acoustic insulators can be used to isolate gyroscopes resonators and oscillators from ambient noise and vibrations[18]. By utilizing acoustic waveguides, one can miniaturize acoustic delay elements which are used in signal processing and oscillators based on delay line. These waveguides can also be used in impedance matching and focusing in acoustic photography like medical ultrasound [4].

By utilizing photonic crystals, power efficiency of solar cells [19, 20, 21, 22] and extraction efficiency of light emitting diodes [23, 24, 25] can be significantly improved. It is possible to embed a layer of PtC beneath miniaturized antennas manufactured on chips to increase the power of radiated waves through reflecting the waves, which would otherwise penetrate into the substrate [26, 27]. PtC cavities can be used to construct ultra-fast low-threshold tiny lasers [28, 29], which find use in the field of cavity quantum electrodynamics [30] as well. Like acoustic waveguides, optical waveguides created in photonic crystals are also widely used [31, 32, 33, 34].

Numerous methods have been developed for PnCs analysis, most of them being in common with those of PtCs. The most well-known methods are plane wave expansion (PWE), finite element (FE) and finite difference time domain (FDTD) [35]. Among these methods, PWE is more popular than the other two since it is much simpler to set up and suitable for periodic structures. It is based on Bloch theorem and approximating periodic quantities by their Fourier series. Its shortcoming arises when the structure is not periodic nor invariant in one or more dimensions or the wave does not propagate uniformly throughout the whole crystal. For instance, PnC slab in not periodic nor invariant along its thickness or surface acoustic waves (SAWs) propagate near the surface of a semi-infinite crystal. In order to utilize this method in these circumstances we have to make some changes either in its standard formulation or in our structure.

A reliable and conventional approach for dealing with PnC slabs is to sandwich the slab between two vacuum layers and create a 3D structure by stacking such sandwiches on top of each other[36]. The obtained structure is then analyzed with full 3D PWE method. Despite its accuracy, this method is time consuming and its convergence is normally non-uniform. Nevertheless, its accuracy is not desirable in solid-solid crystals of large elastic mismatch. To resolve this problem a new formulation has been developed based on inverse rule in contrast to conventional Laurent’s rule to obtain Fourier coefficients [37, 38].

Sigalas and Economoua tried to combine PWE with plate theory to calculate PnC slabs band structure [39]. Although it is a good approach for thin slabs, its error increases as the slab becomes thicker. Tanaka and Tamura proposed a modified PWE method for SAWs analysis [40]. They supposed surface waves to decay exponentially by getting away from the crystal surface. This idea was then utilized to study PnC slabs [41] since they are similar to semi-infinite crystals as they also have free surfaces. It is not an efficient method for crystal slab analysis because it needs a numerical root searching procedure which is time consuming and not straightforward at all.

There exists another variation of the PWE method for PnC slab analysis, proposed by Hsu and Wu, which is based on Mindlin’s plate theory [42]. It has been claimed that this method is much faster than full 3D PWE method but it is restricted to thin slabs and low frequencies. Much research and effort has been put into this domain to overcome the limitations of PWE while avoiding coding complexity of fintie element and boundary element formalisms [43, 44, 45, 46].

In this article we introduce a new modified PWE method which can be employed to find modes of both photonic and phononic modes of phoxonic crystal slabs which are composed of isotropic materials. It is shown to be much faster than the standard full 3D PWE method. Then we investigate optomechanical interactions in an optomechanical slab waveguide that is used to design a recently proposed non-reciprocal device [47] and can be as well utilized as an optical insulator.

II Modified PWE Method

PWE method is based on approximating wave profile by Fourier series. Although this approximation is inevitable for a real periodic crystal, it seems not to be so efficient for a crystal slab along its normal axis. In this section, we introduce functions which can replace Fourier series along the off-plane, or normal component to the slab. We furthermore suppose that the wave profile along the slab normal axis does not change significantly when we carve out some holes or insert some inclusions inside the slab, thus enabling the estimation of wave profile along the normal axis by that of a plain slab. Then, Fourier series is used to expand the waves along the two other orthogonal in-plane axes.

II-A Phononic Crystal Slab

To begin, we first have to find modes of a plain elastic plate. As stated above, the plate is assumed to be made of isotropic material. In this case wave equation is written as [48]

[TABLE]

which after expanding is equivalent to

[TABLE]

In this equation $x_{i}$ and $u_{i}$ are the $i$’th component of position and displacement respectively, $\lambda$ and $\mu$ are Lam$\acute{\text{e}}$ constants and $\rho$ is the mass density of the plate. It is simple to show that two kinds of waves can propagate inside this plate, longitudinal coupled with vertical transverse, and horizontal transverse. We assume the plate normal direction to be along the $z-$axis and the plate to be surrounded between $z=-\text{t}/2$ and $z=\text{t}/2$ surfaces. We also assume the wave to propagate along the $x-$axis. So we can write $\mathbf{u}(\mathbf{r})=\mathbf{u}(x,z)$. In the horizontal transverse wave, displacement is just along the $y-$axis and in the other one it is along $x-$ and $z-$axes. By solving wave equation, horizontal transverse wave is obtained as [48]

[TABLE]

in which $\Omega=\sqrt{(k^{2}+q^{2})c_{2}^{2}}$ is the angular frequency, $c_{2}=\sqrt{\mu/\rho}$ is the phase velocity of transverse wave and $q$ is the wave number. Applying boundary conditions, we conclude,

[TABLE]

The wave equation and boundary conditions give [48]

[TABLE]

for symmetric longitudinal and vertical transverse wave. The anti-symmetric wave becomes [48]

[TABLE]

The dispersion relation for longitudinal and vertical transverse wave is then obtained implicitly as

[TABLE]

where $c_{1}=\sqrt{(\lambda+2\mu)/\rho}$ is the phase velocity of longitudinal wave. Fig. 1 shows dispersion curves of both horizontal transverse and longitudinal coupled with vertical transverse waves for an isotropic Silicon plate.

Now, if we carve out a periodic pattern of holes inside this plate, according to Bloch theorem we may expect displacement wave to be of the form

[TABLE]

where $\tilde{\mathbf{u}}(\mathbf{r}_{xy},z)$ is a periodic function of $x$ and $y$ for every value of $-t/2<z<t/2$, $\bm{q}$ is the Bloch wave vector and $\mathbf{G}$ is a reciprocal lattice vector.

Two types of waves propagate inside a PnC slab: vertical displacement like (VD-like) and horizontal displacement like (HD-like). Displacement of the mid-plane of the slab is perpendicular to slab surface in VD-like and parallel to slab surface in HD-like wave. By adopting wave profile of longitudinal coupled with vertical transverse modes in a plain slab ((3a,b) for HD-like and (4a,b) for VD-like) we now approximate displacement components in VD-like modes as

[TABLE]

and in HD-like modes as

[TABLE]

In these equations we set $k_{1_{\mathbf{G}}}^{2}=\Omega(|\bm{q}+\mathbf{G}|)^{2}/c_{1}^{2}-|\bm{q}+\mathbf{G}|^{2}$ and $k_{2_{\mathbf{G}}}^{2}=\Omega(|\bm{q}+\mathbf{G}|)^{2}/c_{2}^{2}-|\bm{q}+\mathbf{G}|^{2}$, where $\Omega(\bm{q})$ of VD-like mode is obtained from the first anti-symmetric dispersion band and that of HD-like mode from the first symmetric dispersion band of a plain slab with the same material and thickness. Since $\lambda$ and $\mu$ are position dependent, (1) should be written as

[TABLE]

In a compact form it can be written as

[TABLE]

which is a generalized eigenvalue problem. We can find an upper bound for its eigenvalues by replacing $\mathbf{u}(\mathbf{r})$ in

[TABLE]

by its approximate value from (II-A)-(II-A) and searching for optimum values of $A_{\mathbf{G}}$ and $B_{\mathbf{G}}$ coefficients. It can be shown that the optimum value of these coefficients are obtained by solving the following eigenvalue problem,

[TABLE]

In this equation $\mathbb{M}$ and $\mathbb{N}$ are obtained by calculating inner products $\langle\mathbf{u}(\mathbf{r})|\mathbb{L}|\mathbf{u}(\mathbf{r})\rangle$ and $\langle\mathbf{u}(\mathbf{r})|\mathbb{R}|\mathbf{u}(\mathbf{r})\rangle$, respectively, and extracting $A_{\mathbf{G}}$ and $B_{\mathbf{G}}$ coefficients. The VD-like and HD-like band structure of a hexagonal lattice PnC slab calculated utilizing this method, which we call it variational from now on, is compared with the result of 3D PWE method in Fig. 2. It can be seen that both methods give similar bands for lower frequencies (below 20th band). This is a special crystal which, as is investigated in [49], has both photonic and phononic band gaps. Hence it can be used as a basis for phoxonic waveguide. The advantage of our method over 3D PWE method is its faster convergence. This fact is illustrated in Fig. 3. In this figure we have calculated the phase velocities of two modes of the crystal at two high symmetry points of its irreducible Brillouin zone, M*(5)* (fifth mode at M point) and K*(3)* (third mode at K point) with our variational and 3D PWE methods and plotted them versus their calculation times. We began with three plane waves along $x-$ and $y-$axes for both 3D PWE and our method, that is we set the bound of $\Sigma$ expression of Bloch wave to $N=1$. The number of plane waves along $z-$axis for PWE method is equal to $x-$axis at each point of the figure. We then increased the number of plane waves by 2, $(N=N+1)$, at each successive points of the plot to obtain more accurate results. The right most points of PWE method are related to $13$ plane waves $(N=6)$ along each axis while the right most points of our method are related to $19$ plane waves $(N=9)$ along $x-$ and $y-$ axes.

The variational method converges faster because the field trial function along $z-$ axis in our method needs just two unknown parameters to be determined while in the standard 3D PWE method more Fourier series terms (unknown parameters) are needed to obtain a trial function which its similarity to the exact eigen function is comparable with our method. In other words, our method uses a smarter trial function along $z-$ axis than 3D PWE method that with fewer unknown parameters can give an acceptable trial function. The fewer unknown parameters our problem has, the smaller the final matrix which its eigen value should be calculated would be and the faster it is prepared.

II-B Photonic Crystal Slab

In a similar manner, we can go through the same process to find a PtC slab modes. Here, the eigenvalue problem is

[TABLE]

in which $\mathbf{H}(\mathbf{r})$ is the magnetic field, $c$ is the light velocity in vacuum, $\omega$ is the angular frequency of electromagnetic wave and $\mathbb{K}=\nabla\times\left(\eta(\mathbf{r})\nabla\times(\cdot)\right)$, where $\eta(\mathbf{r})=1/\epsilon_{\text{r}}(\mathbf{r})$ is the relative impermeability of environment. The magnetic field inside and around a PtC slab can be written as

[TABLE]

where we have expanded its elements on the coordinate system of wave vectors. Here, this local coordinate system is represented by the unit vectors lying in the mid-plane of the slab and parallel with (perpendicular to) $\bm{\kappa}+\mathbf{G}$ as $\mathbf{e}_{\mathbf{G}_{\|}}$ ($\mathbf{e}_{\mathbf{G}_{\perp}}$). Adopting magnetic field profile around a dielectric slab waveguide, we choose our ansatz magnetic field for the transverse electric like (TE-like) mode as

[TABLE]

where $A_{\mathbf{G}_{\|}}$ and $A_{\mathbf{G}_{\perp}}$ are coefficients to be determined and

[TABLE]

In these equations $\alpha_{\mathbf{G}}=\sqrt{|\bm{\kappa}+\mathbf{G}|^{2}-\omega(|\bm{\kappa}+\mathbf{G}|)^{2}}$ and $\omega(|\bm{\kappa}+\mathbf{G}|)$ is estimated from the first TE dispersion band of a dielectric slab waveguide with the same material and thickness. This estimation is the best one for guided modes, i.e. the modes below the light cone. $k_{\mathbf{G}_{\|}}$ and $k_{\mathbf{G}_{\perp}}$ are also calculated from the following implicit relations.

[TABLE]

where $\epsilon_{r_{\text{eff}}}$ is the effective permittivity of the PtC slab. These functions are chosen in a manner to satisfy $\nabla\cdot\mathbf{H}(\mathbf{r})=0$. For transverse magnetic like (TM-like) modes, trial functions are similar to the above functions except that $f_{\mathbf{G}_{\|}}(z)$ and $f_{\mathbf{G}_{\perp}}(z)$ should be even functions and $f_{\mathbf{G}_{z}}(z)$ should be an odd function. Also, $k_{\mathbf{G}_{\|}}$ and $k_{\mathbf{G}_{\perp}}$ are calculated from

[TABLE]

In this mode $\omega(|\bm{\kappa}+\mathbf{G}|)$ is estimated from the first TM dispersion band of the dielectric slab waveguide mentioned above. To test this method, we calculated the band structure of the crystal slab analyzed in the previous sub-section with both 3D PWE and the proposed methods. Fig. 4a compares the results.

Since both of the methods are based on variational principle, each of them gives a lower frequency for a mode is more accurate in that mode. Inspecting Fig. 4a, reveals that our proposed method is more accurate in plenty of modes. In TE-like mode, our method results in bands which are very similar in shape to that of 3D PWE, while in TM-like mode, shapes of bands differ slightly. However, our method is more accurate in TM-like than TE-like modes.

In Fig. 4b we have compared the mode profiles calculated with variational and 3D PWE methods at two high symmetry points of Brillouin zone. It can be seen the results are very similar.

III Photo-elastic Interaction in Slab Optomechanical Waveguides

Interaction of electromagnetic and elastic waves originates from different mechanisms among them are electrostriction, magnetostriction, radiation pressure, piezoelectricity, photoelasticity, and interface displacement [50]. Electrostriction, magnetostriction and radiation pressure are mechanisms through which electromagnetic wave affects elastic wave. In piezoelectric materials electromagnetic and elastic waves can have mutual effects on each other through piezoelectricity phenomenon. The last two mechanisms, i.e. photoelasticity and interface displacement are responsible for effects of elastic wave on electromagnetic wave [51, 52]. In general, for piezo-electric nano-structured materials, radiation pressure terms which arise from moving boundaries [12, 53, 47] are negligible compared to the bulk electrostrictive forces arising from photo-elastic interactions. For non-piezoelectric media, radiation pressure and electrostriction have been shown to contribute comparably, or at least within the same order of magnitude [12]. A coupling coefficient between electromagnetic and elastic waves can be defined for each of these mechanisms. Different expressions have been presented for the coupling coefficient originated from photoelastic effect [51, 52, 47]. In this section we present a new expression for this coefficient which can give some insight about electromagnetic wave evolution in time in the presence of elastic wave. We use this expression to illustrate an application of our proposed method.

Now suppose we have an optomechanical slab waveguide in which both elastic and electromagnetic wave propagate along $x-$axis. We start with (10) but with a perturbed operator as

[TABLE]

in which $\Delta\mathbb{K}=\nabla\times\left(\Delta\bm{\eta}(\mathbf{r},t)\nabla\times(\cdot)\right)$ and

[TABLE]

In (26), $\Delta\bm{\eta}$ is the change in impermeability tensor caused by strain wave $\bm{\varepsilon}=\left[(\nabla\mathbf{u})^{T}+\nabla\mathbf{u}\right]/2$ and $\bm{P}$ is the photo-elastic tensor. With regards to (II-A) we can write $\Delta\bm{\eta}$ as

[TABLE]

where $\bm{v}(\mathbf{r})$ is a periodic tensor of $x$ and $q$ is the elastic wave number. Since $\mathbb{K}$ is a Hermitian operator, its eigenstates build a complete set on which we can expand $\mathbf{H}^{\prime}(\mathbf{r},t)$ as

[TABLE]

where the sum is over all wave numbers in the irreducible Brillouin zone and all eigenstates corresponding to each value of wave number. Here, we have discretized the wavenumber just for the sake of simplicity. It can be replaced by an integral as $\left(\sum_{\kappa}\longleftrightarrow\frac{a}{\pi}\int_{0}^{\pi/a}\,\mathrm{d}\kappa\right)$, where $a$ is the length of our crystal slab waveguide unit cell along $x$ axis. $\alpha_{n,\kappa}(t)$ are time-dependent constants to be determined. We now insert $\mathbf{H}^{\prime}(\mathbf{r},t)$ from (28) into (25) to obtain

[TABLE]

The second equation is obtained with the assumption, $\dot{\alpha}_{n,\kappa}(t)\ll\omega_{n,\kappa}$. By inner producting both sides of (III) by a specific eigenstate of $\mathbb{K}$ such as $\tilde{\mathbf{H}}_{m,\ell}(\mathbf{r}_{xy},z)e^{-\imath(\ell x-\omega_{m,\ell}t)}$ we get

[TABLE]

In this equation, $V$ is the volume in which we take the integral and should be enlarged to occupy the whole space and $l_{x}$ is the dimension of $V$ along the $x-$axis. We also have supposed the eigenstates to be normalized, i.e.

[TABLE]

where $\delta_{n,m}$ and $\delta_{\kappa,\ell}$ are Kronecker delta functions. By affecting $\Delta\mathbb{K}$, the equation becomes

[TABLE]

where

[TABLE]

is the coupling strength between $(\ell,m)$ and $(\kappa,n)$ modes through photoelastic effect. Replacing $\alpha_{m,\ell}(t)$ by $\beta_{m,\ell}(t)=\alpha_{m,\ell}(t)e^{\imath\omega_{m,\ell}t}$ and doing some simplifications, (III) becomes

[TABLE]

(III) is a system of ordinary differential equations. Suppose the system is initially in $(\kappa_{1},1)$ state, then the stronger the coupling between an arbitrary state and this initial state, the faster the growth of its coefficient in time. The first step in solving this system of equations is to calculate the coupling strength between different states. In the next subsection, we utilize our proposed method in previous sections for this purpose to demonstrate its usefulness in computing the optomechanical coupling of photonic and phononic modes [47].

III-A Application Example

To demonstrate one of the applications of our modified PWE method, we compute the coupling strength due to photoelasticity between different modes of a PtC slab waveguide caused by phononic modes. This phonomenon is at the heart of optomechanical interactions in phoxonic crystals and can be theoretically investigated by coupled-mode theory once the propagating eigenmodes of the slab are known [47]. Conventional PWE method is inappropriate here since it is too much time and memory consuming, and also suffers from non-uniform convergence. It means that the numerical solutions start to oscillate around the values ideally obtained at infinite number of expansion terms instead of exponential convergence [54]. At first, we need to design a simultaneous photonic and phononic (phoxonic) waveguide. For this purpose we use the hexagonal crystal slab analyzed in previous sections. As can be seen in Figs. 2 and 4, this crystal has both photonic and phononic forbidden gaps in its band structures, so it is a good candidate to be used as the base of our waveguide. By inserting a gap in the mentioned crystal along $\Gamma\text{K}$ direction, and following the design instructions given elsewhere [47], we obtain a phoxonic waveguide as is shown in the inset of Fig. 5a.

It can be proved that VD-like elastic waves cannot couple any two optical modes in a phoxonic slab waveguide through photoelasticity, i.e. the coupling strength between any two optical modes through a VD-like elastic wave equals zero. This is actually quite expected from overlap integrals of optomechanical interactions [47]. Hence, we calculated photonic TM-like and phononic HD-like band structures of the designed waveguide using our proposed variational method in this work. Figs. 5a and 5b show these band structures. To compute the optomechanical coupling strength we need to obtain $\Delta\bm{\eta}$. We used the photo-elastic tensor of Silicon reported in [55] for this purpose. The independent photo-elastic tensor elements of Silicon are $p_{1111}=-0.094$, $p_{1122}=0.017$ and $p_{1212}=-0.051$. The coupling strengths $\xi_{S_{1},S_{3}}$ and $\xi_{S_{2},S_{3}}$ through phononic modes $T_{1}$, $T_{2}$ and $T_{3}$ are gathered in Table I. These values are calculated after normalizing phononic waves profiles, i.e.

[TABLE]

The specified modes are marked in Figs. 5a and 5b.

IV Conclusion

A new approach for analysis of photonic and phononic crystal slab was proposed. This approach which is based on PWE method converges much faster than the conventional 3D PWE method and is uniformly convergent. Furthermore, it lets us analyze crystal slabs with large unit cells by a personal computer unlike 3D PWE method which needs a huge amount of memory in such cases. An application example was presented, where we constructed a non-reciprocal photonic device utilizing a phoxonic waveguide and calculated the coupling strength between different modes with the proposed method. Further details and deep discussions on this subject are being intended as a separate article to be published soon [56].