Even-Odd Integration Scheme

Overview: In the following, we introduce a 2-D explicit-implicit integration scheme that can be used to integrate the hyperbolic partial differential equations describing the propagation of the light fields inside a semiconductor waveguide cavity.

Finite Difference Even-Odd Integration Scheme
    1.1  Explicit Integration Step
    1.2  Implicit Integration Step
    1.3  Explicit-Implicit Integration Scheme
    1.4  Integraton Along The Boundary

The first step towards the numerical integration of a differential equation using a finite difference method involves the discretization of the differential operators. We illustrate this step using the hyperbolic partial differential equation:

⎡ ⎣ ∂ ∂t ± β k0 c ϵp ∂ ∂z ⎤ ⎦ ET± = ⎡ ⎣ i 2 ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ∂2 ∂x2 + i ω0 2 δϵa ϵp Γ− β k0 c ϵp α 2 ⎤ ⎦ ET± + i ω0 2 Γ ϵ0 ϵp PT± + F(t)ex,

(1.1)

describing the dynamics of the optical fields counter-propagating inside the wave-guide of a semiconductor amplifier or laser. A derivation of (1.1) starting from Maxwell's equations is presented here.



The differential operators occurring in (1.1) are discretized using central differences [2]:

∂ ∂t ET±     ⇒     ET±(t+1,i,j) − ET±(t−1,i,j) 2∆t , (1.2) ∂ ∂z ET±     ⇒     ET±(t,i+1,j) − ET±(t,i−1,j) 2∆z , (1.3) ∂2 ∂x2 ET±     ⇒     ET±(t,i,j+1) − 2ET±(t,i,j) + ET±(t,i,j−1) (∆x)2 , (1.4)

where t denotes the time point, and the indices i and j represent the position on the numerical grid in longitudinal and transverse direction, respectively (see Fig. 1.1). For the derivative in transverse direction, we have used the approximation: [(∂)/(∂x)] E_T^± ⇒ [(E_T^±(t,i,j+1/2) − E_T^±(t,i,j−1/2))/(∆x)], involving the field at half-steps between grid points. Nevertheless, the second derivative (1.4) is well defined in terms of integer grid number locations. Using the differential operators introduced in (1.2) - (1.4), the numerical approximation of (1.1) can be written as:

ET±(t+1,i,j) − ET±(t−1,i,j) 2∆t = ± β k0 c ϵp ET±(t−1,i+1,j) − ET±(t−1,i−1,j) 2∆z + ⎡ ⎣ i ω0 2 δϵa ϵp Γ− β k0 c ϵp α 2 ⎤ ⎦ ET±(t−1,i,j) + i 2 ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ET±(t−1,i,j+1) − 2ET±(t−1,i,j) + ET±(t−1,i,j−1) (∆x)2 + i ω0 2 Γ ϵ0 ϵp PT±(t−1,i,j)+ F(t−1,i,j)ex, (1.5)

Note: The integration alogrithm only converges if the following inequations are satisfied [2]:

β k0

c ϵp

∆t ∆z

≤ 1

and

1 ω0

c2 ϵp

∆t (∆x)2

≤ 1

(1.11)

From equation (1.5) we can determine E_T^±(t+1,i,j), the electric components of the optical field at time point t+1 and spatial coordinate (i∆z, j∆x):

ET±(t+1,i,j) = ⎧ ⎨ ⎩ 1 + ⎡ ⎣ i ω0 2 δϵa ϵp Γ− β k0 c ϵp α 2 ⎤ ⎦ 2∆t − i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ 2∆t (∆x)2 ⎫ ⎬ ⎭ ET±(t−1,i,j) ± β k0 c ϵp ∆t ∆z [ ET±(t−1,i+1,j) − ET±(t−1,i−1,j) ] + i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ∆t (∆x)2 [ ET±(t−1,i,j+1) + ET±(t−1,i,j−1)] + 2 ∆t ⎡ ⎣ i ω0 2 Γ ϵ0 ϵp PT±(t−1,i,j)+ F(t−1,i,j)ex ⎤ ⎦ . (1.6)

(1.6) is an explicit equation, since the quantities on the right side depend on the previous time step t−1. These variables are known either from the initial boundary conditions or from the previous integration step.

Alternatively, the spatial derivatives introduced in (1.3) and (1.4) can be evaluated at time step t + 1. In this case, the numerical approximation of wave-equation (1.1) is given by:

ET±(t+1,i,j) − ET±(t−1,i,j) 2∆t = ± β k0 c ϵp ET±(t+1,i+1,j) − ET±(t+1,i−1,j) 2∆z + ⎡ ⎣ i ω0 2 δϵa ϵp Γ− β k0 c ϵp α 2 ⎤ ⎦ ET±(t+1,i,j) + i 2 ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ET±(t+1,i,j+1) − 2ET±(t+1,i,j) + ET±(t+1,i,j−1) (∆x)2 + i ω0 2 Γ ϵ0 ϵp PT±(t+1,i,j)+ F(t+1,i,j)ex, (1.7)

From (1.7), the electric component of the optical field at time point t+1 and spatial coordinate (i∆z,j ∆x) can be calculated using:

ET±(t+1,i,j) = ⎧ ⎨ ⎩ ET±(t−1,i,j)± β k0 c ϵp ∆t ∆z [ ET±(t+1,i+1,j) − ET±(t+1,i−1,j) ] + i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ∆t (∆x)2 [ ET±(t+1,i,j+1) + ET±(t+1,i,j−1)] + i ω0 Γ∆t ϵ0 ϵp PT±(t+1,i,j) + 2 ∆t[ F(t+1,i,j)ex] ⎫ ⎬ ⎭ ⎧ ⎨ ⎩ 1 − ⎡ ⎣ i ω0 2 δϵa ϵp Γ− β k0 c ϵp α 2 ⎤ ⎦ 2∆t + i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ 2∆t (∆x)2 ⎫ ⎬ ⎭ −1   . (1.8)

In this case, the right side of (1.8) contains terms that refer to the previous time step: t−1, and the new time step: t+1. In order to determine E_T^±(t+1,i,j) at all grid points, we would have to solve a system of linear equations during each iteration step. An integration algorithm of this type is called `implicit'. It requires additional computational effort, but adds numerical stability to the integration scheme [2].

A simplified implicit integration algorithm, established by Gourlay [1], can be formulated by observing that only terms containing the electric field at the nearest neighbouring grid points appear in (1.8). Marking the points of the numerical grid in the style of a checker board (see Fig. 1.1), allows us to distinguish between `even' and `odd' points.

Figure 1.1: Numerical integration grid. The grid-points marked with black and white fields refer to the even-odd numerical integration algorithm.

Points on black fields are called `even', since the sum of their indices is an even number. Points located on white fields are denoted `odd', since the sum of their indices is an odd number. Formally, these definitions are given by:

(i,j)even = { (i,j) | i+j=2k, k∈N, i ∈{1,…,N}, j ∈{1,…,M} }, (1.9) (i,j)odd = { (i,j) | i+j=2k+1, k∈N, i ∈{1,…,N}, j ∈{1,…,M} } (1.10)

where N and M are the number of numerical grid points in longitudinal and transverse direction, respectively, and N represents the set of natural numbers.

Using (1.8), the optical field at even grid positions can the determined, if the field at odd grid positions (at the same time step) is known. The required field quantities can be calculated using the explicit equation (1.6). In the next integration step, the variables at odd grid positions are calculated first [using the explicit integration equation (1.6)]. Subsequently, the variables at even position can be calculated using (1.8). The integration scheme is visualized in Fig. 1.2

Figure 1.2: Even-odd explicit-implicit integration scheme. The implicit integration step always follows the explicit integration step. The cyan coloured fields on the checker-board indicate the coordinates of the field quantities being updated in the current step.

The algorithm shown in Fig. 1.2 can be used to integrate the optical field at inner grid points. In order to integrate along the longitudinal boundary, we have to take into consideration that optical fields propagating in positive direction are partially reflected at the rear laser facet. Similarly, fields propagating in negative direction are partially reflected at the front facet.

For the integration along the front facet, we introduce an additional term in order to simulate the injection of an optical pulse:

ET+(t+1,1,j) = (1−R1)  Π(t,j) e−i(ωopinj − ωop) − R1 ET−(t,1,j). (1.12)

The term Π(t,j) is defined by the shape and duration of the pulse, ω_op^inj is the central frequency of the optical pulse, and R₁ is the reflectivity of the front facet with respect to the optical field. For the injection of a Gaussian shaped pulse with a duration of ∆t_p, a half-width in transverse direction of ∆w_p, and a maximum intensity of I₀, Π(t,j) is given by:

Π(t,j) = √   I0   exp ⎡ ⎢ ⎣ − 1 2 t − t0 ∆tp (2 √   2 ln2   ) ⎤ ⎥ ⎦ exp ⎡ ⎢ ⎣ − 1 2 j∆x − x0 ∆wp (2 √   2 ln2   ) ⎤ ⎥ ⎦ , (1.13)

where t₀ is the injection time-point and x₀ is the central injection position in transverse direction. The equation for the integration of E_T⁻ along the front facet is given by:

ET−(t+1,1,j) = ⎧ ⎨ ⎩ 1 + ⎡ ⎣ i ω0 2 δϵa ϵp Γ− β k0 c ϵp α 2 ⎤ ⎦ 2∆t − i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ 2∆t (∆x)2 ⎫ ⎬ ⎭ ET−(t−1,1,j) + β k0 c ϵp 2∆t ∆z [ ET−(t−1,2,j) − ET−(t−1,1,j) ] + i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ∆t (∆x)2 [ ET−(t−1,1,j+1) + ET−(t−1,1,j−1)] + 2 ∆t ⎡ ⎣ i ω0 2 Γ ϵ0 ϵp PT−(t−1,1,j)+ F(t−1,1,j)ex ⎤ ⎦ . (1.14)

For the integration of the optical field along the output facet we use the equations:

ET−(t+1,N,j) = − R2 ET+(t,N,j), (1.15)

ET+(t+1,N,j) = ⎧ ⎨ ⎩ 1 + ⎡ ⎣ i ω0 2 δϵa ϵp Γ − β k0 c ϵp α 2 ⎤ ⎦ 2∆t − i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ 2∆t (∆x)2 ⎫ ⎬ ⎭ ET+(t−1,N,j) + β k0 c ϵp 2∆t ∆z [ ET+(t−1,N,j) − ET+(t−1,N−1,j) ] + i ω0 ⎛ ⎝ c2 ϵp ⎞ ⎠ ∆t (∆x)2 [ ET+(t−1,N,j+1) + ET+(t−1,N,j−1) ] + 2 ∆t ⎡ ⎣ i ω0 2 Γ ϵ0 ϵp PT+(t−1,N,j)+ F(t−1,N,j)ex ⎤ ⎦ , (1.16)

where i ∈ {1,2,…,N} and N is the number of numerical grid points in longitudinal direction. The integration of the optical fields along the longitudinal boundaries has been introduced using the explicit equation (1.6). The implicit equation (1.8) can be modified in a similar way to perform the integration along the longitudinal boundaries during the implicit cycle of the even-odd integration algorithm. The transverse boundary conditions are non-critical, if we extend the numerical grid (in transverse direction) well beyond the region of the injection stripe. In this case, the optical fields are concentrated along a central axis in propagation direction, while regions on the left and right side of the injection stripe form the `lossy wings', absorptive regions with a low 2-D carrier density. To integrate along the transverse boundary, one can use von Neumann boundary conditions, specifying the spatial derivative of the variable. Alternatively, we can use Dirichlet boundary conditions, specifying the value of the variable at the boundary.

[1]

A.R. Gourlay and J. L. Morris. Hopscotch difference methods for nonlinear hyperbolic systems. IBM Journal of Research and Development, 16(4):349, 1972.

[2]

Franz J. Vesely. Computational Physics: An Introduction. Kluver Academic Publishers, New York, 2 edition, 2001.