EvenOdd Integration Scheme
Overview: In the following, we introduce a 2D explicitimplicit
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.
Contents
Finite Difference EvenOdd Integration Scheme
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

± 
β
k_{0}


c
ϵ_{p}


∂
∂z
 ⎤ ⎦

E_{T}^{±} =  ⎡ ⎣

i
2 ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


∂^{2}
∂x^{2}

+ 
i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ− 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

E_{T}^{±} + 
i ω_{0}
2


Γ
ϵ_{0} ϵ_{p}

P_{T}^{±} + F(t)e_{x}, 
  (1.1) 
describing the dynamics of the
optical fields counterpropagating inside the waveguide 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]:
 


E_{T}^{±}(t+1,i,j) − E_{T}^{±}(t−1,i,j)
2∆t

, 
  (1.2) 
 


E_{T}^{±}(t,i+1,j) − E_{T}^{±}(t,i−1,j)
2∆z

, 
  (1.3) 
 


E_{T}^{±}(t,i,j+1) − 2E_{T}^{±}(t,i,j) + E_{T}^{±}(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 halfsteps 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:


E_{T}^{±}(t+1,i,j) − E_{T}^{±}(t−1,i,j)
2∆t

= 
 

± 
β
k_{0}


c
ϵ_{p}


E_{T}^{±}(t−1,i+1,j) − E_{T}^{±}(t−1,i−1,j)
2∆z

+  ⎡ ⎣

i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ− 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

E_{T}^{±}(t−1,i,j) 
 

+ 
i
2 ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


E_{T}^{±}(t−1,i,j+1) − 2E_{T}^{±}(t−1,i,j) + E_{T}^{±}(t−1,i,j−1)
(∆x)^{2}


 
 + 
i ω_{0}
2


Γ
ϵ_{0} ϵ_{p}

P_{T}^{±}(t−1,i,j)+ F(t−1,i,j)e_{x}, 
  (1.5) 

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

β
k_{0}


c
ϵ_{p}


∆t
∆z
 ≤ 1 
and 
1
ω_{0}


c^{2}
ϵ_{p}


∆t
(∆x)^{2}
 ≤ 1 

(1.11) 
1.1 Explicit Integration Step
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):

E_{T}^{±}(t+1,i,j) =  ⎧ ⎨
⎩

1 +  ⎡ ⎣

i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ− 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

2∆t − 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


2∆t
(∆x)^{2}

 ⎫ ⎬
⎭

E_{T}^{±}(t−1,i,j) 
 

± 
β
k_{0}


c
ϵ_{p}


∆t
∆z

[ E_{T}^{±}(t−1,i+1,j) − E_{T}^{±}(t−1,i−1,j) ] 
 

+ 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


∆t
(∆x)^{2}

[ E_{T}^{±}(t−1,i,j+1) + E_{T}^{±}(t−1,i,j−1)] 
 
 + 2 ∆t  ⎡ ⎣

i ω_{0}
2


Γ
ϵ_{0} ϵ_{p}

P_{T}^{±}(t−1,i,j)+ F(t−1,i,j)e_{x}  ⎤ ⎦

. 
  (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.
1.2 Implicit 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 waveequation (
1.1)
is given by:


E_{T}^{±}(t+1,i,j) − E_{T}^{±}(t−1,i,j)
2∆t

= 
 

± 
β
k_{0}


c
ϵ_{p}


E_{T}^{±}(t+1,i+1,j) − E_{T}^{±}(t+1,i−1,j)
2∆z

+  ⎡ ⎣

i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ− 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

E_{T}^{±}(t+1,i,j) 
 

+ 
i
2 ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


E_{T}^{±}(t+1,i,j+1) − 2E_{T}^{±}(t+1,i,j) + E_{T}^{±}(t+1,i,j−1)
(∆x)^{2}


 
 + 
i ω_{0}
2


Γ
ϵ_{0} ϵ_{p}

P_{T}^{±}(t+1,i,j)+ F(t+1,i,j)e_{x}, 
  (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:

E_{T}^{±}(t+1,i,j) =  ⎧ ⎨
⎩

E_{T}^{±}(t−1,i,j)± 
β
k_{0}


c
ϵ_{p}


∆t
∆z

[ E_{T}^{±}(t+1,i+1,j) − E_{T}^{±}(t+1,i−1,j) ] 
 

+ 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


∆t
(∆x)^{2}

[ E_{T}^{±}(t+1,i,j+1) + E_{T}^{±}(t+1,i,j−1)] + 
i ω_{0} Γ∆t
ϵ_{0} ϵ_{p}

P_{T}^{±}(t+1,i,j) 
 
 + 2 ∆t[ F(t+1,i,j)e_{x}]  ⎫ ⎬
⎭

 ⎧ ⎨
⎩

1 −  ⎡ ⎣

i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ− 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

2∆t + 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{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].
1.3 ExplicitImplicit Integration Scheme
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 gridpoints marked with black and white fields refer to the evenodd 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)  i+j=2k, k∈N, i ∈{1,…,N}, j ∈{1,…,M} },

  (1.9) 
 

{ (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: Evenodd explicitimplicit integration scheme. The implicit integration step always follows the explicit integration step.
The cyan coloured fields on the checkerboard indicate the coordinates of the field quantities being updated in the current step.
1.4 Integration Along The Boundary
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:
 E_{T}^{+}(t+1,1,j) = (1−R_{1}) Π(t,j) e^{−i(ωopinj − ωop)} − R_{1} E_{T}^{−}(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_{1} 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 halfwidth
in transverse direction of ∆
w_{p}, and a maximum intensity of
I_{0}, Π(
t,
j) is given by:
 Π(t,j) = 
√

I_{0}

exp  ⎡ ⎢
⎣

− 
1
2


t − t_{0}
 ⎤ ⎥
⎦

exp  ⎡ ⎢
⎣

− 
1
2


j∆x − x_{0}
 ⎤ ⎥
⎦

, 
  (1.13) 

where
t_{0} is the injection timepoint and
x_{0} is the central injection position in transverse direction.
The equation for the integration of
E_{T}^{−} along the front facet is given by:

E_{T}^{−}(t+1,1,j) =  ⎧ ⎨
⎩

1 +  ⎡ ⎣

i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ− 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

2∆t − 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


2∆t
(∆x)^{2}

 ⎫ ⎬
⎭

E_{T}^{−}(t−1,1,j) 
 

+ 
β
k_{0}


c
ϵ_{p}


2∆t
∆z

[ E_{T}^{−}(t−1,2,j) −
E_{T}^{−}(t−1,1,j) ] 
 

+ 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


∆t
(∆x)^{2}

[ E_{T}^{−}(t−1,1,j+1) +
E_{T}^{−}(t−1,1,j−1)] 
 
 +
2 ∆t  ⎡ ⎣

i ω_{0}
2


Γ
ϵ_{0} ϵ_{p}

P_{T}^{−}(t−1,1,j)+ F(t−1,1,j)e_{x}  ⎤ ⎦

. 
  (1.14) 

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

E_{T}^{−}(t+1,N,j) = −
R_{2} E_{T}^{+}(t,N,j), 

 (1.15) 


E_{T}^{+}(t+1,N,j) =  ⎧ ⎨
⎩

1 +  ⎡ ⎣

i ω_{0}
2


δϵ_{a}
ϵ_{p}

Γ − 
β
k_{0}


c
ϵ_{p}


α
2
 ⎤ ⎦

2∆t − 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


2∆t
(∆x)^{2}

 ⎫ ⎬
⎭

E_{T}^{+}(t−1,N,j) 
 

+ 
β
k_{0}


c
ϵ_{p}


2∆t
∆z

[ E_{T}^{+}(t−1,N,j) − E_{T}^{+}(t−1,N−1,j) ] 
 

+ 
i
ω_{0}

 ⎛ ⎝

c^{2}
ϵ_{p}
 ⎞ ⎠


∆t
(∆x)^{2}

[ E_{T}^{+}(t−1,N,j+1) + E_{T}^{+}(t−1,N,j−1) ] 
 
 + 2 ∆t  ⎡ ⎣

i ω_{0}
2


Γ
ϵ_{0} ϵ_{p}

P_{T}^{+}(t−1,N,j)+ F(t−1,N,j)e_{x}  ⎤ ⎦

, 
  (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 evenodd integration algorithm.
The transverse boundary conditions are noncritical, 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 2D 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.
Bibliography
 [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.