Polynomial expansion method for full wave three-dimensional analysis of dielectric waveguides and periodic structures

Kamalodin Arik; Mahmood Akbari

doi:10.1364/OE.519283

1. Introduction

The modal analysis of dielectric waveguides within photonic integrated circuits is a pivotal aspect of the design process, shedding light on potential modes of propagation. The quest for an accurate, efficient, and stable method to determine mode shapes and their propagation constants has led to a spectrum of computational techniques, broadly categorized into three groups. Firstly, fully numerical approaches including finite-difference (FD) [1,2], finite-element method (FEM) [3,4], and beam propagation method (BPM) [5,6] have been extensively reviewed [7–9]. While these methods yield accurate results, challenges emerge from inherent errors [8], numerical convergence issues, substantial memory requirements, and prolonged computation times, particularly as problem dimensions increase. Moreover, limitations arise due to solution space discretization or the necessity for absorbing (transparent) boundary conditions [11].

Secondly, approximate approaches, like the effective refractive index method [12], variational method [13], and coupled mode theory [14], offer ease of implementation and faster execution under specific assumptions. However, their precision and versatility may be compromised.

Lastly, semi-analytical methods grounded in mathematical treatments of optical waveguides have gained prominence. Techniques such as the transfer matrix method (TMM) [15,16], Fourier modal method (FMM) [17,18], rigorous coupled-wave analysis (RCWA) [19,20], and plane wave expansion method (PWEM) [21,22] have garnered attention. All these approaches or any other conventional modal method by Fourier expansion involve expanding the electromagnetic fields in a basis of plane waves. The reason for this is that plane waves are complete, can be systematically improved, and do not require any prior knowledge as to the distributions of the fields. However, it was shown that this technique suffers from poor plane-wave convergence [23]. In fact, the discontinuous nature of both the permittivity and the electromagnetic fields severely limits the accuracy of the conventional rigorous approaches which are based on modal expansion, particularly for high dielectric contrasts and higher frequencies. Fluctuation of the truncated series from the actual case is large and the convergence is very slow [23]. It should be noticed that the Fourier series of a periodic function converges spectrally fast with respect to the number of terms in the series, that is, with an algebraic order that increases with the number of available derivatives and exponentially fast for analytic functions. However, when the function in question is nonperiodic, the situation is very different. Regardless of how smooth this function is, convergence is slow and there is a permanent oscillatory overshoot close to the endpoints due to the Gibbs phenomenon [24]. Furthermore, in many cases, it is revealed that the slow convergence rate of the solution is not solely due to the Gibb’s phenomenon and the use of Fourier expansions, but can be also imputed to the slowly convergent eigensolutions and the inadequate formulation of the conventional eigenproblem [25]. The PWEM, as the dominant method in electromagnetics for calculating Bloch waves in periodic dielectric structures [22,26–28], encounters convergence challenges, especially with complex refractive index profiles [29]. For larger 2D or moderately sized 3D problems, its matrix equations exceed conventional memory capacities, as observed in recent studies [21]. Another alternative, the transmission line formulation (TLF) [30], encompasses a broader class of structures like 2D and 3D photonic crystals. This method expands unknown fields through spatial harmonics and models their propagation via voltage and current waves within a network of coupled transmission lines. However, the derived transcendental dispersion equation, akin to other modal methods, may exhibit hyperbolic or sinusoidal behavior near its roots, necessitating complex-function theory for root finding [31,32].

Expansion of functions in fast-converging series is of central importance in approximation theory and is, through the agency of spectral methods, fundamental in designing efficient computational methods for the numerical solution of differential equations [33]. Among which, the Legendre Spectral-Galerkin method occurs in a wide range of practical problems and applications in many fields of mathematics and physics [34] where the smooth and orthogonal known polynomials are particularly appealing owing to their superior properties: 1) they are orthogonal with respect to the unit weight function which makes them preferable in Galerkin methods for differential equations [35]. 2) they have excellent error properties in the approximation of a globally smooth function [35]. There has been a continuing interest in developing explicit bounds for the Legendre coefficients and error estimates for the Legendre series expansion of differentiable functions and many results can be found in the problem of the rate of convergence of Legendre approximation [34,36]. 3) A central dogma of approximation theory asserts that a continuous function on a bounded interval can be approximated arbitrarily closely by the polynomials and the smoother a continuous function, the faster its approximants converge [37]. One of the most remarkable advantages of Legendre approximations is that their accuracy depends solely upon the smoothness of the underlying functions [36]. 4) An explicit relationship between the expansion coefficients of a general derivative of an infinitely differentiable function in terms of those of the function greatly facilitates the setting up of an algebraic system to determine the unknown Legendre coefficients [38].

In our work, we build upon the Legendre polynomial expansion method initially introduced for the precise extraction of electromagnetic eigenmodes in stratified homogeneous waveguides [39]. Similar polynomial expansion has also been applied to analysis of diffraction gratings [40] and analysis of layers with an inhomogeneous refractive index profile [29]. This nonmodal method demonstrates behavior where the dispersion equation behaves as a polynomial, enabling accurate computation of modes with closely aligned propagation constants [41]. Crucially, this approach alleviates the need for an ordinary transcendental dispersion equation of conventional modal expansion methods, which often grapple with numerical overflow and ill-conditioned matrices due to the simultaneous presence of extremely large and small coefficients arising from evanescent modes [42,43], whereas in this method an easy to manipulate algebraic dispersion equation can be obtained without any problem of numerical instability [39,40] not only because the involved matrices are very sparse, but also because the numbers constituting the matrices are neither extremely small nor extremely large, i.e. the condition numbers of the matrices are very good [39], which can be originated from the decay rates properties of the coefficients in the Legendre series expansion. This approach works properly even in those special cases in which other methods usually fail to render numerically stable results [42]. The physical intuition behind the accuracy of the method can be described by the fact that, even though in practice the expansion of the electromagnetic field in Hilbert space spanned by Legendre polynomials is truncated, each polynomial that remains in the calculation contains some information from all the infinite number of the natural modes of the system through their projection. Thus, no modal information on the system is lost in the truncation process [39].

In this manuscript, we present an extension of this nonmodal polynomial expansion-based approach tailored for the analysis of general multilayer photonic crystals comprising homogeneous dielectric layers and one or more biperiodic permittivity layers. Commencing with an equivalent matrix representation of Maxwell’s equations, we derive a governing mathematical relation describing the spatial evolution of coupled space harmonics within the structure, accounting for mixed polarizations. Utilizing Legendre polynomial basis functions, we expand the amplitudes of these space harmonics, transforming a set of linear differential equations into a more manageable algebraic system. Subsequently, we apply appropriate boundary conditions and solve for the unknown expansion coefficients and propagation constants. The resulting field distributions are then readily obtained as linear combinations of the expansion functions in accordance with the Floquet theorem, with the unknown expansion coefficients serving as weights. The proposed method is efficient in many commonly encountered problems and can be easily implemented for analyzing different problems including photonic-crystal band-structure calculation, and guided-mode analysis of different waveguides.

Notably, for 3D nonperiodic waveguides, our method can be applied to a periodically repeated version of the structure, allowing for analysis under analogous circumstances. This equivalence enables exploration of the artificial periodic structure's propagation characteristics in relation to the original configuration.

The structure of this paper is organized as follows: Section 2 elaborates on the Legendre polynomial expansion formulation for analyzing stratified biperiodic waveguides. Section 3 presents three numerical examples showcasing the applicability and effectiveness of our proposed method. Finally, Section 4 encapsulates our conclusions.

2. Formulation

2.1. Maxwell's equations within a biperiodic layer

Figure 1 shows a general planar waveguide composed of a finite number of stratified dielectric layers (say $NL)$ with biperiodic permittivity. We assume that the multilayer stack under investigation is made of nonmagnetic, linear and isotropic dielectric materials. The coordinate system is oriented such that the z-axis is normal to the layers of the structure. Each layer indexed by $l$ (which varies from zero to $NL + 1)$ with thickness ${d_l}$ and extending from $z = {z_{l - 1}}$ to $z = {z_{l - 1}} + {d_l},$ is assumed piecewise invariant in the z-direction. The semi-infinite layers at the top and the bottom of the structure are called cover and substrate, respectively. Formulation of a rigorous semi-analytical method begins with Maxwell’s equations describing the fields inside a single layer of the device. Considering the $l$th layer of the structure, its relative permittivity ${\varepsilon _r}(x,\,y,z)$ is biperiodic in x and y directions with periods of ${\Lambda _x}$ and ${\Lambda _y},$ respectively. It can thus be uniquely expanded using the following two-dimensional Fourier series, i.e., by:

(1a)$${\varepsilon _r}(x,\,y,\,z)\, = \sum\limits_{m ={-} \infty }^\infty {\sum\limits_{n ={-} \infty }^\infty {{\varepsilon _{mn}}(z)\,{e^{ - j(\frac{{2m\pi }}{{{\Lambda _x}}}x + \frac{{2n\pi }}{{{\Lambda _y}}}y)}}} }$$

where the integers m and n in the expansion are the indices of planar gratings comprising the permittivity throughout the unitcell. The z-dependent complex amplitudes ${\varepsilon _{mn}}(z)$ for the $mn$th planar grating are given by:

(1b)$${\varepsilon _{mn}}(z) = \frac{1}{{{\Lambda _x}{\Lambda _y}}}\int_0^{{\Lambda _x}} {\int_0^{{\Lambda _y}} {{\varepsilon _r}(x,\,y,\,z)\,} } {e^{j(\frac{{2m\pi }}{{{\Lambda _x}}}x + \frac{{2n\pi }}{{{\Lambda _y}}}y)}}dxdy$$

Fig. 1. The 3D view of a layered waveguide composed of a stack of layers with biperodic permittivity.

Download Full Size | PDF

This can be accomplished either numerically by building arrays that show the material properties as a function of position and then calculating Fourier coefficients through a fast Fourier transform (FFT) of the arrays or analytically as closed form expressions. Note that a homogeneous layer (a layer composed of a single dielectric) of the structure is a special case of a biperiodic layer.

As the goal of analysis, the propagating modes of the structure whose propagation vectors are parallel to the xy-plane are to be determined. Regarding the periodicity of the structure in the xy-plane, the Bloch theorem requires that the source-free solutions to Maxwell’s equations satisfy the following conditions:

(2a)$$\vec{E}(\vec{r} + \vec{\Lambda }) = \vec{E}(\vec{r})\,{e^{ - j\vec{k}.\vec{\Lambda }}}$$

(2b)$$\vec{H}(\vec{r} + \vec{\Lambda }) = \vec{H}(\vec{r})\,{e^{ - j\vec{k}.\vec{\Lambda }}}$$

which should hold at every position vector $\vec{r} = x\hat{x} + y\hat{y} + z\hat{z}.$ In Eq. (2), $\vec{k} = {k_x}\hat{x} + {k_y}\hat{y}$ denotes the still unknown wavevector for a propagation mode parallel to the xy-plane and $\vec{\Lambda } = {\Lambda _x}\hat{x} + {\Lambda _y}\hat{y}$ is a vector connecting diagonal corners of the unitcell of a layer of the structure. It is straightforward to prove that the following pseudo-Fourier series satisfy Eq. (2):

(3a)$$\vec{E}(\vec{r}) = \sum\limits_{m ={-} \infty }^\infty {\sum\limits_{n ={-} \infty }^\infty {{{\vec{E}}_{mn}}(z)\,{e^{ - j({\alpha _{mn}}x + {\beta _{mn}}y)}}} }$$

(3b)$$\vec{H}(\vec{r}) = \sum\limits_{m ={-} \infty }^\infty {\sum\limits_{n ={-} \infty }^\infty {{{\vec{H}}_{mn}}(z)\,{e^{ - j({\alpha _{mn}}x + {\beta _{mn}}y)}}} }$$

where ${\alpha _{mn}} = {{2m\pi } / {{\Lambda _x} + {k_x}}}$ and ${\beta _{mn}} = {{2n\pi } / {{\Lambda _y} + {k_y}}}.$ In Eq. (3), the integers $m$ and n are the indices of the space harmonics in the expansion. ${\vec{E}_{mn}}(z)\,$ and ${\vec{H}_{mn}}(z)$ are the complex amplitudes of the electric and magnetic fields, respectively, for the $mn$th space harmonic. In order to uniquely characterize a propagation mode, we must determine the space harmonics amplitudes ${\vec{E}_{mn}}(z)\,$ and ${\vec{H}_{mn}}(z)$ along with ${k_x}$ and ${k_y}.$ To be solved on a computer, the summations in Eq. (3) is inevitably truncated to a finite set of space harmonics. Sorting the electric and magnetic field components from $m ={-} M$ to $m = M$ and $n ={-} N$ to $n = N,$ we define $K \times 1$ matrices $[{{E_u}(z)} ]$ and $[{{H_u}(z)} ]$ the entries of which are ${E_{u,mn}}(z)$ and ${H_{u,mn}}(z)$, respectively, with $\textrm{u} \in \{{\textrm{x,}\,\textrm{y,}\,\textrm{z}} \}$ and $K = (2M + 1) \times (2N + 1).$ Now, it is necessary to substitute Eq. (3) into Maxwell’s equations. The pseudo-Fourier series of the field components implies that their spatial derivatives with respect to x (or y) are equivalent to multiplication of their corresponding matrix representation by diagonal matrix $- j[\alpha ]$ (or $- j[\beta ])$ having $- j{\alpha _{mn}}$ (or $- j{\beta _{mn}})$ as its diagonal elements. Furthermore, we can rewrite the constitutive relation $\vec{D}(\vec{r}) = {\varepsilon _0}{\varepsilon _r}(\vec{r})\vec{E}(\vec{r})$ using the expansion coefficients ${\varepsilon _{mn}}(z)$ and ${\vec{E}_{mn}}(z)\,$ as:

(4a)$$\vec{D}(\vec{r}) = {\varepsilon _0}\sum\limits_{m ={-} \infty }^\infty {\sum\limits_{n ={-} \infty }^\infty {{\varepsilon _{mn}}(z) \ast {{\vec{E}}_{mn}}(z)\,{e^{ - j({\alpha _{mn}}x + {\beta _{mn}}y)}}} }$$

in which ${\ast} $ signifies two-dimensional convolution, i.e.

(4b)$${\varepsilon _{mn}}(z) \ast {\vec{E}_{mn}}(z) = \sum\limits_{k ={-} \infty }^\infty {\sum\limits_{l ={-} \infty }^\infty {{\varepsilon _{(m - k)(n - l)}}{{\vec{E}}_{kl}}(z)} }$$

After some algebraic manipulations we arrive at $[{{D_u}(z)} ]= {\varepsilon _0}[{\varepsilon _r}][{{E_u}(z)} ]$ with $u \in \{{x,\,y,\,z} \}$ where the convolution matrix $[{\varepsilon _r}]$ which encompasses the Fourier coefficients is a $\textrm{K} \times \textrm{K}$ square matrix and constructed according to the summations in Eq. (4b). Moreover, the constitutive relation $\vec{B}(\vec{r}) = {\mu _0}\vec{H}(\vec{r})$ converts into matrix form $[{{B_u}(z)} ]= {\mu _0}[{{H_u}(z)} ]$ with $u \in \{{x,\,y,\,z} \}$ as well. Substituting the above matrix representation of the fields and their derivatives along with the constitutive relations in Maxwell’s equations and eliminating the normal components of electric and magnetic fields ($[{{E_z}(z)} ]$ and $[{{H_z}(z)} ])$ will lead to the following equations:

(5a)$$\frac{d}{{dz}}\left( \begin{array}{l} [{{E_x}} ]\\ {[{{E_y}} ]}\end{array} \right) ={-} j\omega [L ]\left( \begin{array}{l} [{{H_y}} ]\\ - [{{H_x}} ]\end{array} \right)$$

(5b)$$\frac{d}{{dz}}\left( \begin{array}{l} [{{H_y}} ]\\ - [{{H_x}} ]\end{array} \right) ={-} j\omega [C ]\left( \begin{array}{l} [{{E_x}} ]\\ {[{{E_y}} ]}\end{array} \right)$$

in which the $2K \times 2K$ matrices $[L ]$ and $[C ]$ are given by [31]:

(6a)$$[L ]= {\mu _0}\left[ \begin{array}{ll} I - [{\bar{\alpha }} ]{[{{\varepsilon_r}} ]^{ - 1}}[{\bar{\alpha }} ]& - [{\bar{\alpha }} ]{[{{\varepsilon_r}} ]^{ - 1}}[{\bar{\beta }} ]\\ \, - [{\bar{\beta }} ]{[{{\varepsilon_r}} ]^{ - 1}}[{\bar{\alpha }} ]&I - [{\bar{\beta }} ]{[{{\varepsilon_r}} ]^{ - 1}}[{\bar{\beta }} ]\end{array} \right]$$

(6b)$$[C ]= {\varepsilon _0}\left[ \begin{array}{lc} [{{\varepsilon_r}} ]- \,{[{\bar{\beta }} ]^2}&[{\bar{\beta }} ][{\bar{\alpha }} ]\\ \,[{\bar{\alpha }} ][{\bar{\beta }} ]&[{{\varepsilon_r}} ]- {[{\bar{\alpha }} ]^2}\, \end{array} \right]$$

in Eq. (6), $[{\bar{\alpha }} ]= {{[\alpha ]} / {{k_0}}},[{\bar{\beta }} ]= {{[\beta ]} / {{k_0}}},$where ${k_0} = \omega \sqrt {{\mu _0}{\varepsilon _0}}$ and I is a $\textrm{K} \times \textrm{K}$ identity matrix. By assuming $[V(z)] = {[{[{{E_x}} ],\,[{{E_y}} ]} ]^T}$ and $[I(z)] = {[{[{{H_y}} ],\, - [{{H_x}} ]} ]^T},$it can be easily shown that the system of Eq. (5) may be reduced to Helmholtz equation as:

(7a)$$\frac{{{d^2}}}{{d{z^2}}}[{V(z)} ]={-} {\omega ^2}[L ][C ][{V(z)} ]$$

in $l$th layer.

2.2. Polynomial expansion of electromagnetic fields inside the biperiodic layer

For any biperiodic permittivity profile, the resulting equation in Eq. (7a) governs the electromagnetic fields that can propagate in the layer and can be reformulated as a set of linear and second-order coupled differential equations as:

(7b)$$\frac{{{d^2}}}{{d{z^2}}}{v_i}(z) = \sum\limits_{j = 1}^{2K} {{k_{i,\,j}}{v_j}(z)} ,\quad i = 1,\,2,\,\ldots ,2K$$

where $[V(z)]$ containing a set of space harmonic amplitudes $\{{{v_i}(z)} \}_{i = 1}^{2K}$ is determined by solving this system of equations, subject to appropriate continuity conditions at the interfaces and ${k_{i,\,j}}$ is the $(i,\,j)$ element of the square matrix $- {\omega ^2}[L ][C ]$ denoted by $[k]$ for simplicity. In the presented method, each space harmonic amplitude is now expanded in a complete space spanned by Legendre polynomials:

(8)$${v_i}(z) = \sum\limits_{n = 0}^{{M_0}} {q_n^i{P_n}(\xi )} ,\quad \xi = \frac{{2z - {z_{l - 1}} - {z_l}}}{{{z_l} - {z_{l - 1}}}},\quad {z_{l - 1}} < z < {z_l}$$

where ${M_0}$ is the number of polynomials retained in the expansion and determined by the required level of accuracy. Note that each of the polynomials remaining in the calculation contains the projection of all electromagnetic eigenmodes of the system. Thus, no modal information of the system is fully lost in the truncation process. ${P_n}(\xi )$ stand for the $n$th order normalized Legendre polynomial and $q_n^i$ are the corresponding unknown expansion coefficients for $i$th space harmonic. After substituting the expansion of Eq. (8) in Eq. (7b) and using normalized Legendre polynomial properties, one obtains:

(9a)$$\sum\limits_{n = 0}^{{M_0} - 2} {r_n^i{P_n}(\xi ) = {{(\frac{{{d_l}}}{2})}^2}\sum\limits_{j = 1}^{2K} {{k_{i,\,j}}\sum\limits_{n = 0}^{{M_0}} {q_n^j{P_n}(\xi )} } } ,\quad i = 1,\,2,\,\ldots ,\,2K\,$$

in Eq. (9a), $r_n^i$ are the Legendre coefficients of the 2nd derivative of ${v_i}(z)$ and can be analytically carried out as:

(9b)$$r_n^i = \left( {n + \frac{1}{2}} \right)\sum\limits_{p = n + 2,\,n + 4,\,\ldots }^{{M_0}} {(p + n + 1)(p - n)q_p^i}$$

Projecting the preceding set of equations on orthogonal Legendre polynomial basis functions, i.e., multiplying Eq. (9a) by ${P_m}(\xi ),\,m = 0,\,1,\,\ldots ,\,{M_0} - 2$ and integrating over the interval $[{ - 1,\,1} ]$ yields the following set of coupled algebraic equations for the $l$th layer:

(10a)$$\left\{ \begin{array}{l} r_m^1 = {(\frac{{{d_l}}}{2})^2}[{{k_{1,1}}q_m^1 + {k_{1,2}}q_m^2 + \ldots + {k_{1,\,2K}}q_m^{2K}} ]\\ r_m^2 = {(\frac{{{d_l}}}{2})^2}[{{k_{2,1}}q_m^1 + {k_{2,2}}q_m^2 + \ldots + {k_{2,\,2K}}q_m^{2K}} ]\\ {\kern1.2cm} .{\kern1.2cm} .\,{\kern1.2cm} .\,\,\,{\kern1.2cm} .\\ {\kern1.2cm} .{\kern1.2cm} .\,{\kern1.2cm} .\,\,\,{\kern1.2cm} .\,\,\,\\ r_m^{2K} = {(\frac{{{d_l}}}{2})^2}[{{k_{2K,1}}q_m^1 + {k_{2K,2}}q_m^2 + \ldots + {k_{2K,\,2K}}q_m^{2K}} ]\end{array}\quad,\,m = 0,\,1,\,\ldots\,{M_0} - 2 \right.$$

Some straightforward algebraic manipulations can be done to recast the constituent factors of the both left and right hand side of Eq. (10a) and have them converted into an equivalent set of linear equations, which can be expressed in the following matrix form:

(10b)$$[{A_L}][q ]= [{A_R}][q ]$$

where $[q ]$ is the column vector of $(({M_0} + 1) \times 2K)$ unknown coefficients with the following entries ${[{q_0^1,q_0^2,\ldots ,\,q_0^{2K},\,q_1^1,q_1^2,\ldots ,\,q_1^{2K},\,\ldots ,\,q_{{M_0}}^1,q_{{M_0}}^2,\ldots ,\,q_{{M_0}}^{2K}} ]^T}.$ Accordingly, the matrix $[{A_R}]$ has $(({M_0} - 1) \times 2K)$ rows and $(({M_0} + 1) \times 2K)$ columns and can be subdivided into $({M_0} - 1)$ submatrices along the main diagonal:

(10c)$$[{A_R}] = {\left( {\frac{{{d_l}}}{2}} \right)^2}\left[ \begin{array}{ccccccc} [k ]&0&0&.&.&0\\ 0&[k ]&0&.&.&0\\ 0&0&.&0&.&0&\\ .&.&.&.&0&.\\ .&.&.&0&[k ]&0\\ 0&0&0&.&0&[k ]\end{array} \right]\,$$

and the matrix $[{A_L}]$ is numerically constructed based on the summation in Eq. (9b). Finally, Eq. (10b) can be reformulated as:

(10d)$${[A]_{(({M_0} - 1) \times 2K) \times (({M_0} + 1) \times 2K)}}{[q ]_{(({M_0} + 1) \times 2K) \times 1}} = 0$$

where, $[A] = [{A_R}] - [{A_L}].$ However, this leaves us with $4K$ missing equations to fully determine the unknown $q_m^i$ coefficients., which are to be found by applying the appropriate boundary conditions in next section.

2.3. Applying boundary conditions

For the $l$th layer sandwiched between two $(l - 1)$th and $(l + 1)$th layers, boundary conditions call for the continuity of the tangential electromagnetic fields at $z = {z_{l - 1}}$ and $z = {z_l}$ as:

(11a)$$[{{V_{l - 1}}({z_{l - 1}})} ]= [{{V_l}({z_{l - 1}})} ]$$

(11b)$$[{{I_{l - 1}}({z_{l - 1}})} ]= [{{I_l}({z_{l - 1}})} ]$$

(11c)$$[{{V_l}({z_l})} ]= [{{V_{l + 1}}({z_l})} ]$$

(11d)$$[{{I_l}({z_l})} ]= [{{I_{l + 1}}({z_l})} ]$$

It should be noted that to compute $[{V(z)} ]$ (or $[{I(z)} ])$ at $z = {z_{l - 1}}$ and $z = {z_l},$ one needs normalized Legendre polynomials (or their first derivatives) at the ends of the orthogonality interval, i.e. $\xi ={-} 1$ and $\xi ={+} 1,$ respectively, which can be taken into account as:

(12)$${P_m}({\pm} 1) = {({\pm} 1)^m},\quad\,{\left. {\frac{{d{P_m}(\xi )}}{{d\xi }}} \right|_{\xi ={\pm} 1}} = {({\pm} 1)^{m + 1}}\frac{{m(m + 1)}}{2}$$

Now, the continuity equations in Eq. (11) are projected on Legendre basis polynomials. Eliminating $[{I(z)} ]$ from Eq. (11b) and Eq. (11d) by means of Eq. (5a), and substituting for space harmonic amplitudes, four set of algebraic equations are obtained as:

(13a)$$\begin{array}{l} {\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l - 1)}}{P_n}( - 1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l - 1)}}{P_n}( - 1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l - 1)}}{P_n}( - 1)} } \right]^T} = \\ {\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l)}}{P_n}( - 1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l)}}{P_n}( - 1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l)}}{P_n}( - 1)} } \right]^T} \end{array}$$

(13b)$$\begin{array}{l} \frac{{ - 2}}{{j\omega }}\frac{{{{[{{L_{l - 1}}} ]}^{ - 1}}}}{{{d_{l - 1}}}}{\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l - 1)}}{P_n}^\prime ( - 1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l - 1)}}{P_n}^\prime ( - 1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l - 1)}}{P_n}^\prime ( - 1)} } \right]^T} = \\ \frac{{ - 2}}{{j\omega }}\frac{{{{[{{L_l}} ]}^{ - 1}}}}{{{d_l}}}{\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l)}}{P_n}^\prime ( - 1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l)}}{P_n}^\prime ( - 1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l)}}{P_n}^\prime ( - 1)} } \right]^T} \end{array}$$

(13c)$$\begin{array}{l} {\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l)}}{P_n}(1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l)}}{P_n}(1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l)}}{P_n}(1)} } \right]^T} = \\ {\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l + 1)}}{P_n}(1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l + 1)}}{P_n}(1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l + 1)}}{P_n}(1)} } \right]^T} \end{array}$$

(13d)$$\begin{array}{l} \frac{{ - 2}}{{j\omega }}\frac{{{{[{{L_l}} ]}^{ - 1}}}}{{{d_l}}}{\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l)}}{P_n}^\prime (1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l)}}{P_n}^\prime (1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l)}}{P_n}^\prime (1)} } \right]^T} = \\ \frac{{ - 2}}{{j\omega }}\frac{{{{[{{L_{l + 1}}} ]}^{ - 1}}}}{{{d_{l + 1}}}}{\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(l + 1)}}{P_n}^\prime (1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(l + 1)}}{P_n}^\prime (1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(l + 1)}}{P_n}^\prime (1)} } \right]^T} \end{array}$$

Substitution of Eq. (13a) in Eq. (13b) results in:

(14a)$${[{B_1}]_{2K \times (({M_0} + 1) \times 2K)}}{[q ]_{(({M_0} + 1) \times 2K) \times 1}} = 0$$

and similarly, substitution of Eq. (13c) in Eq. (13d) renders:

(14b)$${[{B_2}]_{2K \times (({M_0} + 1) \times 2K)}}{[q ]_{(({M_0} + 1) \times 2K) \times 1}} = 0$$

where, the matrices $[{B_1}]$ and $[{B_2}]$ can be numerically calculated. Equation (14) provide the missing equations and along with Eq. (10d) form a complete set of equations, through which the unknown coefficients $q_m^i$ are obtained. However, boundary conditions at the bottom and the top of a layered structure differ from those in Eq. (11). Simulating the open nature of the cross section, we apply transparent boundary conditions at $z = {z_0}$ and $z = {z_N}$ as:

(15a)$$[{{V_1}({z_0})} ]= [{{Z_s}} ][{{I_1}({z_0})} ]$$

(15b)$$[{{V_N}({z_N})} ]= [{{Z_c}} ][{{I_N}({z_N})} ]$$

where, $[{{Z_i}} ]= \sqrt {{{{L_i}} / {{C_i}}}}$ with $i \in \{{s,\,c} \}$ corresponding to substrate and cover are numerically obtained. The preceding procedure for Eq. (11) is employed to project the resultant equations in Eq. (15) into the Hilbert space spanned by Legendre polynomials as:

(16a)$$\begin{array}{l} {\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(1)}}{P_n}( - 1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(1)}}{P_n}( - 1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(1)}}{P_n}( - 1)} } \right]^T} = [{{Z_s}} ]\times \\ \frac{{ - 2}}{{j\omega }}\frac{{{{[{{L_1}} ]}^{ - 1}}}}{{{d_1}}}{\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(1)}}{P_n}^\prime ( - 1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(1)}}{P_n}^\prime ( - 1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(1)}}{P_n}^\prime ( - 1)} } \right]^T} \end{array}$$

(16b)$$\begin{array}{l} {\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(N)}}{P_n}(1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(N)}}{P_n}(1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(N)}}{P_n}(1)} } \right]^T} = [{{Z_c}} ]\times \\ \frac{{ - 2}}{{j\omega }}\frac{{{{[{{L_N}} ]}^{ - 1}}}}{{{d_N}}}{\left[ {\sum\limits_{n = 0}^{{M_0}} {{{({q_n^1} )}_{(N)}}{P_n}^\prime (1),\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^2} )}_{(N)}}{P_n}^\prime (1)} } ,\cdots ,\,\sum\limits_{n = 0}^{{M_0}} {{{({q_n^{2K}} )}_{(N)}}{P_n}^\prime (1)} } \right]^T} \end{array}$$

which can be, respectively, arranged in the following matrix forms as,

(17a)$${[{B_s}]_{2K \times (({M_0} + 1) \times 2K)}}{[q ]_{(({M_0} + 1) \times 2K) \times 1}} = 0$$

(17b)$${[{B_c}]_{2K \times (({M_0} + 1) \times 2K)}}{[q ]_{(({M_0} + 1) \times 2K) \times 1}} = 0$$

where the matrices $[{B_s}]$ and $[{B_c}]$ can be numerically constructed. At last, we add Eq. (17) to the previous equations that result from the Helmholtz equation and the boundary conditions to balance the number of unknown coefficients and independent equations. In this fashion, a set of $(NL \times ({M_0} + 1) \times 2K)$ homogeneous equations is derived. Evidently for nontrivial solutions, we set the determinant of the coefficient matrix of the resulting set of unforced equations to zero for a given $\vec{k}$ at a certain angular frequency $\omega .$ Consequently, the dispersion diagram of the multilayer stack under investigation and as well as the corresponding electromagnetic fields are obtained.

3. Numerical results

In this section, numerically accurate results for two commonly used photonic structures will be illustrated to assess the capability of the proposed method.

As the first numerical example, consider the structure consisting of a single row of parallel and infinitely long dielectric rods of generally rectangular cross section in air, described in [44] and shown in the inset of Fig. 2. The lattice constant $a = 0.1\,\mathrm{\mu}\textrm{m}$ denotes the spacing between the rods with an index of refraction of ${\varepsilon _r} = 13$ and a dimension of $0.3a \times 0.3a.$

Fig. 2. Calculated band diagram for the TE mode in a dielectric-rod waveguide structure for ${M_0} = 5$ and $M = 5.$

Download Full Size | PDF

For this calculation, $N = 0$ (because there is no periodicity in y-direction) has been set. Furthermore, a convergence test reveals that the truncations at ${M_0} \ge 2$ and $M \ge 5$ are necessary for accurate results. The dispersion diagram, i.e., normalized frequency versus normalized wavevector, of the TE guided mode is plotted in Fig. 2 in which the accuracy of our proposed method has been demonstrated through results that closely match those obtained from FDTD method presented in [44]. Moreover, Fig. 3 illustrates the relative magnitudes of the electromagnetic field components for TE guided mode of the investigated dielectric-rod waveguide at free space wavelength of $600\,\textrm{nm}$ by use of our method. As expected, the dominant field components are ${E_y}, {H_z}$ and ${H_x}.$

Fig. 3. Calculated Relative magnitude distribution of electric and magnetic field components of the TE mode for the waveguide shown in inset of Fig. 1 for ${M_0} = 5$ and $M = 5,$ when the wavelength of light is $600\,\textrm{nm}.$

Download Full Size | PDF

The above example possesses an inherently biperiodic layer suspended in air. However, using the idea of periodic repetition [10,21], we have generalized the developed method to the analysis of some kinds of non-periodic structures. As second example, the proposed method is applied to calculate the natural propagation modes of a typical optical waveguide in x-direction on a standard 220nm silicon-on insulator platform with a silica $(Si{O_2})$ substrate and cladding layers of thickness of $2\,\mathrm{\mu}\textrm{m}$ and silicon $(Si)$ core of finite width $W = 500\,\textrm{nm,}$ shown in Fig. 4(a). It should be noted that the refractive indices of $Si$ and $Si{O_2}$ respectively, are taken as $3.47$ and $1.44$ at the wavelength of ${\lambda _0} = 1.55\,\mathrm{\mu m,}$ where their material dispersions are also considered throughout this paper.

Fig. 4. (a) The 3D view of a strip waveguide embedded in $Si{O_2}$ background. (b) Original waveguide after being repeated with period ${L_y}$ along the y-axis.

Download Full Size | PDF

Figure 4(b) shows the original waveguide being repeated along the y-axis with a period of ${L_y}.$ For adequately large repetition periods, the solutions of the periodically repeated waveguide converge to those of the original problem. Therefore, we concentrate our attention on the analysis of the biperiodic structure, illustrated in Fig. 4(b). To this end, we set ${L_y} = 5\,W$ to minimize the coupling between the waveguides and $M = 0$ (because there is no periodicity in x-direction). In most practical cases, including this example, we found that ${M_0} = 10$ ensures very accurate numerical results. The truncation order $N$ is also determined through a convergence test, where we demonstrate the convergence rate of the effective index ${n_{eff}}$ of the structure for the lowest TE-like and TM-like mode solutions versus N at the free space wavelength of $1.55\,\mathrm{\mu}\textrm{m},$ and calculated to be $N > 30$ and $N > 20$, respectively.

By a 3D analysis ($M = 0,$ $N = 40$ and ${M_0} = 10$ have been set) of the structure of the second example, the effective indices of refraction for the lowest TE-like and TM-like modes have been calculated over the broad wavelength range of $0.2\,\mathrm{\mu}\textrm{m} - 2\,\mathrm{\mu}\textrm{m}$ as illustrated in Fig. 5. In justification of our obtained results, the standard finite-difference frequency-domain method (FDFD) using Lumerical Mode Solutions is employed. An excellent agreement is observed between the two results. The field distribution of this structure has been also calculated and is visualized in Figs. 6 and 7 for fundamental TE-like and TM-like modes, respectively. The simulations are conducted via our method at the free space wavelength of $1.55\,\mathrm{\mu}\textrm{m}\textrm{.}$ As can be seen, the dominant components for each mode are maximum in the core and decays exponentially towards the vicinity.

Fig. 5. Dispersion curve for the lowest (a) TE-like and (b) TM-like mode of the structure shown in Fig. 4(b)

Download Full Size | PDF

Fig. 6. Field distribution for the lowest TE-like guided mode of the structure shown in Fig. (4b) for $W = 500\,\textrm{nm}$ and ${\lambda _0} = 1.55\,\mathrm{\mu}\textrm{m}\textrm{.}$

Download Full Size | PDF

Fig. 7. Field distribution for the lowest TM-like guided mode of the structure shown in Fig. (4b) for $W = 500\,\textrm{nm}$ and ${\lambda _0} = 1.55\,\mathrm{\mu}\textrm{m}.$

Download Full Size | PDF

The above typical examples were provided only to justify the validity of the proposed approach. As a more complicated and practical example of application, the presented method was used to determine the guided modes of a subwavelength grating (SWG)-based multimode waveguide for which the convergence of the obtained results and the overall computational complexity is numerically studied. Figure 8(a) shows schematically the structure under investigation widely studied in literature [45,46]; a periodic array of rectangular silicon pillars by assuming longitudinal period ${\mathrm{\Lambda }_\textrm{x}}\textrm{ = 180}\,\textrm{nm,}$ duty cycle $\textrm{D = 0}\textrm{.5,}$ width ${\textrm{W}_{\textrm{MMI}}}\textrm{ = 3}\textrm{.5}\,\mathrm{\mu}\textrm{m}$ and the standard height of $\textrm{h = 220}\,\textrm{nm}$ embedded in $\textrm{Si}{\textrm{O}_\textrm{2}}$ background. According to the periodic repetition concept, we can repeat the structure shown in Fig. 8(a) in y-direction with ${\mathrm{\Lambda }_\textrm{y}}\textrm{ = 3}{\textrm{W}_{\textrm{MMI}}}$ to obtain a biperiodic structure in y- and x-directions (what has been done in Fig. 8(b)). Using the technique described in the present paper, four first guided modes of this waveguide are solved at $\lambda = 1.55\,\mathrm{\mu}\textrm{m}$ for which the results are presented in Fig. 9, where the convergence rate of the effective index $n_{eff}^i$ corresponding to i th mode is plotted versus both the number of retained space harmonics along x $\textrm{(P = 2M + 1)}$ and y $\textrm{(Q = 2N + 1)}$ directions. It can be seen that after $\textrm{P = 5}\,\textrm{(5,}\,21\textrm{,}\,17\textrm{)}$ and $\textrm{Q = 7}\,\textrm{(11,}\,15\textrm{,}\,23\textrm{)}$ space harmonic terms are used for the first (second, third, fourth) mode the effective index converges to its final value. It should be noticed that in this calculation ${M_0},$ i.e., the number of Legendre basis functions, is fixed at 10 which is high enough to ensure the stability of the numerical results. Additionally, the computation times corresponding to the proposed polynomial-based approach and the PWEM algorithm (with parameters P, Q and R as the number of retained space harmonics along x, y and z directions, implemented and described in details in recent paper [21]) to achieve an accuracy of 1% in the final values of effective indices of four first guided modes of the structure shown in Fig. 8(a) at $\lambda = 1.55\,\mu m$ are given in Table 1. All the numerical results are obtained by using MATLAB on a PC (Intel Core i7-6700K CPU at 4.00GHz and 48 GB of RAM). In comparison with Fourier-based analytical techniques, e.g., plane wave expansion by using the Fourier expansion in an appropriately defined supercell, our proposed method is applied only along the transverse direction (y) and the longitudinal variation of the fields is analytically incorporated into V(z) function introduced in Eq. (7). The size of the involved matrices is therefore smaller than the size of the supercell matrices to which three-dimensional Fourier expansion is applied. As solving eigenvalue problems for matrices with order n typically requires O (${n^3}$) operations, reduction in the size of the involved matrices usually plays the decisive role in determining the overall run time. This difference is especially pronounced for calculation of dispersion diagram when a large number of frequency points is to be considered. Table 1 clearly shows that the proposed method can outperform the PWEM. However, to further demonstrate the superiority of the proposed method, the achieved relative error in computing the effective indices of four first eigenmodes is plotted in Fig. 10 versus the elapsed run time of our method (blue solid line) and the PWEM (red solid line). In this comparison, the total number of space harmonics corresponding to the methods was kept large enough for each mode to ensure the error of less than 1%. The reference (final) values for the calculation of the relative error are conducted in Table 1. The accuracy of the proposed method and the PWEM is improved by increasing the total number of space harmonics. On the other hand, as the total number of space harmonics increases, the computation time is drastically increased, whereas the proposed method yields converged results efficiently and the relative error of using PWEM happens to be worse than that of ours calculated by using the proposed method. This figure clearly demonstrates that our method outpaces the PWEM in achieving the lower error level and, in addition, shows that our method shares good convergence properties, far better than those of PWEM, due to the fact that Legendre polynomial expansions of smooth functions (here is V(z) function introduced in Eq. (7)) converge rapidly. This point was further explained in Section 1, where the properties of Legendre polynomials and the Gibb’s phenomenon problem were briefly discussed. Moreover, Figs. 11(a)–(11(d)) illustrate the normalized spatial distribution of ${E_y}$ for the first, second, third and fourth modes of the SWG waveguide, respectively. To generate these spatial distributions, the method described in section 2 has been used for P = 11, Q = 17, ${M_0} = 10$ and $\lambda = 1.55\,\mathrm{\mu}\textrm{m}.$

Fig. 8. (a) The 3D view of a SWG multimode waveguide composed of arrays of rectangular silicon pillars embedded in $\textrm{Si}{\textrm{O}_\textrm{2}}$ background. For this structure the following parameters are assumed: ${\Lambda _x} = 180\,nm, {W_{MMI}} = 3.5\,\mu m,$ $h = 220\,nm$ and $D = 0.5$ (b) Original waveguide after being repeated with period ${\Lambda _y}$ along the y-axis.

Download Full Size | PDF

Fig. 9. Effective indices of the four first guided modes of the SWG multimode waveguide shown in Fig. 8(a) at ${\lambda _0} = 1.55\,\mathrm{\mu}\textrm{m}$ versus the retained space harmonics P and Q in applying our method with ${M_0} = 10.$

Download Full Size | PDF

Fig. 10. Computational efficiency of the proposed method and the PWEM for four first guided modes of the SWG multimode waveguide shown in Fig. 8(a) at ${\lambda _0} = 1.55\,\mathrm{\mu}\textrm{m}\textrm{.}$

Download Full Size | PDF

Fig. 11. Spatial distributions of ${E_y}$ for first, second, third and fourth modes of SWG waveguide shown in Fig. 8(a) at ${\lambda _0} = 1.55\,\mathrm{\mu}\textrm{m}$ visualized in (a), (b), (c) and (d), respectively.

Download Full Size | PDF

Table 1. Required Run Times of the Proposed Method and the PWEM to Achieve an Accuracy of 1% in the Reference (Final) Values of Effective Indices of Four First Guided Modes of the SWG Multimode Waveguide at $\mathrm{\lambda =\ 1}\textrm{.55}\,\mathrm{\mu}\textrm{m}$

View Table

4. Conclusion

In this paper, we present a novel method for extracting guided and leaky modes within layered biperiodic waveguides, employing the Legendre polynomial expansion of electromagnetic fields. By solving the generic Helmholtz equation (Eq. (7a)) that governs the behavior of mixed polarizations, subject to appropriate boundary conditions (Eqs. (11) and (15)), within the complete space defined by Legendre polynomial basis functions, our approach yields algebraic dispersion equations (Eqs. (10d), (14) and (17)) rather than transcendental ones. This characteristic facilitates easier manipulation and unconditional stability, ensuring accurate results with minimal computational resources. Conversely, traditional modal expansion methods face numerical instability when retaining a large number of spatial harmonics in the analysis. The validity of our proposed approach is confirmed through comparison with results derived from finite difference methods and PWEM, showcasing strong agreement in all provided examples. Notably, our method demonstrates efficacy in analyzing the hybrid guided modes of general 3D nonperiodic waveguides through the concept of periodic repetition, as evidenced in the third example. Furthermore, the extensibility of this technique to waveguides comprised of anisotropic layered media holds promise for future applications.

Our findings underscore the robustness and applicability of the Legendre polynomial expansion method, offering a stable, accurate and efficient alternative for analyzing layered biperiodic waveguides. This method's potential for broader utilization across various waveguide configurations, including anisotropic layered media, suggests promising avenues for further exploration and advancement in this field.

Disclosures

The authors declare no conflicts of interest.

Data availability

Simulation details and associated data are available free of charge from authors upon reasonable request.

References

1. H. Mosallaei and Y. Rahmat-Samii, “Periodic bandgap and effective dielectric materials in electromagnetics: characterization and applications in nanocavities and waveguides,” IEEE Trans. Antennas Propagat. 51(3), 549–563 (2003). [CrossRef]

2. H. Y. Yang, “Finite difference analysis of 2-D photonic crystals,” IEEE Trans. Microwave Theory Techn. 44(12), 2688–2695 (1996). [CrossRef]

3. G. Pelosi, A. Cocchi, and A. Monorchio, “A hybrid FEM-based procedure for the scattering from photonic crystals illuminated by a Gaussian beam,” IEEE Trans. Antennas Propagat. 48(6), 973–980 (2000). [CrossRef]

4. B. P. Hiett, J. M. Generowicz, S. J. Cox, et al., “Application of finite element methods to photonic crystal modelling. IEE Proceedings-Science,” IEE Proc.: Sci., Meas. Technol. 149(5), 293–296 (2002). [CrossRef]

5. M. D. Feit and J. A. Fleck, “Light propagation in graded-index optical fibers,” Appl. Opt. 17(24), 3990 (1978). [CrossRef]

6. W. P. Huang and C. L. Xu, “Simulation of three-dimensional optical waveguides by a full-vector beam propagation method,” IEEE J. Quantum Electron. 29(10), 2639–2649 (1993). [CrossRef]

7. C. Vassallo, “1993–1995 Optical mode solvers,” Opt. Quant. Electron 29(2), 95–114 (1997). [CrossRef]

8. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quant. Electron 26(3), S113–S134 (1994). [CrossRef]

9. R. Scarmozzino, A. Gopinath, R. Pregla, et al., “Numerical techniques for modeling guided-wave photonic devices,” IEEE J. Select. Topics Quantum Electron. 6(1), 150–162 (2000). [CrossRef]

10. M. Akbari, K. Schünemann, and H. Burkhard, “A full-wave three-dimensional analysis of open dielectric waveguides and periodic structures using an equivalent network model,” Optical and quantum electronics. 32, 991–1004 (2000). [CrossRef]

11. M. A. Boroujeni and M. Shahabadi, “Modal analysis of multilayer planar lossy anisotropic optical waveguides,” J. Opt. A: Pure Appl. Opt. 8(10), 856–863 (2006). [CrossRef]

12. P. C. Kendall, M. J. Adams, S. Ritchie, et al., “Theory for calculating approximate values for the propagation constants of an optical rib waveguide by weighting the refractive indices. IEE Proceedings A-Physical Science, Measurement and Instrumentation,” Management and Education-Reviews. 8(134), 699–702 (1987).

13. R. Mittra, Y. L. Hou, and V. Jamnejad, “Analysis of open dielectric waveguides using mode-matching technique and variational methods,” IEEE Trans. Microwave Theory Techn. 28(1), 36–43 (1980). [CrossRef]

14. C. Vassallo, “Optical waveguide concepts,” Optical Wave Sciences and Technology. 1, 5 (1991).

15. JB. Pendry, “Photonic band structures,” J. Mod. Opt. 41(2), 209–229 (1994). [CrossRef]

16. ZY Li and LL. Lin, “Photonic band structures solved by a plane-wave-based transfer-matrix method,” Phys. Rev. E 67(4), 046607 (2003). [CrossRef]

17. E Noponen and J. Turunen, “Eigenmode method for electromagnetic synthesis of diffractive elements with three-dimensional profiles,” J. Opt. Soc. Am. A 11(9), 2494 (1994). [CrossRef]

18. P. Lalanne, “Improved formulation of the coupled-wave method for two-dimensional gratings,” J. Opt. Soc. Am. A 14(7), 1592 (1997). [CrossRef]

19. L. R. Lewis and A. Hessel, “Propagation characteristics of periodic arrays of dielectric slabs,” IEEE Trans. Microwave Theory Techn. 19(3), 276–286 (1971). [CrossRef]

20. ST Peng, HL Bertoni, and T. Tamir, “Analysis of periodic thin-film structures with rectangular profiles,” Opt. Commun. 10(1), 91–94 (1974). [CrossRef]

21. K Arik, M Akbari, and A. Khavasi, “Design and simulation of a dynamically tunable 1× 2 power splitter using MZI configuration in the telecom regime,” Opt. Commun. 2(4), 956 (2023). [CrossRef]

22. RC. Rumpf, “Engineering the dispersion and anisotropy of periodic electromagnetic structures,” InSolid State Physics 66, 213–300 (2015). [CrossRef]

23. HS Sözüer, JW Haus, and R. Inguva, “Photonic bands: Convergence problems with the plane-wave method,” Phys. Rev. B 45(24), 13962–13972 (1992). [CrossRef]

24. M. Webb, V. Coppé, and D. Huybrechs, “Pointwise and uniform convergence of Fourier extensions,” Constructive approximation. 52(1), 139–175 (2020). [CrossRef]

25. P Lalanne and GM. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A 13(4), 779 (1996). [CrossRef]

26. RC Rumpf and JJ. Pazos, “Optimization of planar self-collimating photonic crystals,” J. Opt. Soc. Am. A 30(7), 1297–1304 (2013). [CrossRef]

27. S. Guo and S. Albin, “Simple plane wave implementation for photonic crystal calculations,” Opt. Express 11(2), 167–175 (2003). [CrossRef]

28. S. Johnson and J. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis,” Opt. Express 8(3), 173–190 (2001). [CrossRef]

29. M. Chamanzar, K. Mehrany, and B. Rashidian, “Legendre polynomial expansion for analysis of linear one-dimensional inhomogeneous optical structures and photonic crystals,” J. Opt. Soc. Am. B 23(5), 969 (2006). [CrossRef]

30. M. Shahabadi, S. Atakaramians, and N. Hojjat, “Transmission line formulation for the full-wave analysis of two-dimensional dielectric photonic crystals,” IEE Proc.: Sci., Meas. Technol. 151(5), 327–334 (2004). [CrossRef]

31. L. M. Delves and J. N. Lyness, “A numerical method for locating the zeros of an analytic function,” Mathematics of computation. 21(100), 543–560 (1967). [CrossRef]

32. L. C. Botten, M. S. Craig, and R. C. McPhedran, “Complex zeros of analytic functions,” Comput. Phys. Commun. 29(3), 245–259 (1983). [CrossRef]

33. D Gottlieb and SA. Orszag, “Numerical analysis of spectral methods: theory and applications,” Society for Industrial and Applied Mathematics; (1977).

34. H. Wang and S. Xiang, “On the convergence rates of Legendre approximation,” Math. Comp. 81(278), 861–877 (2012). [CrossRef]

35. H. Wang, “How much faster does the best polynomial approximation converge than Legendre projection?” Numer. Math. 147(2), 481–503 (2021). [CrossRef]

36. H. Wang, “New error bounds for Legendre approximations of differentiable functions,” J. Fourier Anal. Appl. 29(4), 42 (2023). [CrossRef]

37. LN. Trefethen, Approximation Theory and Approximation Practice, Extended Edition. (Society for Industrial and Applied Mathematics; 2019).

38. T. N. Phillips, “On the Legendre coefficients of a general-order derivative of an infinitely differentiable function,” IMA J. Numer. Anal. 8(4), 455–459 (1988). [CrossRef]

39. K. Mehrany and B. Rashidian, “Polynomial expansion for extraction of electromagnetic eigenmodes in layered structures,” J. Opt. Soc. Am. B 20(12), 2434 (2003). [CrossRef]

40. A. Khavasi, K. Mehrany, and B. Rashidian, “Three-dimensional diffraction analysis of gratings based on Legendre expansion of electromagnetic fields,” J. Opt. Soc. Am. B 24(10), 2676 (2007). [CrossRef]

41. V. Lancellotti and R. Orta, “Guided waves in layered cubic media: Convergence study of a polynomial expansion approach,” J. Acoust. Soc. Am. 104(5), 2638–2644 (1998). [CrossRef]

42. M. Chamanzar, K. Mehrany, and B. Rashidian, “Planar diffraction analysis of homogeneous and longitudinally inhomogeneous gratings based on Legendre expansion of electromagnetic fields,” IEEE Trans. Antennas Propagat. 54(12), 3686–3694 (2006). [CrossRef]

43. A. Khavasi, A. K. Jahromi, and K. Mehrany, “Longitudinal Legendre polynomial expansion of electromagnetic fields for analysis of arbitrary-shaped gratings,” J. Opt. Soc. Am. A 25(7), 1564 (2008). [CrossRef]

44. S. Fan, J. N. Winn, A. Devenyi, et al., “Guided and defect modes in periodic dielectric waveguides,” J. Opt. Soc. Am. B 12(7), 1267 (1995). [CrossRef]

45. P Cheben, R Halir, DR. Schmid, et al., “Subwavelength integrated photonics,” Nature 560(7720), 565–572 (2018). [CrossRef]

46. Robert Halir, Alejandro Ortega-Moñux, Daniel Benedikovic, et al., “Subwavelength-grating metamaterial structures for silicon photonic devices,” Proc. IEEE 106(12), 2144–2157 (2018). [CrossRef]

	Required Run Times (s)		Reference Values
Modes	Our Method	PWEM	Our Method	PWEM
1st Mode	8	37	2.1118	2.1111
2nd Mode	14	42	2.0518	2.0515
3rd Mode	32	72	1.9501	1.9509
4th Mode	41	93	1.8080	1.8087

Polynomial expansion method for full wave three-dimensional analysis of dielectric waveguides and periodic structures

Abstract

1. Introduction

2. Formulation

2.1. Maxwell's equations within a biperiodic layer

2.2. Polynomial expansion of electromagnetic fields inside the biperiodic layer

2.3. Applying boundary conditions

3. Numerical results

4. Conclusion

Disclosures

Data availability

References

Data availability

Cited By

Figures (11)

Tables (1)

Equations (38)

Optics Express