Efficient smoothed finite element time domain analysis for photonic devices

Khaled S. R. Atia; A. M. Heikal; S. S. A. Obayya

doi:10.1364/OE.23.022199

1. Introduction

Finite difference method (FDM) is one of the most widely used numerical tools for optical circuits analysis [1–6]. Although simple and efficient, FDM suffers from the insufficient discretization of irregularly-shaped geometries (staircasing problem). Alternatively, the finite element method (FEM) is able to divide the computational domain into triangular elements and to accurately model complex geometries. Several FEM analyses rely on the second order triangular element (T6) due to the poor accuracy of the first order triangular element (T3) [7–10]. However, for the conventional T6 element, the additional nodes located midway the distance between the pairs of triangle vertexes cannot be adjusted to set on the boundary of a curved structure. This problem was solved by using curvilinear elements [11]. Compared with conventional elements, curvilinear element approach converges faster for curved structures. However, the T3 element is often preferable because of its computational simplicity and easiness for adaptive mesh refinement. Moreover, most of mesh generation software’s rely on first order triangular elements.

As a consequence, several research efforts have been dedicated to increase the accuracy of the T3 element. Recently, a new technique termed ”smoothed Galerkin” has been introduced [12]. This new technique relies on weakening the weak form (W²) as in meshfree methods [13]. Beside approximating the field, W² form approximates also the field derivatives via the gradient smoothing technique [14]. The high accuracy and the less sensitivity to mesh distortion are the most important features of W² form. The accuracy of any method relies on Galerkin form can be improved by using the smoothed Galekrin form. Based on the smoothed Galerkin technique, a novel class of so-called smoothed finite element method (SFEM) has been introduced for various mechanics problems [15]. The shape function derivatives at element boundaries, in the conventional FEM, are numerically discontinuous which causes FEM poor performance. This disadvantage in FEM has been smartly overcome in SFEM by approximating the derivatives at the element boundaries to be continuous.

Time domain analysis is indispensable for design and actual performance prediction of photonic devices. Relying on the FDM, a powerful and versatile technique called finite difference time domain (FDTD) has been introduced [1]. However, besides the already mentioned staircasing problem, FDTD experiences stability problems where the time step size should satisfy the so-called Courant-Friedrich-Levy condition [1] which makes the FDTD greatly sensitive to the step size choice. Such condition reduces the time step size used and consequently increases the overall computational time. Various variations of the FDTD have been proposed to overcome the stability problem. Based on the time domain beam propagation method (TD-BPM), a FDTD-BPM has been introduced. TD-BPM removes the fast varying field component leaving only the slowly varying field to be simulated [16,17]. TD-BPM depends on simulating the wave equation. Also a complex envelop alternating direction implicit technique (CE-ADI-FDTD) [5] and split step technique (SS-FDTD) [18] have been proposed to employ large time step size. Unlike the FDTD-BPM, CE-ADI-FDTD and SS-FDTD rely on simulating the coupled Maxwell’s equations which require to solve 3 coupled equations at each time step (more computational cost). Alternatively, The FEM provides a genuine solution to the staircasing problem. A variety of finite element time domain (FETD) techniques has been introduced over the past decade. Relying on the edge basis functions, a vectorial solution to the time domain Maxwell’s equations has been proposed [19]. Unfortunately, it yields a FDTD like conditionally stable algorithm. On the other hand, a solution to the curl-curl wave equation which solves one of Maxwell’s equations field vectors has been presented [20]. This method is unconditionally stable, however, it produces a large system of equations that needs to be solved at each time step. Therefore, it is considered computationally expensive. Dissimilar to the vectorial shape function, the nodal shape function can be used to solve only one field component beside it produces a well-scaled matrices [7, 21]. Hence, the nodal shape functions is considered computationally efficient. Based on the nodal shape function, a solution to the time domain Maxwell’s equations is made through the application of TD-BPM which permits the use of large time step size. FETD methods can be categorized into explicit and implicit schemes. In the explicit scheme, a mass lumping technique is used and the slowly varying wave equation is temporally discretized using central finite difference formula [21]. In this scheme, no linear system of equations is solved at each time step. However, one of the major disadvantages of this scheme is that it yields a conditionally stable time-marching algorithm where the time step size should satisfy the Courant-Friedrich-Levy condition to ensure the stability of the algorithm. Moreover, the mass lumping technique used severely affects the accuracy of the method. The second scheme discretizes the slowly varying wave equation in time domain yielding a linear system of equations to be solved at each time step. The advantage of the second scheme is that the produced time-marching algorithm is unconditionally stable which allows the use of moderate time step size. This scheme can be implemented by either discretizing the slowly varying wave equation using the unconditional stable Newmark-Beta scheme [22] or by padé approximation [7]. The latter algorithm is less memory demanding as, unlike Newmark-Beta method, it computes the field in the next time step using only the previous one.

In this paper, SFEM is newly-proposed to analyse photonic circuits in time domain through utilising the TD-BPM. The presented method uses Padé approximation and is abbreviated as SFETD-BPM-T3. The introduced method is, to the best of authors knowledge, the first SFEM for time domain analysis. The rest of this paper is organized as follows. The time domain analysis and the SFEM formulation is presented in Section 2. In Section 3, the numerical dispersion of the proposed method is assessed through a single mode slab waveguide. Accuracy and convergence of SFETD-BPM-T3 are demonstrated through computing the reflection of an optical grating structure. Moreover, the efficiency of handling complex geometries while maintaining high accuracy is investigated by a photonic crystal (PC) cavity. Furthermore, a new formulation for plasma material based on padé approximation is introduced and investigated by calculating the reflection and transmission of plasma material. Additionally, the material dispersion is used to enhance the quality factor of PC cavity. Finally, conclusions are drawn in Section 4.

2. Analysis

2.1. Time domain equations

Consider the two-dimensional (2D) computational domain (Ω) shown in Fig. 1, where Ω is on the yz-plane and the variation on x direction is neglected. By considering the perfectly matched layer (PML) boundary conditions reported in [7], the 2D wave equation under slowly varying approximation can be written in the following form

\nabla [D] \nabla^{T} ϕ + s \frac{ω_{o}^{2} q}{c^{2}} ϕ - 2 j s \frac{ω_{o} q}{c^{2}} \frac{\partial ϕ}{\partial t} - s \frac{q}{c^{2}} \frac{\partial^{2} ϕ}{\partial t^{2}} = 0,

with

[D] = [\begin{matrix} p \frac{s_{y}^{2}}{s} & 0 \\ 0 & p \frac{s_{z}^{2}}{s} \end{matrix}],

\nabla = [\begin{matrix} \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \end{matrix}],

ϕ = E_{x}, p = 1, q = n^{2}, for TE modes

ϕ = H_{x}, p = 1 / n^{2}, q = 1, for TM modes,

s = {\begin{array}{l} 1 - j \frac{3 c}{2 ω_{o} n d} {(\frac{ρ}{d})}^{2} \ln (\frac{1}{R}), & in PML region \\ 1, & in non-PML region, \end{array}

where E_x and H_x are the slowly varying electric and magnetic fields in x direction, respectively, ω_o is the carrier center frequency, c is the speed of light in free space, n is the material refractive index distribution, t is the time, d is the PML thickness, R is the theoretical reflection coefficient at the boundary of the computational domain, and ρ is the distance from the beginning of the PML. s, s_y, and s_z are the PML parameters [7] where s_y and s_z are defined in Table 1.

Fig. 1 The computational domain (Ω) divided into triangular elements (solid black lines). On top of the element mesh, smoothing domains is constructed (dashed red lines) to calculate the smoothed stiffness matrix. The filled circles (•) represent the field nodes and the empty circles (○) represent the centroid of the element.

Download Full Size | PDF

Table 1. The PML parameters s_y and s_z definition in each region.

View Table | View all tables in this article

2.2. Galerkin weak form

Applying Galerkin weak form to Eq. (1) gives

\begin{array}{l} - \iint_{Ω} {B} [D] {B}^{T} {ϕ} d Ω + \frac{ω_{o}^{2}}{c^{2}} \iint_{Ω} sq {N} {N}^{T} {ϕ} d Ω \\ - 2 j \frac{ω_{o}}{c^{2}} \iint_{Ω} sq {N} {N}^{T} \frac{\partial {ϕ}}{\partial t} d Ω - \frac{1}{c^{2}} \iint_{Ω} sq {N} {N}^{T} \frac{\partial^{2} {ϕ}}{\partial t^{2}} d Ω = 0, \end{array}

with

{B} = {\nabla N} .

By discretizing Ω into N_e first order triangular elements and letting

[M] = \iint_{Ω} sq {N} {N}^{T} d Ω = \sum_{N_{e}} \iint_{Ω_{e}} sq {N} {N}^{T} d Ω_{e},

and

[K] = - \iint_{Ω} sq {B} [D] {B}^{T} d Ω = - \sum_{N_{e}} \iint_{Ω_{e}} sq {B} [D] {B}^{T} d Ω_{e},

we have

- \frac{1}{c^{2}} [M] \frac{d^{2} {ϕ}}{d t^{2}} - 2 j \frac{ω_{o}}{c^{2}} [M] \frac{d {ϕ}}{d t} + ([K] + \frac{ω_{o}^{2}}{c^{2}} [M]) {ϕ} = {0},

where {N} is a vector contains the FEM-T3 nodal shape functions for all the nodes in Ω, Ω_e is the element domain, {ϕ} is the global electric or magnetic field vector, and {0} is a null vector. [M] and [K] are the standard FEM mass and stiffness matrices, respectively.

2.3. Smoothed Galerkin procedure

On top of the FEM element mesh, the computational domain Ω is divided into N_s smoothing domains $Ω_{k}^{s}$ . These smoothing domains are created by simply connecting the centroid of the element to the two end points of each edge as shown in Fig. 1. The reason for choosing the smoothing domains over the edges is that we need to overcome the discontinuity of the shape function derivatives over the edges. We should address that the FEM poor accuracy is caused by the discontinuity of the shape function derivatives. Such discontinuity increases the numerical error as the wave equation assumes the continuity of the function and its derivatives. In our method, a smoothed stiffness matrix [K̄] is computed over the smoothing instead of the conventional stiffness matrix [K]. Since the smoothing domain is shared between two elements, an ”inter-element” interaction between the two adjacent elements should be considered. Hence, [K̄] can be expressed as

[\bar{K}] = - \iint_{Ω} sq {B} [D] {B}^{T} d Ω = - \sum_{N_{s}} \iint_{Ω_{k}^{s}} {B} [D] {B}^{T} d Ω_{k}^{s},

According to the gradient smoothing technique [13], {B} is approximated using the following equation

{B} = {\nabla N} \approx \iint_{Ω_{k}^{s}} \hat{W} {\nabla N} d Ω_{k}^{s} .

Integrating Eq. (13) by parts gives

{B} \approx - \iint_{Ω_{k}^{s}} \nabla \hat{W} {N} d Ω_{k}^{s} + \int_{Γ_{k}^{s}} \hat{W} {N} \vec{u} d Γ_{k}^{s},

where

Γ_{k}^{s}

is the boundary of

Ω_{k}^{s}

, u⃗ is the outward unit normal vector, and Ŵ is the smoothing function. Ŵ should satisfy the following conditions; (1) it should be positive over

Ω_{k}^{s}

and zero otherwise, (2) it should be centered at the edge, and (3) it should have a unity integral over

Ω_{k}^{s}

. Using such shape function derivative approximation leads to a special form of Galerkin weak form termed ”smoothed Galerkin form”. In this paper we choose (

\hat{W} = 1 / A_{k}^{s}

) for simplicity where

A_{k}^{s}

is the area of the smoothing domain

Ω_{k}^{s}

. The smoothing function can be seen as an averaging function. Since the derivative of the FEM-T3 shape function is constant (in the element domain), using a higher order polynomial function as a smoothing function wouldn’t affect the accuracy. However, if the shape function derivative is of a higher order polynomial, a smoothing function of the same order should be used to get better performance [12,13,15,23]. It is worth mentioning that the smoothing function causes no additional burden as it doesn’t affect the number of degree of freedoms. It only describes how the shape function derivative is smoothed at the edges. By choosing a constant Ŵ in our case, the surface integral in Eq. (14) vanishes leaving only the line integral over

Γ_{k}^{s}

. This line integral can be written in the following form

{B} = \frac{1}{A_{k}^{s}} \sum_{i}^{n_{s}} \int_{Γ_{k, i}^{s}} {\hat{N}} d Γ_{k, i}^{s} [\begin{matrix} u_{y, i} & u_{z, i} \end{matrix}],

where n_s is the total number of the smoothing domain boundary segments, u_y,i and u_z,i are the components of the outward unit normal vector on the boundary segment

Γ_{k, i}^{s}

in y and z directions, respectively. n_s equals 4 for inner smoothing domains and 3 for the boundary smoothing domain. The line integral of the shape function can be evaluated by Gauss quadrature technique [15]. By using Eq. (15), we assume the derivatives of the shape functions are continuous over

Ω_{k}^{s}

. This assumption allows the method to produce good accuracy for coarse meshes as it is shown in the results section. Also, the local smoothed stiffness matrix is now 4 × 4 because the number of nodes related to the smoothing domain is 4 (i.e. nodes 1,2,3, and 4 shown in Fig. 1) while the conventional FEM-T6 and FEM-T3 local stiffness matrix are 6 × 6 and 3 × 3, respectively. Therefore, it is expected that the SFEM-T3 matrix bandwidth would be slightly larger than FEM-T3 matrix bandwidth and yet smaller that FEM-T6. Since prediction of the matrix bandwidth is mainly dependent on the mesh used and nodes ordering, we can elucidate the increase in matrix bandwidth by considering a node shared between 6 elements as in the structured mesh shown in Fig. 2. The contribution of this node extends to 6 nearby nodes for FEM-T3 case and to 18 nearby nodes for FEM-T6 (the blue points shown in Fig. 2(a) and 2(b), respectively) while for SFEM-T3, the contribution extends to 12 nearby nodes (the red points added to the blue points shown in Fig. 2(a)). Hence, the estimated matrix bandwidth for FEM-T3, FEM-T6, and the presented method is 7, 19, and 13,respectively. Such bandwidth increase is the price for the significant improvement of the accuracy levels of the proposed method as it is shown in the results section. Another advantage added by the presented method is the mesh distortion immunity. For the distorted two adjacent elements shown in Fig. 3, the conventional stiffness matrix is proportional to the inverse of the element area as described by Eq. (10). Therefore, the value of the stiffness matrix related to point 4 is larger than the related to point 3. On contrary, the smoothed stiffness matrix is proportional to the inverse of the area of the smoothed domain (as described by Eq. (12)) which is shared between the two adjacent elements. Therefore, the values of the stiffness matrix related to the points 4 and 3 are of the same level. This property makes the smoothed matrix computation faster than the conventional one as it is shown in the results section.

Fig. 2 Schematic diagram elucidates the increase in matrix bandwidth of the proposed method over FEM-T3 and FEM-T6.

Download Full Size | PDF

Fig. 3 Two adjacent elements describe the mesh distortion immunity.

Download Full Size | PDF

As [D] is defined for each element and {B} is constant over each smoothing domain, we can rewrite Eq. (12) in the following form

[\bar{K}] = - \sum_{N_{s}} {B} [(A_{k, 1}^{s} {[D]}_{1}) + (A_{k, 2}^{s} {[D]}_{2})] {B}^{T},

where

A_{k, 2}^{s}

and

A_{k, 1}^{s}

are the areas of

Ω_{k, 1}^{s}

and

Ω_{k, 2}^{s}

depicted in Fig. 1, respectively. Moreover, [D]₁ and [D]₂ are the definitions of [D] over the elements 1 and 2 sharing the same smoothing domain

Ω_{k}^{s}

, respectively.

2.4. Time domain discretization

To treat wideband optical pulses, Padé approximation is applied to Eq. (11) leading to the following equation

- 2 j \frac{ω_{o}}{c^{2}} [\hat{M}] \frac{d {ϕ}}{d t} + ([\bar{K}] + \frac{ω_{o}^{2}}{c^{2}} [M]) {ϕ} = 0,

where [M^] = [M] for narrowband optical pulses, and for wideband optical pulses [M^] is given by [7]

[\hat{M}] = [M] = \frac{c^{2}}{4 ω_{o}^{2}} ([\bar{K}] - \frac{ω_{o}^{2}}{c^{2}} [M]) .

Finally, the slowly varying field {ϕ} is discretized in time t using Crank-Nicholson algorithm conducting the following recursive relation

{[R]}_{i} {ϕ}_{i + 1} = {[Q]}_{i} {ϕ}_{i},

with

{[R]}_{i} = - 2 j \frac{ω_{o}}{c^{2}} {[\hat{M}]}_{i} + 0.5 Δ t ({[\bar{K}]}_{i} + \frac{ω_{o}^{2}}{c^{2}} {[M]}_{i}),

{[Q]}_{i} = - 2 j \frac{ω_{o}}{c^{2}} {[\hat{M}]}_{i} - 0.5 Δ t ({[\bar{K}]}_{i} + \frac{ω_{o}^{2}}{c^{2}} {[M]}_{i}),

where Δt is the time step size, {ϕ}_i and {ϕ}_i+1 are the field vectors at ith time step and (i +1)th time step, respectively. The previous system of equations is iteratively solved using Bi-CGSTAB algorithm [24]. Unlike the well-known FDTD method [1], this time discretization technique has been previously reported as unconditionally stable technique [7]. The proposed method can be easily extended to three-dimensional problems by smoothing the derivatives over the faces of tetrahedral elements.

3. Simulation results

3.1. Numerical dispersion

The numerical dispersion of the proposed method is assessed through calculating the pulse velocity inside a single-mode waveguide. Therefore, the single-mode slab waveguide shown in Fig. 4(a) is investigated. During the simulation, the time step size Δt and the central wavelength λ_o are fixed to 1 fs and 1.55μm, respectively. The waveguide is excited by a pulse with a Gaussian profile in z direction and a transverse profile ϕ_o(y) corresponding to the fundamental mode of the planner waveguide. Therefore, the initial field at t = 0 is taken as

ϕ (y, z, t = 0) = ϕ_{o} (y) \exp [- {(\frac{z - z_{o}}{W_{o}})}^{2}],

where β is the propagation constant of the fundamental mode, z_o is the initial pulse position, and W_o is the spot size. During the simulation, W_o and z_o are taken as 2.0μm and 20.0μm, respectively. The pulse velocity is calculated by measuring the distance travelled by the optical pulse after 100 fs. Figure 4(b) shows the pulse velocity calculated at different longitudinal step size Δz using the proposed method, FETD-BPM-T3, and FETD-BPM-T6. As may be shown in Fig. 4(b), the calculated pulse velocity using the proposed method converges better than FETD-BPM-T3 and FETD-BPM-T6 to the exact solution. Moreover, it can be noticed that the proposed method is capable of producing accurate results even for large Δz; thanks to the gradient smoothing technique which relaxes the discontinuity in shape function derivatives. Figure 4(b) also gives an indication that the proposed method is more stable than FETD-BPM-T3 and FETD-BPM-T6 with respect to Δz.

Fig. 4 (a) Schematic diagram of a single mode slab waveguide. (b) Numerical dispersion calculated using the proposed method, FETD-BPM-T6 and FETD-BPM-T3 [7].

Download Full Size | PDF

3.2. Optical grating

In order to show the accuracy and the high convergence rate of the proposed method in a more complicated problem, the optical grating structure shown in Fig. 5(a) is considered. The structure is excited as described in the previous example. In this study, 220 time steps are used. The converged results are obtained by the proposed method when the computational domain is divided into 46240 structured non-uniform T3 elements. This discretization produces 23835 unknowns. For a fair comparison, the computational domain is discretized into 5780 non-uniform T6 elements to produce the same number of unknowns. During the simulation run, the reflected pulse is monitored inside the waveguide at point P as shown in Fig. 5(a). By normalizing the Fast Fourier Transform of the reflected data to the input pulse spectrum, the reflection characteristics of the input pulse can be obtained. Figure 5(c) and 5(d) show the TE and TM modes reflection characteristics, respectively. An excellent agreement may be seen between the SFETD-BPM-T3 and FETD-BPM-T6 [7]. Moreover, it may be noticed that the results obtained using FETD-BPM-T3, under the same assumptions, still not converged. The computational time of SFETD-BPM is 829 s while the FETD-BPM-T3 is 1481 s which is approximately twice the SFETD-BPM-T3 computational time. The computational time for FETD-BPM-T6 is 2025 s. Table 2 includes a comparison between the presented method, FETD-BPM-T3, and FETD-BPM-T6 in the computational time required at a certain level of accuracy. The simulations are carried out on HP Z620 Workstation (2 GHz). Figure 5(b) shows the wavelength of maximum reflection calculated using different number of nodes for the TE case. It may be shown that the presented method outperforms the FETD-BPM-T3 in the convergence rate and the accuracy. Moreover, the accuracy of the proposed method is comparable to the FETD-BPM-T6. Figure 5(b) also indicates that the proposed method can produce accurate results even for small number of unknowns (coarse mesh).

Fig. 5 (a) Schematic diagram of the optical grating structure. (b) The wavelength of maximum reflectivity computed using SFETD-BPM-T3, FETD-BPM-T6 and FETD-BPM-T3. (c) Optical gratting Reflection characteristics for TE mode. (d) Optical gratting Reflection characteristics for TM mode.

Download Full Size | PDF

Table 2. Comparison between SFETD-BPM-T3, FETD-BPM-T3, and FETD-BPM-T6 in computing the reflection characteristics of TE-mode and TM-mode.

View Table | View all tables in this article

3.3. Photonic crystal cavity

Modelling complex structures such as PCs requires a large number of elements to represent the structure accurately. In this case, using a higher order element would produces more unknowns. To illustrate this issue, a 5 × 5 PC cavity shown in Fig. 6(a) is studied [5, 25]. The computational domain is discretized to 24604 elements to accurately model the PCs. This discretization produces 12383 unknowns in the case of T3 element and in the case of T6 element it produces 49369 unknowns which is much higher than first order case. The cavity considered consists of dielectric rods surrounded by air with refractive index n_rods = 3.4. The cavity is excited by a point source placed at the center of the cavity. The source has a Gaussian profile in time, and it is defined as

ϕ (x = x_{o}, y = y_{o}, t) = e^{- {(\frac{t - t_{o}}{T_{o}})}^{2}},

where x_o and y_o are the coordinates of the center of the cavity. t_o and T_o are the delay and the width of the Gaussian profile in time, respectively. In this simulation, the wavelength is fixed to 1.5μm, t_o and T_o are set to 90 fs and 30 fs, respectively. Figure 6(b) shows the time variation of E field monitored for 1024 time step (time step = 1 fs) at the center of the cavity. The field starts to grow up till the excitation vanishes. After that, the field decays exponentially and the resonance effect is presented. The normalized resonance frequency (F = a/λ) can be evaluated from Fig. 6(c). Figure 6(d) shows the localized cavity mode formed in the cavity after 1024 fs computed using the proposed method. The quality factor (Q-factor) can then be computed from the ratio of the stored energy to the lost energy after one cycle T = 5.16 fs for the resonance wavelength λ_c = 1.55μm. The Q-factor calculated using FETD-BPM-T6 equals ∼ 183.86 and equals ∼ 179 using FETD-BPM-T3. On the other side, the Q-factor calculated using the proposed method equals ∼ 183.9. By choosing (F_o = 0.37844) reported in [25] as a reference, we can compute F relative error in order to illustrate the convergence and the accuracy of the proposed method. Table 3 shows that the presented method outperforms FETD-BPM-T3 and FETD-BPM-T6 in calculating F of PC cavity. Table 3 also includes the computational time consumed by each method. It maybe noticed that our method is computationally efficient than FETD-BPM-T3 and FETD-BPM-T6.

Fig. 6 (a) Schematic diagram of 5 × 5 PC cavity. (b) Time variation of the envelop of the electric field monitored at the center of the cavity. (c) The spectral distribution of the resonant mode energy for the square cavity computed using the proposed method. (d) The electric field profile of the resonant mode inside the square cavity computed using the proposed method.

Download Full Size | PDF

Table 3. Comparison between the SFETD-BPM-T3, FETD-BPM-T6, and FETD-BPM-T3 in calculating the resonance frequency (F = a/λ) of the PC cavity.

View Table | View all tables in this article

3.4. Dispersive material

Here, the presented method is extended to simulate the propagation in dispersive materials. There are many models that have been proposed for the dispersion effects of the medium [26]. In this example, the plasma model is investigated and TE case is considered. The relative permittivity according to plasma model is given by

ε_{r} (ω) = ε_{\infty} + \frac{ω_{p}^{2}}{j ω (j ω + ν_{c})},

where ω_p and ν_c are the plasma and the damping frequencies, respectively. ε_∞ is the relative permittivity at high frequency. Substituting the plasma model into the slowly varying wave equation Eq. (1) leads to

\nabla D \nabla^{T} E_{x} + s \frac{ω_{o}^{2}}{c^{2}} (ε_{\infty} - \frac{ω_{p}^{2}}{ω_{o}^{2}}) E_{x} - 2 j s \frac{ω_{o} ε_{\infty}}{c^{2}} \frac{\partial E_{x}}{\partial t} - s \frac{ε_{\infty}}{c^{2}} \frac{\partial^{2} E_{x}}{\partial t^{2}} = - \frac{ν_{c} ω_{p}^{2}}{c^{2}} ψ,

where ψ = e^−ν_ctU(t) ⊗ E_x, U(t) is the unit step function and ⊗ denotes convolution. The resultant convolution arises from transforming the plasma model into time domain. The convolution is treated by linearly approximating E_x and analytically solve the other term. Consequently, a recursive relation for ψ can be obtained

{ψ}_{i + 1} = e^{- (ν_{c} + j ω_{o}) Δ t} {ψ}_{i} + \frac{1 - e^{- j ω_{o} Δ t}}{2 j ω_{o}} [{E_{x}}_{i + 1} + e^{- ν_{c} Δ t} {E_{x}}_{i}] .

Discretizing the resultant wave equation using SFEM and applying padé approximation to treat wide band optical pulses, we have

- 2 j \frac{ω_{o}}{c^{2}} [{\hat{M}}_{d}] \frac{d {E_{x}}}{d t} + ([K] + \frac{ω_{o}^{2}}{c^{2}} [G]) {E_{x}} = - \frac{ω_{p}^{2} ν_{c}}{c^{2}} [W] {Ψ},

with

[G] = \sum_{N_{e}} \iint_{Ω_{e}} s (ε_{\infty} - \frac{ω_{p}^{2}}{ω_{o}^{2}}) g {N} {N}^{T} d Ω_{e},

[{\hat{M}}_{d}] = [M] - \frac{c^{2}}{4 ω_{o}^{2}} ([K] - \frac{ω_{o}^{2}}{c^{2}} [G]),

[W] = \sum_{N_{e}} \iint_{Ω_{e}} {N} sg {N}^{T} d Ω_{e}, and g = {\begin{array}{l} 1, & in plasma material, \\ 0, & in non-plasma material . \end{array}

E_x is discretized in time domain using Crank-Nicholson formulation. The developed formulation is more memory-efficient than the previously introduced formulation [10] because it requires to store the field values at the previous time step only. It is also stable for the used PML [7, 8]. The developed formulation is used to simulate a plane wave source placed in air at 80mm from a plasma material of 20mm thickness as shown in Fig. 7. The computational domain is truncated using PML in the longitudinal direction and a perfect magnetic conductor boundary condition is used in the transverse direction. The plasma material is characterized by ω_p = 2π × 28.7 GHz, ν_c = 20 GHz and ε_∞ = 1. Figures 8(a) and 8(b) show the reflection and transmission of the plasma slab, respectively, compared with the analytical solution [27]. It may be seen an excellent agreement between the simulation results using the proposed method and the analytical solution.

Fig. 7 Schematic diagram of the computational domain used to calculate the reflection and the transmission of plasma material.

Download Full Size | PDF

Fig. 8 (a) Reflection and (b) transmission of a plane wave incidents normally from air to 20 mm slab of plasma material.

Download Full Size | PDF

We have to mention that this formulation is general and can be applied for TM cases also (as most of plasmonics and metamaterials application) but in this case E_y must be used instead of E_x. Of course the boundary conditions in y direction must be satisfied. The boundary conditions in this method can be imposed by rewriting {B} and [K̄] in the following forms [15]

{B} = \frac{1}{A_{k}^{s}} (\sum_{i = 1}^{2} \int_{Γ_{k, i}^{s}} {\hat{N}} d Γ_{k, i}^{s} [\begin{matrix} ε_{r 1} u_{y, i} & u_{z, i} \end{matrix}] + \sum_{i = 3}^{4} \int_{Γ_{k, i}^{s}} {\hat{N}} d Γ_{k, i}^{s} [\begin{matrix} ε_{r 2} u_{y, i} & u_{z, i} \end{matrix}]),

[\bar{K}] = - \sum_{N_{s}} {B} [\begin{matrix} A_{k}^{s} & 0 \\ 0 & \frac{A_{k, 1}^{s}}{ε_{r 1}} + \frac{A_{k, 2}^{s}}{ε_{r 2}} \end{matrix}] {B}^{T},

where ε_r1 and ε_r2 are the relative permittivity in

Ω_{k, 1}^{s}

and

Ω_{k, 2}^{s}

, respectively.

3.5. Enhanced cavity

In this section, we investigate the effect of material dispersion on the Q-factor of PC cavities. We start first by computing the quality factor of a 7 × 7 PC square cavity consists of dielectric rods of refractive index n = 3.4 surrounded by air. The Q-factor, for this case, is 1422. Then, we compute the Q-factor of a cavity formed by replacing 4 rods in the second ring by rods of diameter di = 0.1a surrounded by silver so that the total diameter becomes d = 0.2a as depicted in Fig. 9(a). The silver is described using plasma model with ε_∞ = 1, ω_p = 2π × 2074.8 THz and ν_c = 105.2 THz. The Q-factor, in this case, is 4039. This enhancement can be possibly explained as follows. The energy lost in the cavity is affected by the high permittivity of the silver particles. Therefore, the decay rate is reduced as shown in Fig. 9(b). Next, when 4 rods in the inner ring are replaced instead of the second ring rods, in the same manner, the quality factor drops to 247. In this case, silver losses severely increase the total energy lost and consequently reduce the quality factor. Finally, we find that the normalised resonance frequency for the enhanced cavity is 0.3787 while for the conventional cavity is 0.379.

Fig. 9 (a) Schematic diagram of 7 × 7 PC cavity modified by replacing 4 rods in the second ring by rods of diameter di = 0.1a surrounded by silver so that the total diameter becomes d = 0.2a. (b) Time variation of the envelop of the electric field monitored at the center of the cavity for the conventional 7 × 7 PC cavity, 7 × 7 PC cavity with 4 coated rods in the second ring and 7 × 7 PC cavity with 4 coated rods in the first ring.

Download Full Size | PDF

4. Conclusion

An efficient SFEM has been successfully developed. Unlike the conventional FEM, the field derivatives across element interfaces are smoothed through accurately forming ”inter-element” matrices. The proposed method treats the numerical discontinuity in FEM shape function derivatives and hence boosts the accuracy. The proposed method is adopted to TD-BPM to analyse photonic devices in time domain. The SFETD-BPM outperformed the FETD-BPM in calculating the pulse velocity inside a single-mode slab waveguide. Moreover, the optical grating results indicate that our method has a convergence rate comparable to the second order FEM. Furthermore, the presented method has shown a remarkable performance (time and accuracy) in computing the resonance frequency of PC cavities. Moreover, in the regard of extending our work to include dispersive materials, a new dispersive material formulation has been developed and evaluated. The newly-developed formulation has been used to enhance the quality factor of PC cavity.

References and links

1. A. taflove, Computational Electrodynamics: The Finite Difference Time Domain Method (Artech, 1995)

2. M. F. O. Hameed, S. S. A. Obayya, and H. A. El-Mikati, “Passive polarization converters based on photonic crystal fiber with L-shaped core region,” J. Lightwave Technol. 30(3), 283–289 (2012). [CrossRef]

3. F. F. K. Hussain, A. M. Heikal, M. F. O. Hameed, J. El-Azab, W. S. Abdelaziz, and S. S. A. Obayya, “Dispersion characteristics of asymmetric channel plasmon polariton waveguides,” IEEE J. Quantum Electron. 50(6), 474–482 (2014). [CrossRef]

4. A. M. Heikal, M. F. O. Hameed, and S. S. A. Obayya, “Coupling characteristic of a novel hybrid long-range plasmonic waveguide including bends,” IEEE J. Quantum Electron. 49(8), 621–627 (2013). [CrossRef]

5. D. Pinto and S. S. A. Obayya, “Improved complex-envelope alternating-direction-implicit finite-difference-time-domain method for photonic-bandgap cavities,” J. Lightwave Technol. 25(1), 440–447 (2007). [CrossRef]

6. M. Muhammad, S. S. A. Obayya, A. M. Heikal, and M. F. O. Hameed, “Porous core photonic crystal fibre with metal-coated central hole for terahertz applications,” IET Optoelectronics 9, 37–42 (2015). [CrossRef]

7. M. Koshiba, Y. Tsuji, and M. Hikari, “Time-domain beam propagation method and its application to photonic crystal circuits,” J. Lightwave Technol. 18(1), 102–110 (2000). [CrossRef]

8. V. F. R.- Esquerre, M. Koshiba, and H. E. H.- Figueroa, “Finite element analysis of photonic crystal cavities:time and frequency domain,” J. Lightwave Technol. 23(3), 1514–1521 (2005). [CrossRef]

9. T. Fujisawa and M. Koshiba, “Time-domain beam propagation method for nonlinear optical propagation analysis and its application to photonic crystal circuits,” J. Lightwave Technol. 22(2), 684–691 (2004). [CrossRef]

10. V. F. R.- Esquerre, M. Koshiba, and H. E. H.- Figueroa, “Frequency-dependent envelope finite element time domain analysis of dispersion materials,” Microwave and Opt. Technol. Lett. 44(1), 13–16 (2004). [CrossRef]

11. A. Niiyama, M. Koshiba, and Y. Tsuji, “An efficient scalar finite element formulation for nonlinear optical channel waveguides,” J. Lightwave Technol. 13(9), 1919–1925 (1995). [CrossRef]

12. G. R. Liu, “A generalized gradient smoothing technique and the smoothed bilinear form for galerkin formulation of a wide class of computational methods,” International Journal of Computational Methods 5(2), 199–236 (2008). [CrossRef]

13. G. R. Liu, Meshfree Methods: Moving Beyond the Finite Element Method (CRC Press, 2010).

14. R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems (Cambridge, 2002) [CrossRef]

15. G. R. Liu and N. T. Trung, Smoothed Finite Element Methods (CRC Press, 2010) [CrossRef]

16. P.-L. Liu, Q. Zhao, and F.-S. Choa, “Slow-wave finite-difference beam propagation method,” IEEE Photon. Technol. Lett. 7(8), 890–892 (1995). [CrossRef]

17. G. H. Jin, J. Harari, J. P. Vilcot, and D. Decoster, “An improved time domain beam propagation method for integrated optics components,” IEEE Photon. Technol. Lett. 9(3), 117–122 (1997). [CrossRef]

18. J. Lee and B. Fornberg, “A split step approach for the 3-D Maxwell’s equations,” J. Comput. Appl. Math. 158(2), 485–505 (2003). [CrossRef]

19. M. Movahhedi and A. Abdipour, “Alternating direction implicit formulation for the finite element time domain method,” IEEE Trans. Microwave Theory Technology 55(6), 1322–1331 (2007). [CrossRef]

20. J. F. Lee, “WETD-A finite element time-domain approach for solving Maxwell’s equations,” IEEE Microwave and Guided wave Letters 4(1), 11–13 (1994). [CrossRef]

21. S. S. A. Obayya, “Efficient finite-element-based time-domain beam propagation analysis of Optical integrated circuits,” IEEE J. Quantum Electron. 40(5), 591–595 (2004). [CrossRef]

22. V. F. R.- Esquerre and H. E. H.- Figueroa, “Novel time-domain step-by-step scheme for integrated optical applications,” IEEE Photon. Technol. Lett. 13(4), 311–313 (2001). [CrossRef]

23. G. R. Liu and M. B. Liu, Smoothed Particle Hydrodynamics - A Meshfree Method (World Scientific: Singapore, 2003).

24. H. A. Van der Vorst, “Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems,” SIAM J. Sci. and Stat. 13(2), 631–644 (1992).

25. J. K. Hwang, S. B. Hyun, H. Y. Ryu, and Yong-Hee Lee, “Resonant modes of two-dimensional photonic bandgap cavities determined by the finite-element method and by use of the anisotropic perfectly matched layer boundary condition,” J. Opt. Soc. Am. B 15(8), 2316–2324 (1998). [CrossRef]

26. I. Ahmed and E. Li, “Time domain simulation of dispersive materials from microwave to optical frequencies,” Proceedings of IEEE 7th International Conference on Emerging Technologies (ICET) (IEEE, 2011) 1–5.

27. S.J. Orfanidis, “Electromagnetic Waves and Antenna,” http://www.ece.rutgers.edu/orfanidi/ewa.

	perpendicular to y axis	perpendicular to z axis	PML corners
s_y	1	s	1
s_z	s	1	1

	SFETD-BPM-T3	FETD-BPM-T6	FETD-BPM-T3	Ref [25]
Number of elements	24604	24604	24604	—–
Number of unknowns	12383	49369	12383	—–
CPU time (s)	987	32239	1132	—–
F = a/λ	0.37877	0.37887	0.37940	0.37844
Error(F − F_o)/F_o %	0.0528%	0.0922%	0.253%	0 %

	perpendicular to y axis	perpendicular to z axis	PML corners
s_y	1	s	1
s_z	s	1	1

	SFETD-BPM-T3	FETD-BPM-T6	FETD-BPM-T3	Ref [25]
Number of elements	24604	24604	24604	—–
Number of unknowns	12383	49369	12383	—–
CPU time (s)	987	32239	1132	—–
F = a/λ	0.37877	0.37887	0.37940	0.37844
Error(F − F_o)/F_o %	0.0528%	0.0922%	0.253%	0 %

Efficient smoothed finite element time domain analysis for photonic devices

Abstract

1. Introduction

2. Analysis

2.1. Time domain equations

2.2. Galerkin weak form

2.3. Smoothed Galerkin procedure

2.4. Time domain discretization

3. Simulation results

3.1. Numerical dispersion

3.2. Optical grating

3.3. Photonic crystal cavity

3.4. Dispersive material

3.5. Enhanced cavity

4. Conclusion

References and links

Cited By

Figures (9)

Tables (3)

Equations (32)

Optics Express

	SFETD-BPM-T3	FETD-BPM-T3	FETD-BPM-T6
Number of elements	46240	46240	5780
Number of unknowns	23835	23835	23835
CPU time (s)	829	1481	2025