## Abstract

We report a general computational model of complex material media for electrodynamics simulation using the Finite-Difference Time-Domain (FDTD) method. It is based on a multi-level multi-electron quantum system with electron dynamics governed by Pauli Exclusion Principle, state filling, and dynamical Fermi-Dirac Thermalization, enabling it to treat various solid-state, molecular, or atomic media. The formulation is valid at near or far off resonance as well as at high intensity. We show its FDTD application to a semiconductor in which the carriers’ intraband and interband dynamics, energy band filling, and thermal processes were all incorporated for the first time. The FDTD model is sufficiently complex and yet computationally efficient, enabling it to simulate nanophotonic devices with complex electromagnetic structures requiring simultaneous solution of the mediumfield dynamics in space and time. Applications to direct-gap semiconductors, ultrafast optical phenomena, and multimode microdisk lasers are illustrated.

© 2006 Optical Society of America

## 1. Introduction

Numerical computation for electrodynamics capable of providing full spatial-temporal solution for the electromagnetic field interacting with material media using Finite Difference Time Domain (FDTD) [Refs. 1, 2] method has attracted much interest in photonic and optoelectronic device simulation, especially in situations where conventional device simulation methods are inadequate or cannot cope with the complex device geometry. An example is the case of photonic bandgap or micocavity structures with active semiconductor medium as part of the structure to provide gain or ultrafast optical nonlinearities. In this case, the spatial variations of the medium’s refractive index, gain, or absorption could be dependent transiently on the varying lightwave intensity while the mode profile that determines the intensity could in turn be affected by the medium’s refractive index, gain, absorption, and the optical nano-structure. Another example is the case of a laser cavity with multiple longitudinal or transverse modes in which multi-frequency lasing, instability, or laser pulsing behaviors could develop. These behaviors could be dependent on the carrier band-filling effect and the spatial or spectral hole burning of the carrier distribution, which in turn could be affected by the transient and local interaction of the electromagnetic field mode with the gain medium.

In actual implementation, the dynamics of the active media in the devices must be modeled with reasonable physical realism in order for the potentials of the FDTD method to be fully realized. The medium involved typically has complicated internal electron dynamics. For example, the dynamics of multiple electrons could be involved and the quantum energy level structure for the electrons in the medium must be accounted for. Furthermore, these devices often operate at room temperature, and the thermal process of the electrons must be considered. It would be of substantial interest to develop a FDTD computational model for complex dynamical media that is sophisticated enough to encompass the essential physics of such media and yet computationally efficient. We will show that with use of a multi-level multi-electron quantum system having the electron dynamics governed by Pauli Exclusion Principle, state filling, and dynamical Fermi-Dirac Thermalization, the essential physics of many complex media, including semiconductor, can be encompassed in FDTD simulation. We will refer to this FDTD computational model as dynamical-thermal-electron quantum-medium FDTD (DTEQM-FDTD). The dynamical-thermal-electron nature of the model refers to the explicit inclusion of thermal processes dynamically for the electrons. The quantum nature of the model refers to the explicit inclusion of energy level structure and the Pauli Exclusion Principle. We will show how a simple thermal hopping model dictated by the Pauli Exclusion Principle automatically resulted in Fermi Dirac Statistics for the electron state filling in the steady state. The broad applicability of the model methodology relies on the fact that the interaction of electromagnetic field with material media is via a collection of oscillating electric dipoles governed by the electron dynamics and the DTEQM-FDTD model shows how the complex electron dynamics in multiple energy levels can be efficiently treated.

The model methodology shall be generally applicable to a wide range of solid-state, molecular, or atomic media, by employing the treatment of electron dynamics described here with the appropriate medium-specific effective carrier rate equations and effective multi-energy-level structure needed for obtaining the medium polarization density. One of the most complex material media for photonic devices is direct bandgap semiconductor. We will illustrate the powerful capability of the model methodology by specifically developing it for modeling semiconductor medium as an example to show how the effective carrier rate equations and multi-energy-level structures are formulated, making them computationally efficient. Its applications to other solid-state, molecular, or atomic media will involve similar approaches and in many cases will be simpler than the semiconductor example.

In applying to semiconductor, we show that this model could take into account the transient intraband and interband electron dynamics, the semiconductor band structure, and the carrier thermal equilibrium process. The model automatically incorporates the required band filling effect. It also incorporates the typical nonlinear optical effects associated with carrier dynamics. Our approach enables the model to treat thermally activated carrier scattering process under transient excitation spatial-temporally. These capabilities enable the model to treat sophisticated semiconductor optoelectronic and nanophotonic devices having complex device geometries and medium dynamics with full spatial-temporal solutions at the microscopic level. We will show that although the model is complex enough to include all the essential physical effects, it is simple enough to achieve fast computational time.

In this paper, we present the basic formulation of this powerful general model with applications to direct bandgap semiconductor, ultrafast optical phenomena, and multimode microdisk semiconductor laser as illustrations, which also validate the computational stability and efficiency of the model.

## 2. General approach

In the past, various methods have been proposed to model active medium in FDTD simulation. They can be grouped into two main categories: the macroscopic approach [Ref. 3] and the microscopic approach [Ref. 4]. The macroscopic approach models the medium gain by introducing a negative frequency-dependent conductivity term into Maxwell equations and solve for the imaginary part of the polarization, while the microscopic approach explores rate equations. The rate equation approach in [Ref. 4] employed a modified classical electron model to describe the gain and absorption dynamics of an atomic system with a single electron. Recently, we have introduced a quantum mechanical model of a 4-level 2-electron atomic system with the incorporation of the Pauli Exclusion Principle for FDTD simulation [Ref. 5–8]. The employment of two electrons enabled us to provide a simplified model for semiconductor, which took into account a simple picture of electron-hole pumping dynamics from the lower valence band to the upper conduction band. Our quantum-medium model [Refs. 5–8] with four energy levels and two electrons governed by the Pauli Exclusion Principle is an advancement over the previous models. However, the model is still too simple to properly encompass the complex physical effects in the medium such as the semiconductor energy band structure, band-filling effect with Fermi Dirac statistics, carrier induced gain and refractive index change, and carrier relaxation to thermal equilibrium after excitation. As a result, much of the transient and nonlinear behaviors in the medium were not included.

The typical approach to modeling carriers in semiconductor band structure involves solving Bloch equations at many energy states in the momentum space (k-states) [Ref. 9, 10]. In FDTD simulation, the structure to be simulated is first discretized spatially, and then the electromagnetic field at each spatial point is updated at each time step, making FDTD an intrinsically numerically intensive method. Now in addition to all that, if we use the typical approach of semiconductor modeling, then for each grid, we will have to update the carrier distribution function in many k-states at each time step [Ref. 11], making the computational time forbiddingly long. The original k-states are quantized states given by k_{x}=2*π*n_{x}/L_{x}, which is dependent on the size of the medium L_{x}. As these k-states are energy broadened by the spontaneous decay time, we may combine the original k-states into energy-broadened k-states that are interspaced by the spontaneous decay width. With a typical spontaneous decay time of ~ 3ns for direct bandgap III-V semiconductors, there will be ~10^{5} energy broadened k-states within a typical 50nm wavelength spectral width of interest at the optical communication wavelength of *λ*=1550nm This is still a large number for FDTD simulation.

In our 4-level 2-electron FDTD model treatment [Ref. 5–8], it is pointed out that due to the short dipole de-phasing time of excited dipoles in semiconductor (~100fsec), comparing to the spontaneous decay lifetime of 3nsec, one can further represent many energy-broadened k states by one effective k state, provided the electron-scattering process that affect the dipole phase is lumped into an effective dipole dephasing time. For, example, a dipole dephasing time (T_{2}) of ~50fsec will result in a spectral broadening of Δ*λ*~50nm (FWHM) at *λ*=1550nm. Thus, only a few such dipoles are needed to cover a wavelength range of 100-300nm, which is sufficient for various optoelectronic and photonic device applications.

It turns out that the dipole dephasing is due mainly to the thermally activated carrier scattering. At room temperature, the thermal energy of *k*_{B}*T* (*k*_{B}
is the Boltzman constant and *T* is the Kelvin temperature) corresponds to a FWHM energy broadening of *ħ*Δ*ω* =18meV (Δ*λ*=35nm at *λ*=1550nm). For a solid-state medium with strong electron interaction, the dipole dephasing may correspond to the thermal energy, giving a thermally-induced dipole dephasing time of T_{2}=2/Δ*ω*=75fsec at 300K. Thus, although the effective values for T_{2} in different media may depend on the actual dephasing mechanisms (the T_{2} value can be experimentally obtained from the homogeneously-broadened linewidth), they will be in the order of ~100fsec [Ref. 9, 10]. Hence, the fast dipole dephasing is generally applicable to all solid-state media with strong electron interactions. This means that most solid-state media can also be modeled with only a few dipoles to cover a wide wavelength range. While illustrated for the case of semiconductor, the DTEQM-FDTD modeling approach developed here can be adapted for modeling a wide range of solid-state, molecular, and atomic media.

To set up the model, let us first consider a simple case for which a semiconductor medium is under optical excitation with above-bandgap energy, as shown in Fig. 1(a). The medium is at ground state initially. Under excitation of the photons, the electron at the corresponding energy (position 1) will undergo interband excitation to the conduction band (position 2), leaving a hole at the valence band. The electron at position 2 and the hole at position 1 in Fig. 1(a) will then undergo intraband decay to the band edge positions 3 and 0, respectively, through phonon-assisted processes. The hole position moving from 1 to 0 in the valence band is equivalent to electron moving from position 0 to 1. Subsequently, the electron and hole will recombine via radiative or nonradiative decay and the medium will return to its original ground state. Current injection via p-i-n junction can likewise be modeled by pumping electrons electrically from the lowest valence-band level to the topmost conduction-band level. In this example, we can see that the broad dynamics for the electron transition can be described by having two electrons resting in four energy levels, as shown in Fig. 1(b). To ensure that the electron will not decay to the lower level unless there is an energy-state vacancy, Pauli Exclusion term was implemented. This forms the basis of our earlier 4-Level 2-Electron Model with Pauli Exclusion Principle [Refs. 5–8].

To more accurately model semiconductor medium dynamics, additional energy levels and electron dynamics are needed. The main idea of our multi-level multi-electron model discussed recently [Ref. 12] is to divide the conduction and valence band states to several groups, and then represent each group by a single dipole with broadened width, as shown in Fig. 2.

By doing so, the semiconductor medium dynamics included are: interband carrier radiative or nonradiative decays, intraband carrier relaxation, and interband transition (gain or absorption). The semiconductor density of states is accounted for in the multiple energy level structure. Implementation of Pauli Exclusion Principle for the transitions between any two levels will result in carrier band-filling effect but would not account for finite temperature.

In our model, finite temperature is accounted for by a temperature-dependent intraband carrier-hopping dynamics that ensures Fermi-Dirac thermal distribution for the electron/hole occupation in the steady state, which we refer to as dynamical Fermi-Dirac thermalization. As will be seen below, this multi-level multi-electron model gives the expected band-filling behavior for a direct bandgap semiconductor. In a more sophisticated situation not shown here, by using temperature as a spatial-temporally varying parameter, this approach will enable the FDTD model to treat thermally activated carrier scattering process and carrier heating or cooling under transient excitation in both space and time, while solving for the interaction of the electromagnetic field modes with the medium. The spatial movement of carriers is also not shown here but can be built on top of the FDTD model [Refs. 9, 10].

Thus, the basic DTEQM-FDTD model developed here, when combined with a lattice temperature diffusion model and/or a spatial carrier diffusion model, will result in a powerful device model with sophisticated spatial-temporal solutions for both the medium and electromagnetic field. When applied to semiconductor, we will simply refer to the DTEQM-FDTD model as “dynamical-semiconductor-medium FDTD” (DSM-FDTD) model.

## 3. The dynamical-semiconductor-medium (DSM) FDTD model

The energy-level structure for the DSM-FDTD model considered below is shown in Fig. 3. The continuous semiconductor band structure is simplified to discrete levels. The conduction band levels are labeled as *i*_*c* and the valence band levels are labeled as *i*_*ν*, where *i*=1,2, *M* and *M* is the total number of levels used in the model. *E*
_{i} represents the transition energy between levels *i*_*c* and *i*_*ν* and obeys *E*
_{i} =*ħω*_{i}
=2*πħc*/*λ*_{i}
. Levels *i*_*c* and *i*_*ν* are used to represent the conduction and valence band states, respectively, with optical transition energy between ${E}_{i-}^{B}$ = *E*
_{i} -(*E*
_{i} - *E*
_{i-1})/2 and ${E}_{i+}^{B}$ = *E*
_{i} +(*E*
_{i+1} - *E*
_{i})/2 , so the densities of states per
unit volume (volume densities of states) ${N}_{\mathit{\text{Ci}}}^{0}$ and ${N}_{\mathit{\text{Vi}}}^{0}$ in *i*_*c* and *i*_*ν*, respectively, will be the sum of all the volume densities of states with energy in the bracketed interval [${E}_{i-}^{B}$, ${E}_{i+}^{B}$] as ${E}_{(c\mathit{,}v)i}^{0}$ = *∫*
^{${E}_{i+}^{\mathbf{B}}$}
_{${E}_{i-}^{B}$}
*ρ*
_{c,ν}(*E*)*dE* , where *ρ*_{c,ν}
(*E*) is the number of states per unit volume per unit energy at the optical transition energy *E* for conduction band (C) and valance band (V).

If we take level *i*_*c* as example, the semiconductor carrier dynamics we considered in the model include: (1) interband radiative and nonradiative decays between conduction band levels *i*_*c* and *i*_*ν* with decay time *τ*
_{i}; (2) interband transition (stimulated decay or absorption) governed by dipole matrix element; (3) intraband down-transition (phonon-assisted decay) from *i*_*c* to *i*-1_*c* with transition time *τ*
_{(i,i-1)C}; (4) intraband up-transition (thermally activated) from *i*-1_*c* to *i*_*c* with transition time *τ*
_{(i-1,i)C}; (5) Intraband down-transition from *i*+1_*c* to *i*_*c* with transition time *τ*
_{(i+1,i)C}; (6) intraband up-transition (thermally activated) from *i*_*c* to *i*+1_*c* with transition time *τ*
_{(i,i+1)C}. Note that intraband transitions are between adjacent levels with times labeled by *τ*
_{(i,f)C,V}, where indices *i* and *f* label the initial and final levels, respectively.

## 4. Quantum derivation of basic set of equations with minimal approximation

To obtain the basic set of equations governing the multi-level system quantum mechanically, let us first consider the optical interaction that occurred between interband levels. All the carrier dynamics below are described by electrons, including the hole carrier dynamics in the valence band. However, the electrons in the conduction and valence bands are characterized by different effective masses and intraband transition rates. The equations governing the transitions can be derived as follows.

We first start with the following exact minimal-coupling non-relativistic Hamiltonian describing atom-field interactions for a collection of *N*
_{0} single-electron atoms [Ref. 13]:

$$+\frac{1}{2}{\int}_{{V}_{Q}}{d}^{3}\mathbf{r}\sum _{\mathit{m\sigma}}\frac{1}{2}\left[\frac{d\left({\Omega}_{m}{\epsilon}_{m}\right)}{d{\Omega}_{m}}{\hat{\mathbf{E}}}_{\mathit{m\sigma}}^{2}(\mathbf{r},t)+{\mu}_{0}{\hat{\mathbf{H}}}_{\mathit{m\sigma}}^{2}\left(\mathbf{r},t\right)\right]+{\hat{H}}_{\mathrm{ph},}$$

where **r̂**_{ej} is the position operator of the electron in the *j*
^{th} atom, **p̂**_{ej} is the canonical momentum operator of the electron in the *j*
^{th} atom, *m*
_{e} is the electron mass, *e* = -|*e*| is the electron charge. The second term $\sum _{j=1}^{{N}_{0}}}{\hat{H}}_{\mathit{Cj}$is the Coulomb interaction energy. The last term *Ĥ*_{ph} is a dipole dephasing reservoir Hamiltonian similar to that given in [Ref. 13]. The dynamics of the heavy nucleus has been neglected. **Ê**_{mσ} and **Ĥ**_{mσ} are the mode amplitudes of the electric and magnetic field operators and the field operators are given by:

$$\hat{\mathbf{E}}({\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right),t)=\sum _{\mathit{m\sigma}}i{\Omega}_{m}{g}_{m}{\mathbf{e}}_{\mathit{m\sigma}}\left[{\hat{a}}_{\mathit{m\sigma}}\left(t\right){e}^{i{\mathbf{k}}_{m}\xb7{\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right)}-{\hat{a}}_{\mathit{m\sigma}}^{\u2020}\left(t\right){e}^{-i{\mathbf{k}}_{m}\xb7{\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right)}\right]=\sum _{\mathit{m\sigma}}{\hat{\mathbf{E}}}_{\mathit{m\sigma}}({\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right),t),$$

$${\mu}_{0}\hat{\mathbf{H}}({\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right),t)=\sum _{\mathit{m\sigma}}i{g}_{m}\left({\mathbf{K}}_{\mathit{m\sigma}}\times {\mathbf{e}}_{\mathit{m\sigma}}\right)\left[{\hat{a}}_{\mathit{m\sigma}}\left(t\right){e}^{i{\mathbf{k}}_{m}\xb7{\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right)}+{\hat{a}}_{\mathit{m\sigma}}^{\u2020}\left(t\right){e}^{-i{\mathbf{k}}_{m}\xb7{\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right)}\right]=\sum _{\mathit{m\sigma}}{\hat{\mathbf{H}}}_{\mathit{m\sigma}}({\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right),t),$$

where *g*
_{m} = [*ħ*
*μ*
_{m}/(2*ε*
_{0}
*V*_{Q}
*n*_{m}
Ω_{m}
*c*)]^{1/2} in which *n*_{m}
= ${\epsilon}_{m}^{1/2}$ is the refractive index of the embedded medium at frequency Ω_{m}, Ω_{m} = |*k*_{mσ}
|*c*/*n*_{m}
is the angular frequency of *m*
^{th} field mode, and *ν*_{m}
= *∂*Q_{m}/*∂*
*k*_{m}
is the group velocity. *V*_{Q}
is the volume of quantization. The constant *g*
_{m} given is valid for quantized field operators in a dispersive dielectric medium [Ref. 14]. The operator *â*_{mσ}(*t*) and ${\widehat{a}}_{m\sigma}^{\u2020}$(*t*) are the annihilation and creation operators of the *m*
^{th} field mode with polarization **σ** (**σ**=1, 2), and **e**
_{mσ} is the unit vector representing the **σ** polarization of mode **m**. The mode number **m** = {*m*_{x}
, *m*_{y}
, *m*_{z}
} can be indexed by its *m*_{x}
, *m*_{y}
, *m*_{z}
components with **k**
_{m} = *k*
_{mx}
**e**
_{x} + *k*
_{my}
**e**
_{y} + *k*
_{mz}
**e**
_{z} and |*k*_{mx}
| = 2*πm*_{x}
/*L*_{x}
, |*k*_{my}
| = 2*πm*_{y}
/*L*_{y}
, |*k*_{mz}
| = 2*πm*_{z}
/*L*_{z}
, so that *L*_{x}
, *L*_{y}
, *L*_{z}
are the *x, y,* and *z* dimensions for an arbitrarily large volume of quantization *V*_{Q}
(*V*_{Q}
= *L*_{x}
*L*_{y}
*L*_{z}
). The mode numbers take on positive and negative integer values as *m*_{x,y,z}
∈ {0, ±1, ±2, ±3, …}. Using Heisenberg equation of motion, we obtain:

$$\frac{de{\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right)}{dt}=\frac{i}{\mathit{\u0127}}[H\left(t\right),e{\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right)]=\frac{e}{{m}_{e}}\left[{\hat{\mathbf{P}}}_{\mathit{ej}}\left(t\right)-e\hat{\mathbf{A}}({\hat{\mathbf{r}}}_{\mathit{ej}}\left(t\right),t)\right]=\frac{d{\hat{\mathbf{\mu}}}_{j}(t)}{dt},$$

where **μ̂**_{j} (*t*) ≡ *e*(**r̂**_{ej} (*t*) - **r**
_{nj}) is the atomic dipole moment operator (**r**
_{nj} is the nucleus position). This gives:

To show how to incorporate three-dimensional (3D) polarization vector, we consider a simple case where a two-level system has three upper levels |*E _{us}* 〉

_{j}with

*s*∊ {

*x, y, z*} and one lower level |

*E*

_{g}〈

_{j}.

*n̂*

_{gj}(

*t*) and

*n̂*

*(*

_{usj}*t*) are the atomic number operators in ground (

*g*) level and upper (

*u*

_{s}) levels, respectively [Ref. 15]. We have the atomic energy down transition operator

*V̂*

*= |*

_{gusj}*E*

_{g}〉

_{jj}〈

*E*|, the atomic energy up transition operator

_{us}*V̂*

^{†}

*= |*

_{gusj}*E*〉

_{us}_{jj}〈

*E*

_{g}|, and the dipole matrix element

**μ**

*≡*

_{gusj}_{j}〈

*u*

_{s}|

**μ̂**

_{ej}|

*g*〉

_{j}. For simplicity, we will represent

*u*

_{s}by

*s*and drop subscribes

*g*in

*gu*

_{s}, so that

*V̂*

_{gusj}≡

*V̂*

_{sj},

*V̂*

^{†}

_{gusj}≡ ${\widehat{V}}_{\mathit{\text{sj}}}^{\u2020}$ and

**μ**

_{gusj}≡

**μ**

_{sj}. Under electric dipole approximation,

**Â**(

**r̂**

_{ej}(

*t*)

*t*) is replaced by

**Â**(

**r**

_{nj},

*t*), where

**r**

_{nj}is the classical position of the nucleus [Ref. 16].

The Hamiltonian can then be expressed in terms of these atomic and field operators (referred to as second quantization), which becomes *H* = *Ĥ*_{Atom} + *Ĥ*_{Field} + *Ĥ*_{AF} with

where *ω*_{a}
=*ω*_{s}
-*μ*_{g}
. The 3D atomic dipole moment operator vector **μ̂**_{j}(*t*)can be expressed in *V̂*
_{sj}
- and ${\widehat{V}}_{\mathit{\text{sj}}}^{\u2020}$ as ${\hat{\mathbf{\mu}}}_{j}\left(t\right)={\displaystyle \sum _{s}}\left[{\mathbf{\mu}}_{\mathit{sj}}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)+{\mathbf{\mu}}_{\mathit{sj}}^{*}{\hat{V}}_{\mathit{sj}}\left(t\right)\right]={\displaystyle \sum _{s}}{\hat{\mathbf{\mu}}}_{\mathit{sj}}\left(t\right)$. From the second quantized Hamiltonian, we can derive the following equations for the atomic operators using the Heisenberg equation of motion:

$$\frac{\text{d}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)}{dt}=i{\omega}_{a}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)+\frac{{\omega}_{a}{\mu}_{\mathit{sj}}^{*}}{\mathit{\u0127}}\left({\hat{n}}_{\mathit{sj}}\left(t\right)-{\hat{n}}_{\mathit{gj}}\left(t\right)\right){\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t),$$

$$\frac{\text{d}{\hat{n}}_{\mathit{sj}}\left(t\right)}{dt}=-\frac{{\omega}_{a}}{\mathit{\u0127}}\left[{\mu}_{\mathit{sj}}^{*}{\hat{V}}_{\mathit{sj}}\left(t\right)+{\mu}_{\mathit{sj}}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)\right]{\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t)=-\frac{{\omega}_{a}}{\mathit{\u0127}}{\hat{\mu}}_{\mathit{sj}}\left(t\right){\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t),$$

$$\frac{\text{d}{\hat{n}}_{\mathit{sj}}\left(t\right)}{dt}=\sum _{s}\frac{{\omega}_{a}}{\mathit{\u0127}}\left[{\mu}_{\mathit{sj}}^{*}{\hat{V}}_{\mathit{sj}}\left(t\right)+{\mu}_{\mathit{sj}}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)\right]{\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t)=\sum _{s}\frac{{\omega}_{a}}{\mathit{\u0127}}{\hat{\mu}}_{\mathit{sj}}\left(t\right){\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t),$$

where ${\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t)={\displaystyle \sum _{s\in \left\{x,y,z\right\}}}\hat{\mathbf{A}}({\mathbf{r}}_{\mathit{nj}},t)\xb7{\mathbf{e}}_{s}$. For a single electron system, we have the completeness relation ${\hat{n}}_{\mathit{gj}}\left(t\right)+{\displaystyle \sum _{s}}{\hat{n}}_{\mathit{sj}}\left(t\right)=1.$.

## 5. Multi-electron treatment

In the case of a system with multiple (2 or more) electrons, following the usual multi-electron treatment, the atomic number operators *n̂*_{gj} and *n̂*_{sj} shall be expressed in terms of the
Fermionic electron creation operators (${\widehat{c}}_{g}^{\u2020}$ , ${\widehat{c}}_{s}^{\u2020}$) and annihilation operators (*ĉ*_{g} , *ĉ*_{s}) as *n̂*_{gj}(*t*) = ${\widehat{c}}_{\mathit{\text{gj}}}^{\u2020}$(*t*)*ĉ*_{gj}(*t*) and *n̂*_{sj}(*t*) = ${\widehat{c}}_{\mathit{\text{sj}}}^{\u2020}$(*t*)*ĉ*_{sj}(*t*) [Refs. 9, 10], where these operators obey the equal-time Fermion anti-commutation relations:

$$\left\{{\hat{c}}_{{k}_{1}j}^{\u2020}\left(t\right),{\hat{c}}_{{k}_{2}j}^{\u2020}\left(t\right)\right\}=\left\{{\hat{c}}_{{k}_{1}j}\left(t\right),{\hat{c}}_{{k}_{2}j}\left(t\right)\right\}=0,$$

where *k*
_{1}, *k*
_{2}∊ {*g, s*}. The atomic transition operators are then replaced by *V̂*_{sj}(*t*) = ${\widehat{c}}_{\mathit{\text{gj}}}^{\u2020}$(*t*)*ĉ*_{sj}(*t*) and ${\widehat{V}}_{\mathit{\text{sj}}}^{\u2020}$(*t*) = ${\widehat{c}}_{\mathit{\text{sj}}}^{\u2020}$(*t*)*ĉ*_{gj}(*t*) [Ref. 17]. This procedure is valid in the free carrier limit in which the many-body Coulomb interaction effects associated with multiple electrons and other many-body effects are neglected [Ref. 9]. The many-body effects often manifested themselves as effective shift in bandgap energy, dipole relaxation behavior (e. g. Markovian vs. non-Markovian) [Ref. 11], density-dependent broadening and saturation, excitation-induced dephasing, or scattering of carrier and polarization [Ref. 9]. These effects could be incorporated phenomenologically by effectively fitting the relevant parameters and in some cases by modifying the parameters in the DTEQM-FDTD model dynamically and using the time-dependent carrier density distributions as state parameters. The Hamiltonians in Eqs. (5) and (6) then become:

Using the usual derivation of quantum Langevin equations of motion from the Hamiltonian [Ref. 18] by tracing over the thermal field reservoir states, which is equivalent to solving (using Eqs. (8)) for the atom-field operator evolutions to the first order under Markov’s (memory-free) approximation, we will obtain the spontaneous decay terms. It turns out that with use of the Fermion anticommutation relations of the electron creation and annihilation operators, the decay terms for the electron upper-level population operator *n̂*_{sj}(*t*) will be in the form of *n̂*_{sj}(*t*)(1 - *n̂*_{gj}(*t*)), which gives the Pauli Exclusion Principle [Refs. 5–8]. This is because the electron transition rate term proportional to *n̂*_{sj}(*t*)(1-*n̂*_{gj}.(*t*)) will reduce to zero when the lower level population *n̂*_{gj}(*t*) becomes 1 or fully occupied. The equations of motion for the atomic operators then become the followings:

$$\frac{\text{d}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)}{dt}=i{\omega}_{a}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)-{\gamma}_{\mathit{Vs}}{\hat{V}}_{\mathit{sj}}^{\u2020}\left(t\right)+\frac{{\omega}_{a}{\mu}_{\mathit{sj}}^{*}}{\mathit{\u0127}}\left({\hat{n}}_{\mathit{sj}}\left(t\right)-{\hat{n}}_{\mathit{gj}}\left(t\right)\right){\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t)+{\hat{\Gamma}}_{{V}_{\mathit{sj}}}^{\u2020}\left(t\right),$$

$$\frac{\text{d}{\hat{n}}_{\mathit{sj}}\left(t\right)}{dt}=-{\gamma}_{\mathit{Ns}}{\hat{n}}_{\mathit{sj}}\left(t\right)\left[1-{\hat{n}}_{\mathit{gj}}\left(t\right)\right]-\frac{{\omega}_{a}}{\mathit{\u0127}}{\hat{\mu}}_{\mathit{sj}}\left(t\right){\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t)+{\hat{\Gamma}}_{{n}_{\mathit{sj}}}\left(t\right),$$

$$\frac{\text{d}{\hat{n}}_{\mathit{gj}}\left(t\right)}{dt}={\displaystyle \sum _{s}}{\gamma}_{\mathit{Ns}}{\hat{n}}_{\mathit{sj}}\left(t\right)\left[1-{\hat{n}}_{\mathit{gj}}\left(t\right)\right]+{\displaystyle \sum _{s}}\frac{{\omega}_{a}}{\mathit{\u0127}}{\hat{\mu}}_{\mathit{sj}}\left(t\right){\hat{A}}_{s}({\mathbf{r}}_{\mathit{nj}},t)-{\hat{\Gamma}}_{{n}_{\mathit{sj}}}\left(t\right),$$

where $\widehat{\Gamma}$
’s are the Langevin noise operators with zero mean [Ref. 18], *γ*
_{Ns} is the decay rate for upper-level population *n̂*_{sj}(*t*), and *γ*
_{Vs} is the dipole dephasing rate for operators *V̂*_{sj}(*t*) and ${\widehat{V}}_{\mathit{\text{sj}}}^{\u2020}$(*t*) that constitute the atomic dipole moment operator vector **μ̂**_{j}(*t*). Note that *γ*
_{Vs} = *γ*
_{Ns}/2 + *γ*
_{ph} where *γ*
_{ph} is the additional dephasing rate due to *Ĥ*_{ph} [Ref. 13]. Taking the mean values of these equations will give us a set of mean-valued equations for the electron variables, which are used in the next section. The top two equations of Eqs. (11) can be used to derive a second-order differential equation for the dipole [Ref. 5–8]. If we specialize to only the two-level transition between one of the *n̂*_{sj}(*t*) and *n̂*_{gj}(*t*) for simplicity, we will obtain the following equation:

$$=-2\frac{{\omega}_{a}}{\mathit{\u0127}}{\mid {\mu}_{\mathit{sj}}\mid}^{2}\left[{\hat{n}}_{\mathit{sj}}\left(t\right)-{\hat{n}}_{\mathit{gj}}\left(t\right)\right]{\hat{E}}_{s}({\mathbf{r}}_{\mathit{nj}}.t),$$

where we have dropped small terms by assuming that Ω_{m} ≫ *γ*
_{Vs}, *γ*
_{Ns} The ${\widehat{A}}_{s}^{2}$ term in Eq. (12) will affect Rabi oscillation but can be neglected at low intensity. Eq. (12) and the bottom two equations of Eqs. (11) together forms the complete set of equations for the medium variables.

## 6. Summary of medium-field equations for the multi-level multi-electron FDTD model

Using the equations derived in sections 4 and 5 to describe the interband transitions, we further include the intraband transitions phenomenologically in the rate equations with the corresponding Pauli Exclusion decay terms [Refs. 5–8] and obtain the full set of mean-valued equations describing the multi-level multi-electron model below. For simplicity, we consider only one polarization direction (z) from now on and consider only *s*=*z* in all equations.

To begin with, we start with the Maxwell equations in which the electromagnetic field is coupled to the macroscopic polarization density **P**(**r**,*t*).

Let the z-component of **P**(**r**,*t*) be given by *P*_{z}
(**r**,*t*) =**e**
_{z}·**P**(**r**,*t*). The macroscopic polarization density *P*_{z}
(**r**,*t*) represents the total dipole moment per unit volume. In FDTD, the spatial region is subdivided into volume *δV* that is small compared to the wavelength of interest. The atomic dipole moment **μ**
_{zj}(*t*) at **r**
_{nj} multiplied by the total number of dipoles *N*
_{0} (in volume *δV*) divided by *δV* will give *P*_{z}
(**r**,*t*). Let the spatial region centered at **r**
_{nj} with volume *δV* be denoted as (**r**
_{nj}, *δV*), we can then express *P*_{z}
(**r**, *t*) as follows:

where *N*
_{dip-i} (**r**) is the dipole volume density given by the number of dipoles *N*
_{0} divided by volume *δV* . We will assign one macroscopic polarization density variable *P*_{z}
(**r**,*t*) for each of the transition energy *E*
_{i} between the levels *i*_*c* and *i*_*ν*. Thus *N*
_{0} will include the dipoles describing transitions centered at transition energy *E*
_{i} within energy width *δE* and volume *δV*, so *N*
_{dip-i}(**r**) will be the volume density of dipoles for level *i* within energy width *δE*. We may thus further label *P*_{z}
(**r**,*t*) with subscript *i* as *P*_{iz}
(**r**,*t*). Let us denote the volume densities of states for level *i* in the conduction band and valence band by ${N}_{\mathit{\text{Ci}}}^{0}$(**r**) and ${N}_{\mathit{\text{Vi}}}^{0}$(**r**) , respectively. We can then express the carrier volume densities *N*
_{Ci} (**r**, *t*) and *N*
_{Vi} (**r**, *t*) (within energy width *δE*) at levels *i*_*c* and *i*_*ν*, respectively, in terms of *n*
_{zj}(*t*) and *n*
_{gj}(*t*) as *N*
_{Ci} (**r**, *t*) = *n*
_{zj}(*t*) · ${N}_{\mathit{\text{Ci}}}^{0}$(**r**) and *N*
_{Vi}(**r**, *t*) = *n*
_{gj}(*t*) · ${N}_{\mathit{\text{Vi}}}^{0}$(**r**) . In general, ${N}_{\mathit{\text{Ci}}}^{0}$(**r**) and ${N}_{\mathit{\text{Vi}}}^{0}$(**r**) are not equal, but we will show below that they are equal for a simple parabolic band case. Furthermore, in a more complex semiconductor medium, not all the states are involved in the dipole transitions. Let ${N}_{\mathit{\text{Ci}}}^{A}$(**r**) and ${N}_{\mathit{\text{Vi}}}^{A}$(**r**) be the active states involving in the dipole transitions having similar dipole transition strengths, the volume density of dipoles *N*
_{dip-i}(**r**) will be given by the larger one of ${N}_{\mathit{\text{Ci}}}^{A}$(**r**) and ${N}_{\mathit{\text{Vi}}}^{A}$ (**r**) [Ref. 19].

Using these relations, the polarization equation can be obtained from Eq. (12) by multiplying both sides of the equation with *N*
_{dip-i}(**r**
_{nj}) so that the microscopic dipole moment variable *μ*
_{zj}(*t*) is replaced by the macroscopic polarization *P*_{z}
(**r**
_{nj} ,*t*), labeled as *P*_{iz}
(**r**
_{nj},*t*) (*i*=1, 2, *M* is the level number). Generalizing it to any position **r** ∊ (**r**
_{nj}, *δV*) and we obtain:

$$=\frac{2{\omega}_{\mathit{ai}}}{\mathit{\u0127}}{\mid {\mu}_{\mathit{zi}}\mid}^{2}\left[\frac{{N}_{\mathit{dip}-i}\left(\mathbf{r}\right)}{{N}_{\mathit{Vi}}^{0}\left(\mathbf{r}\right)}{N}_{\mathit{Vi}}(\mathbf{r},t)-\frac{{N}_{\mathit{dip}-i}\left(\mathbf{r}\right)}{{N}_{\mathit{Ci}}^{0}\left(\mathbf{r}\right)}{N}_{\mathit{Ci}}(\mathbf{r},t)\right]{E}_{z}(\mathbf{r},t),$$

where *μ*
_{zi}, is the z-dipole moment matrix element _{j}〈*E*
*i*_*c*|*μ̂*_{zj}|*E*
*i*_*ν*〉_{j} between levels *i*_*c* and *i*_*ν* and is given by |*μ*
_{zi}|^{2} =(3*πħε*
_{0}
*c*
^{3})/(${\omega}_{\mathit{\text{ai}}}^{3}$
*τ*
_{i}) . The rate equations for *N*
_{Ci}(**r**,*t*) and *N*
_{Vi} (**r**, *t*) in the conduction and valence bands are then given by:

$$\frac{d{N}_{\mathit{Vi}}(\mathbf{r},t)}{dt}=\Delta {N}_{i}(\mathbf{r},t)+\Delta {N}_{\left(i+1,i\right)V}(\mathbf{r},t)-\Delta {N}_{(i,i-1)V}(\mathbf{r},t)-{W}_{\mathrm{pump}}(\mathbf{r},t),$$

where the Δ*N* terms on the right hand side describe:

- 1. Intraband transition for conduction band in which Δ
*N*_{(i,i-1)c}(**r**,*t*) is the number of electrons per unit volume transferred from conduction band level*i*_*c*to*i*-1_*c*given by:$$\Delta {N}_{\left(i,i-1\right)C}(\mathbf{r},t)=\frac{{N}_{\mathit{Ci}}(\mathbf{r},t)\left[1-{N}_{C\left(i-1\right)}(\mathbf{r},t)/{N}_{C\left(i-1\right)}^{0}\left(\mathbf{r}\right)\right]}{{\tau}_{\left(i,i-1\right)C}}-\frac{{N}_{\mathit{C}\left(i-1\right)}(\mathbf{r},t)\left[1-{N}_{Ci}(\mathbf{r},t)/{N}_{Ci}^{0}\left(\mathbf{r}\right)\right]}{{\tau}_{(i-1,i)C}}$$ - 2. Intraband transition for valence band in which Δ
*N*_{(i,i-1)V}is the number of electrons per unit volume transferred from valence band level*i*_*ν*to*i*-1_*ν*and is given by:$$\Delta {N}_{\left(i,i-1\right)V}(\mathbf{r},t)=\frac{{N}_{\mathit{Vi}}(\mathbf{r},t)\left[1-{N}_{V\left(i-1\right)}(\mathbf{r},t)/{N}_{V\left(i-1\right)}^{0}\left(\mathbf{r}\right)\right]}{{\tau}_{\left(i,i-1\right)V}}-\frac{{N}_{\mathit{V}\left(i-1\right)}(\mathbf{r},t)\left[1-{N}_{\mathrm{Vi}}(\mathbf{r},t)/{N}_{\mathit{Vi}}^{0}\left(\mathbf{r}\right)\right]}{{\tau}_{(i-1,i)V}}$$ - 3. Interband driven transition (gain or absorption) and spontaneous decay in which Δ
*N*_{i}is the number of electrons per unit volume transferred from level*i*_*c*to*i*_*ν*and is given by:

In Eq. (16) we have a pumping term *W*
_{pump} on the right hand side, which allows us to describe electrical pumping. If the current density going into an active material of thickness *d* is *J*(**r**,*t*)(A/m^{2}), then *W*
_{pump} is given by W_{pump}(**r**,*t*) = *J*(**r**,*t*)/(*ed*), where *e*=1.6×10^{-19}C is the electron charge. The pumping terms are applied to only the largest *E*
_{i} level to simulate current injection in semiconductor. Optical pumping can be simulated by introducing an optical pumping beam, which will excite carriers from *N*
_{Vi}(**r**,*t*)to *N*
_{Ci}(**r**,*t*)via the first term in Eq. (19).

Note that Eqs. (16)–(19) conserve the total number of electrons in the conduction band plus the valence band so that the “total density” *N*
_{T}(**r**,*t*) ≡ *N*
_{CT}(**r**,*t*) + *N*
_{VT}(**r**,*t*) = *N*
_{T}(**r**,0) for electrons is time-independent, where ${N}_{\mathit{CT}}(\mathbf{r},t)={\displaystyle \sum _{i}}{N}_{\mathit{Ci}}(\mathbf{r},t)$ and ${N}_{\mathit{VT}}(\mathbf{r},t)={\displaystyle \sum _{i}}{N}_{\mathit{Vi}}(\mathbf{r},t).$. The field driven transition term in Eq. (19) is in the *A*·*P* form (instead of the −*E*(*dP* / *dt*) form) valid at far off resonance [Ref. 5–8]. The ${A}_{z}^{2}$ term in Eq. (15) is important for high field intensity case.

## 7. Discussion on number density

In this section, we discuss the procedure to obtain the corresponding volume density of states in each of the discrete energy levels. To obtain the correct volume density of states, we sum up all the available states within an energy width δE at the particular energy level in the original band structure. In the multi-level model, energy level *i* will encompass all the optical transitions between energy levels ${E}_{i-}^{B}$ =*E*
_{i}-(*E*
_{i}-*E*
_{i-1})/2 and ${E}_{i+}^{B}$=*E*
_{i}+(*E*
_{i+1}-*E*
_{i})/2 as shown in Fig. 3. For a parabolic band structure, the volume density of states at each energy level is calculated from the number of states per unit energy per unit volume *ρ* (*E*):

where ΔE is the energy measured from the band edges and are denoted as ΔE_{Ci} and ΔE_{Vi} for the conduction band and valence band, respectively. In a standard parabolic band structure, ΔE_{Ci} and ΔE_{Vi} are given by:

where *m*_{C}
and *m*_{V}
are the electron and hole effective masses, respectively, and *E*_{G}
is the bandgap energy. Hence, the volume density of states for energy level *i*, which is the number of states per unit volume with interband transition energy between ${E}_{i-}^{B}$ and ${E}_{i+}^{B}$ , will be

In this case, we have ${N}_{\mathit{\text{Ci}}}^{0}$(**r**) = ${N}_{\mathit{\text{Vi}}}^{0}$(**r**) = ${N}_{\mathit{\text{Ci}}}^{A}$(**r**) = ${N}_{\mathit{\text{Vi}}}^{A}$(**r**) and hence can set the volume density of dipoles *N*
_{dip-i}(*r*
**) = ${N}_{\mathit{\text{Vi}}}^{0}$( r). It is interesting to see that in this case the volume density of states for level i_c is equal to the volume density of states for level i_ν. However, because of the difference in the intraband relaxation rates for the electrons in the conduction band and valence band due to their different effective masses [Ref. 20], the number of electrons and holes in the corresponding levels i_c and i-ν may not be equal.**

**8. Modeling the Fermi-Dirac Thermalization Dynamics**

**At finite temperature, non-equilibrium electrons will transfer between intraband states through phonon-assisted processes to achieve a quasi-equilibrium distribution, which is referred to as the conduction band or valence band Fermi-Dirac thermalization. This effect is often modeled by calculating the quasi-equilibrium Fermi energy levels for both electrons and holes. In a time-domain numerical method, the Fermi energy levels will need to be updated temporally, resulting in complicated calculations. In our model, we will show that the Fermi-Dirac thermalization can be accurately achieved by assigning the right ratio between the upward and downward intraband transition rates between two neighboring levels. The upward transition mimics thermally excited carrier hopping. It turns out that the ratios between the upward and downward intraband transition rates are independent on the Fermi energy levels, which greatly simplify the simulation. To show that, let us illustrate using the intraband transition rate equation between the conduction band levels i_c and i-1_c given by Eq. (17) as follows:**

**$$\Delta {N}_{\left(i,i-1\right)C}(\mathbf{r},t)=\frac{{N}_{\mathit{Ci}}(\mathbf{r},t)\left[1-{N}_{C\left(i-1\right)}(\mathbf{r},t)/{N}_{C\left(i-1\right)}^{0}\left(\mathbf{r}\right)\right]}{{\tau}_{\left(i,i-1\right)C}}-\frac{{N}_{\mathit{C}\left(i-1\right)}(\mathbf{r},t)\left[1-{N}_{Ci}(\mathbf{r},t)/{N}_{Ci}^{0}\left(\mathbf{r}\right)\right]}{{\tau}_{(i-1,i)C}}$$**

**When the intraband transition between those two levels reaches steady state, we can set Δ N
_{(i,i-1)C} = 0. The ratio for the upward transition rate τ
_{(i-1,i)C}) and downward transition rate τ
_{(i,i-1)C)}between levels i_c and i-1_c is then given by:**

**$$\frac{{\tau}_{\left(i-1,i\right)C}}{{\tau}_{\left(i,i-1\right)C}}=\frac{{N}_{C\left(i-1\right)\_S}\left(\mathbf{r}\right)\left[1-{N}_{\mathit{Ci}\_S}\left(\mathbf{r}\right)/{N}_{\mathit{Ci}}^{0}\left(\mathbf{r}\right)\right]}{{N}_{\mathit{Ci}\_S}\left(\mathbf{r}\right)\left[1-{N}_{\mathit{C}\left(i-1\right)\_S}\left(\mathbf{r}\right)/{N}_{\mathit{C}\left(i-1\right)}^{0}\left(\mathbf{r}\right)\right]},$$**

**where N
_{Ci_S} and N
_{C(i-1,)_s} are the steady-state carrier densities in levels i_c and i-1_c that shall obey the Fermi-Dirac thermal distribution function as follows:**

**$${N}_{\mathit{Ci}\_S}=f\left(\Delta {E}_{\mathit{Ci}}\right)\xb7{N}_{\mathit{Ci}}^{0}\left(\mathbf{r}\right)=\frac{1}{{e}^{\left[\left(\frac{{E}_{\mathit{Ci}}-{E}_{\mathit{FC}}}{{k}_{B}T}\right)\right]}+1}\xb7{N}_{\mathit{Ci}}^{0}\left(\mathbf{r}\right),$$**

**where E
_{Ci} is the energy of level i_c and E
_{FC} is the Fermi energy for the conduction band, both measured with respect to conduction band edge. This gives:**

**$$\frac{{\tau}_{\left(i-1,i\right)C}}{{\tau}_{\left(i,i-1\right)C}}=\frac{\frac{1}{{e}^{(\left({E}_{C\left(i-1\right)}-{E}_{{F}_{C}}\right)/{k}_{B}T)}+1}\xb7{N}_{C\left(i-1\right)}^{0}\left(\mathbf{r}\right)\frac{{e}^{(\left({E}_{\mathit{Ci}}-{E}_{{F}_{C}}\right)/{k}_{B}T)}}{{e}^{(\left({E}_{\mathit{Ci}}-{E}_{{F}_{C}}\right)/{k}_{B}T)}+1}}{\frac{1}{{e}^{(\left({E}_{\mathit{Ci}}-{E}_{{F}_{C}}\right)/{k}_{B}T)}+1}\xb7{N}_{\mathit{Ci}}^{0}\left(\mathbf{r}\right)\frac{{e}^{(\left({E}_{C\left(i-1\right)}-{E}_{{F}_{C}}\right)/{k}_{B}T)}}{{e}^{(\left({E}_{C\left(i-1\right)}-{E}_{{F}_{C}}\right)/{k}_{B}T)}+1}}=\frac{{N}_{C\left(i-1\right)}^{0}\left(\mathbf{r}\right)}{{N}_{\mathit{Ci}}^{0}\left(\mathbf{r}\right)}{e}^{\left({E}_{\mathit{Ci}}-{E}_{C\left(i-1\right)}\right)/{k}_{B}T}.$$**

**Similar relations can be obtained for the valence band case. Thus, we see that the ratio between any adjacent pair of the up and down intraband transition times is only dependent on the energy difference between the two levels involved, the temperature, and the ratio between the volume densities of states for the two levels. It is independent on the Fermi energy levels E
_{FC} and E
_{FV}. This leads to a significant simplification of the simulation as there will be no need to update the Fermi energy level at each spatial grid point. The temperature T describes the local crystal lattice temperature of the medium and, if necessary, can be treated as a spatial-temporal parameter T(r,t) determined separately by thermal diffusion equation. This thermalization also applies to interband but large energy gap makes its contribution negligible.**

**9. FDTD implementation**

**In order to implement FDTD simulation, we divide the spatial region of interest into discrete spatial grid points [Ref. 1]. We assume the dipole to be all in the z-direction. In this case only E
_{z} will be interacting with the medium. We use ν, u, w to label the spatial grid points in the x, y, z coordinate directions, respectively. The set of discrete difference equations for the carrier’s volume-density variables are then given by:**

**$${N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}={N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n-1}+2\Delta t\left(-\Delta {N}_{i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}-\Delta {N}_{\left(i,i-1\right)C}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+\Delta {N}_{\left(i+1,i\right)C}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+{W}_{\mathrm{pump}}\right),$$**

$${N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}={N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n-1}+2\Delta t\left(\Delta {N}_{i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}-\Delta {N}_{\left(i,i-1\right)V}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+\Delta {N}_{\left(i+1,i\right)V}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+{W}_{\mathrm{pump}}\right),$$

$$\Delta {N}_{i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}=\frac{{\omega}_{\mathit{ai}}}{\mathit{\u0127}}{A}_{z}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\bullet {P}_{\mathit{iz}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+\frac{{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}(1-{N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{\mathit{iV}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w})}{{\tau}_{i}},$$

$$\Delta {N}_{\left(i,i-1\right)C}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}=\frac{{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{C}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{\left(i-1\right)\mathit{C}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{\left(i-1,i\right)C}}-\frac{{N}_{\mathit{C}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{C}i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{i\mathit{C}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{\left(i,i-1\right)C}},$$

$$\Delta {N}_{\left(i,i-1\right)V}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}=\frac{{N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{V}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{\left(i-1\right)\mathit{V}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{(i,i-1)V}}-\frac{{N}_{\mathit{V}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{V}i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{i\mathit{V}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{\left(i-1,i\right)V}},$$

$$\sum _{i=1}^{M}}{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}+{\displaystyle \sum _{i=1}^{M}}{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}={\displaystyle \sum _{i=1}^{M}}{N}_{\mathit{iC}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}={\displaystyle \sum _{i=1}^{M}}{N}_{\mathit{iV}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}.$$

$${N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}={N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n-1}+2\Delta t\left(\Delta {N}_{i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}-\Delta {N}_{\left(i,i-1\right)V}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+\Delta {N}_{\left(i+1,i\right)V}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+{W}_{\mathrm{pump}}\right),$$

$$\Delta {N}_{i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}=\frac{{\omega}_{\mathit{ai}}}{\mathit{\u0127}}{A}_{z}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\bullet {P}_{\mathit{iz}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}+\frac{{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}(1-{N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{\mathit{iV}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w})}{{\tau}_{i}},$$

$$\Delta {N}_{\left(i,i-1\right)C}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}=\frac{{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{C}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{\left(i-1\right)\mathit{C}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{\left(i-1,i\right)C}}-\frac{{N}_{\mathit{C}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{C}i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{i\mathit{C}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{\left(i,i-1\right)C}},$$

$$\Delta {N}_{\left(i,i-1\right)V}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}=\frac{{N}_{\mathit{Vi}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{V}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{\left(i-1\right)\mathit{V}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{(i,i-1)V}}-\frac{{N}_{\mathit{V}\left(i-1\right)}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}[1-{N}_{\mathit{V}i}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n}\u2044{N}_{i\mathit{V}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}]}{{\tau}_{\left(i-1,i\right)V}},$$

$$\sum _{i=1}^{M}}{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}+{\displaystyle \sum _{i=1}^{M}}{N}_{\mathit{Ci}}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}^{n+1}={\displaystyle \sum _{i=1}^{M}}{N}_{\mathit{iC}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}={\displaystyle \sum _{i=1}^{M}}{N}_{\mathit{iV}}^{0}{|}_{u-\frac{1}{2},v+\frac{1}{2},w}.$$

**For the macroscopic polarization, we have:**

**$${P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n+1}$$**

$$=\frac{4-2\Delta {t}^{2}\left({\omega}_{\mathit{ai}}^{2}+4\frac{{\omega}_{\mathit{ai}}^{2}}{{\mathit{\u0127}}^{2}}{\mid {\mu}_{i}\mid}^{2}{A}_{z}^{2}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}\right)}{2+\mathrm{\Delta t\xb7}{\gamma}_{i}}{P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}+\frac{\Delta \mathit{t\xb7}{\gamma}_{i}-2}{2+\mathrm{\Delta t\xb7}{\gamma}_{i}}{P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n-1}$$

$$-\frac{4\Delta {t}^{2}{\omega}_{\mathit{ai}}}{\mathit{\u0127}\left(2+\mathrm{\Delta t\xb7}{\gamma}_{i}\right)}{\mid {\mu}_{i}\mid}^{2}\left({N}_{\mathit{Ci}}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}-{N}_{z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}\right){E}_{z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}.$$

$$=\frac{4-2\Delta {t}^{2}\left({\omega}_{\mathit{ai}}^{2}+4\frac{{\omega}_{\mathit{ai}}^{2}}{{\mathit{\u0127}}^{2}}{\mid {\mu}_{i}\mid}^{2}{A}_{z}^{2}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}\right)}{2+\mathrm{\Delta t\xb7}{\gamma}_{i}}{P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}+\frac{\Delta \mathit{t\xb7}{\gamma}_{i}-2}{2+\mathrm{\Delta t\xb7}{\gamma}_{i}}{P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n-1}$$

$$-\frac{4\Delta {t}^{2}{\omega}_{\mathit{ai}}}{\mathit{\u0127}\left(2+\mathrm{\Delta t\xb7}{\gamma}_{i}\right)}{\mid {\mu}_{i}\mid}^{2}\left({N}_{\mathit{Ci}}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}-{N}_{z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}\right){E}_{z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}.$$

**The update of the electrical field components in z direction will be:**

**$${E}_{z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n+1}={E}_{z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n+1}-\frac{1}{\epsilon}{\displaystyle \sum _{i=1}^{M}}\left({P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n+1}-{P}_{i,z}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n}\right)$$**

$$+\frac{\Delta t}{\epsilon \Delta x}\left({H}_{y}{\mid}_{u,\nu +\frac{1}{2},w}^{n+\frac{1}{2}}-{H}_{y}{\mid}_{u-1,,\nu +\frac{1}{2},w}^{n+\frac{1}{2}}\right)-\frac{\Delta t}{\epsilon \Delta y}\left({H}_{x}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n+\frac{1}{2}}-{H}_{x}{\mid}_{u-\frac{1}{2},\nu ,w}^{n+\frac{1}{2}}\right).$$

$$+\frac{\Delta t}{\epsilon \Delta x}\left({H}_{y}{\mid}_{u,\nu +\frac{1}{2},w}^{n+\frac{1}{2}}-{H}_{y}{\mid}_{u-1,,\nu +\frac{1}{2},w}^{n+\frac{1}{2}}\right)-\frac{\Delta t}{\epsilon \Delta y}\left({H}_{x}{\mid}_{u-\frac{1}{2},\nu +\frac{1}{2},w}^{n+\frac{1}{2}}-{H}_{x}{\mid}_{u-\frac{1}{2},\nu ,w}^{n+\frac{1}{2}}\right).$$

**All other electric and magnetic field components ( E_{x}
, E_{y}
, H_{x}
, H_{y}
, H_{z}
) follow the usual updating scheme in Yee’s algorithm [Ref. 1].**

**10. FDTD simulation examples**

**In order to validate the DTEQM approach for FDTD computation of complex media, we will give several simulation examples using the model discussed above. We assume a semiconductor bulk medium with bandgap wavelength of 1550nm and model it with the simple parabolic band case. The effective masses for the conduction and valence bands are 0.046 m_{e}
and 0.36m_{e}
, respectively, with m_{e}
being the free electron mass. The energy levels E
_{i}’s are spaced by constant wavelength spacing Δλ. If we use 5 energy levels for conduction and valence band and let Δλ=50nm, then the optical transition wavelengths for the discrete levels will be 1525nm, 1475nm, 1425nm, 1375nm, and 1325nm. The interband decay rate τ
_{i} for typical direct-gap semiconductor bulk medium and quantum well of interest ranges from hundreds of picoseconds to nanoseconds [Ref. 21]. In the simulation below we use τ
_{i}=1nsec.**

**For illustration purpose, we set the downward intraband transition time τ
_{(i,i-1)C} for conduction band electrons to be about one picosecond and set the downward intraband transition time τ
_{(i-1, i)V} for valence band electrons to be about 100fs, which are within the range of values given in the literature [Ref. 20]. The upward intraband transitions for conduction and valence bands are then set to follow the ratio given in Eq. (26). The initial random distribution of carriers will relax to the quasi-steady-state Fermi-Dirac distribution within the time scale given by the intraband transition rates. The dipole dephasing time is set to be ~50fsec. As those medium time-constants are several orders of magnitude larger than the optical period, it will not affect the choice of the FDTD time step (typically 1–2 orders of magnitude smaller than the optical period).**

**10.1 Application to stead- state absorption spectrum – band filling effect**

**To show the band-filling effect, we used the semiconductor medium parameters discussed above and studied its steady-state absorption spectrum. We simulated the same medium using both 5 energy levels and 10 energy levels. The total valence band volume density of states ${N}_{\mathit{\text{VT}}}^{0}$ obtained by summing over all levels in the valence band are kept to be the same at ${N}_{\mathit{\text{VT}}}^{0}$~7×10 ^{23} m^{-3} to give reasonable bulk absorption coefficient. We electrically pumped a 2 μm long semiconductor medium along a 1 μm wide waveguide to different carrier densities. After a time period of constant pumping, the carrier density in the medium will reach steady state. A weak and short optical pulse (300fs FWHM) is then launched into the semiconductor waveguide and propagated through the medium as a probing signal to obtain medium information such as the absorption spectrum.**

**The absorption spectrum is calculated from the Fourier Transform of the output pulse comparing to the input pulse. The result is plotted in Fig. 4, where Fig. 4(a) shows the spectrum when 5 energy levels are used for the conduction and valence bands and Fig. 4(b) shows the spectrum when 10 energy levels are used with Δ λ=25nm. Comparing the result of Fig. 4 with the typical semiconductor gain/absorption spectrum [Ref. 22], we see the expected gain spectrum broadening and the shift in the peak wavelength of the gain spectrum with increasing carrier density. The gain spectrum broadening is a result of band filling and is the main reason why it is it difficult to get high inversion in semiconductor amplifiers. Although the 10-level case gives much smoother spectrum, the 5-level case can already show the essential effect covering a broad wavelength range from 1300nm to 1600nm. In many device applications, the spectral coverage does not have to be that broad and fewer levels covering ~ 100 nm can be used.**

**10.2 Application to ultrafast transient response**

**The absorption of a short laser pulse with an optical frequency above the bandgap in a semiconductor creates free electrons and holes with an initial energy distribution that is essentially determined by the optical spectrum of the laser pulse. Within a very short time scale, this non-thermal energy distribution is transformed by carrier-carrier scattering into a quasi-equilibrium Fermi-Dirac distribution [Refs. 9, 10]. This process of spectral hole burning and subsequent thermalization is a representative example of ultrafast phenomena in semiconductor [Ref. 23]. Figure 5 shows the normalized number density in the discrete energy levels as a function of time when the semiconductor waveguide is pumped by a 300fs Gaussian pulse with a peak intensity of 200MW/cm ^{2} and a center wavelength of 1425nm. We see that the electrons in the valence band relax to quasi equilibrium much faster than the electrons in the conduction band, exemplified by the steeper slope of the valence band electron relaxation curves. The transient evolution of the entire absorption spectrum can also be obtained [Ref. 12].**

**10.3 Application to simulation of multimode microdisk laser**

**Next we apply the DSM-FDTD model to simulate the lasing behavior of a microdisk semiconductor laser [Refs. 24–26] using 2-dimensional FDTD with the 5-level model on a
2GHz cpu Pentium PC. The laser gain medium is a semiconductor bulk medium given by the example above in Fig. 4(a). Let us consider a 2-μm diameter [Fig. 6(a)] microdisk for which the effective waveguiding refractive index of the disk is n=2.7 and the disk is held in air (n=1) by a pillar at the center of the disk. Here we assume the optical mode has 40% overlapping with the active medium in the vertical direction.**

**We pumped the disk with different electrical current densities and plotted the optical intensity inside the disk as a function of the injection current density (Figs. 6(b) and 6(c)). As spontaneous emission noise is not yet included in the theory, the medium is hit with an optical pulse to initiate lasing. After that, self-sustained lasing behavior can be achieved at above threshold. In our simulation, we ran for ~2000 cavity round trips to achieve the steady state and used a FDTD spatial grid resolution of ~ λ/30. The laser simulated showed a threshold current density of ~400A/cm^{2}, which corresponded to a carrier density of 2.5×10^{23} m^{-3}. This carrier density is above the transparent carrier density of 2×10^{23} m^{-3} shown in Fig. 4. The lasing spectrum inside the cavity is plotted at four different current levels above threshold as shown in Fig. 7(a). At current just above threshold, only one lasing mode is present. At pumping current density 48kA/cm^{2}, the second mode starts to appear. The two modes are the TM0_{8, 1} and TM0_{9, 1} whispering gallery modes of the microdisk cavity. Note that the lasing spectral width shown here is not the real laser linewidth as noise is not included in the current model. The “linewidth”here is actually transform-limited linewidth and is resulted from the limited time period in which we performed Fourier transform on the laser cavity field.**

**Another interesting effect is that as the pumping level increases, the carrier level increases inside the microdisk due to band filing. This leads to a change in the refractive index (Δn~0.003) of the microdisk cavity, which then leads to a change in the lasing wavelength. Fig. 7(b) shows the slight shift in the lasing wavelength as the injection current level increases. The intensity-induced lasing wavelength change as well as the multi-mode lasing effect would be difficult to simulate without the sophisticated semiconductor model introduced here.**

**11. Computational overhead and FDTD medium model extension**

**Perhaps the most important criteria in evaluating the feasibility of a FDTD semiconductor model are its computational complexity, which resulted in long computational time. The DSM-FDTD model illustrated in this paper gives a semiconductor model with minimal computational complexity that can still take into account the essential semiconductor carrier dynamics such as band filling, carrier transient effects, and carrier induced refractive index change. When using five energy-level pairs for the conduction band and valence band, the total computational time for a simulation where the entire geometry is filled with the semiconductor medium is 8 times of the same code without the medium. In a typical geometry, the active medium may occupy only a small fraction of the entire area of simulation, which can be on the order of 10%. In that case, the computational overhead is only ~ 1x. This makes the model useful for a wide range of FDTD simulations.**

**Although the current model has taken into account most of the essential semiconductor carrier dynamics, there are other effects that could potentially be included in the model depending on the particular device simulation requirements. For example, carrier diffusion could be included into the FDTD equation straight-forwardly by adding a spatial carrier transportation term in Eq. (27). Noise term can be included to give the spontaneous emission linewidth and the lasing linewidth. Carrier heating and cooling can be more accurately described by allowing the carrier temperature to evolve in time. Also, non-radiative and higher-order Auger recombination effects can be included [Ref. 9, 10]. Furthermore, quantum well and dot and their carrier capture dynamics can be modeled. However, one should only include the most essential effects sufficient for simulating the device behaviors of interest to reduce the computational burden.**

**12. Discussion and summary**

**In summary, we report a new computational model of material media capable of modeling the nanostructure and electronic dynamics of sophisticated active materials often needed in photonic device simulations. The media that can be modeled include solid-state and semiconductor type media, as well as molecular and atomic type media. This model is computationally efficient for incorporating into the FDTD electrodynamics simulation.**

**The model is based on a multi-energy-level multi-electron quantum system in which the electron dynamics is governed by the Pauli Exclusion Principle and the dynamical Fermi-Dirac Thermalization. The medium is described by a set of rate equations derived quantum mechanically. The formulation is based on the basic minimal-coupling Hamiltonian extended to incorporate multiple electrons via second quantization using Fermion creation and annihilation operators. The set of rate equations and dipole equation describing the model are obtained without the usual rotating-wave approximations and can be applied to the regimes at near or far off-resonance as well as high field intensity.**

**We refer to this model as the Dynamical-Thermal-Electron Quantum-Medium FDTD (DTEQM-FDTD) model. It is built on top of our earlier 4-level 2-electron model with Pauli Exclusion Principle [Ref. 5–8] but extended to multiple levels and multiple electrons with the important inclusion of dynamical Fermi Dirac thermalization. We show that the Fermi-Dirac thermalization can be incorporated via a temperature-dependent carrier hopping process, which mimics thermal carrier excitations. This dynamical process enables the simulation of carrier decay from non-thermal equilibrium after excitation to a quasi equilibrium carrier distribution governed by the quasi Fermi-Dirac distribution.**

**In application to semiconductor, this DTEQM-FDTD model takes into account the transient intraband and interband electron dynamics, the semiconductor band structure, and carrier thermal equilibrium process for the first time in FDTD simulation, and is referred to as the Dynamical-Semiconductor-Medium FDTD (DSM-FDTD) model. The DSM-FDTD model automatically incorporates energy-state filling effect. It also incorporates the typical nonlinear optical effects associated with carrier dynamics and thermally activated carrier scattering process under transient excitation spatial-temporally. The model also allows separate electron dynamics in the conduction and valence bands. These capabilities empower the model to treat sophisticated optoelectronic and nanophotonic devices having complex geometries with full spatial-temporal solutions at the microscopic level under electrical or optical excitation. A further extension of the FDTD model to include spatial diffusion of carriers, lattice temperature heating or cooling, and carrier dependent medium parameter shifts due to many-body effects will make it a highly powerful optoelectronic and photonic device simulator.**

**Most importantly, we show that the FDTD model is sophisticated enough to incorporate the essential multi-physical effects in complex media and yet is simple enough to achieve fast computational time. We illustrated the application of this powerful new model to FDTD computation with the simulation of the entire gain and absorption spectra of a direct-bandgap semiconductor medium, showing the carrier band filling effects with Fermi-Dirac statistics. We then illustrated the application of the FDTD model to simulating spectral hole burning of carriers under a strong optical pulse with subsequent decay to thermal equilibrium representative of ultrafast phenomena in semiconductor. To illustrate its applications to photonic devices, we simulated a microdisk laser in which a second lasing mode is excited due to gain bandwidth broadening at high medium excitation. We also show the shift in the lasing frequency with increased excitation due to carrier-induced refractive index change.**

**Acknowledgments**

**The authors would like to thank Dr. Seongsik Chang for his contribution to the earlier model and useful discussion.**

**The work is supported by NSF under Award No. ECS-0501589, by NSF MRSEC program under grant DMR–0076097, and by the NASA Institute for Nanoelectronics and Computing under Award No. NCC 2–1363.**

**References and links**

**1. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in Isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**2. **S. D. Gedney, “An anisotropic perfectly matched Layer-Absorbing Medium for the truncation of FDTD Lattice,” IEEE Trans. Antennas Propag. **44**, 1630–1639 (1996), and references therein. [CrossRef]

**3. **M. Okoniewski, M. Mrozowski, and M. A. Stuchly, “Simple treatment of multi-term dispersion in FDTD,” IEEE Microwave Guid. Wave Lett. **7**, 121–123 (1997), and references therein. [CrossRef]

**4. **A. S. Nagra and R. A. York, “FDTD analysis of wave propagation in nonlinear absorbing and gain media,” IEEE Trans. Antennas Propag. **46**, 334–340 (1998). [CrossRef]

**5. **Y. Huang, “Simulation of semiconductor material using FDTD method,” Master Thesis, Northwestern University, June 2002. https://depot.northwestern.edu/yhu234/publish/YYHMS.pdf

**6. **S. Chang, Y. Huang, G. Chang, and S. T. Ho, “THz all-optical shutter based on semiconductor transparency switching by two optical p-pulses,” OSA Annual Meeting, TuY3, Long Beach, CA, 2001.

**7. **S. T. Ho, research notes, 1998–1999.

**8. **Y. Huang, “Simulation of semiconductor structure using FDTD method”, presented to the Physics Department at Northwestern University, 15 Jan. 2002.

**9. **W. W. Chow, S. Koch, and M. Sargent III, *Semiconductor-Laser Physics*, (Springer Verlag, Berlin, 1994). [CrossRef]

**10. **J. Piprek, *Optoelectronic Devices: Advanced Simulation and Analysis*, (Springer Verlag, New York, 2005). [CrossRef]

**11. **S. Park, “Development of InGaAsP/InP single-mode lasers using microring resonators for photonic integrated circuits,” PhD Thesis, Northwestern University, Dec. 2000, and references therein.

**12. **Y. Huang and S. T. Ho, “A numerically efficient semiconductor model with Fermi-Dirac thermalization dynamics (band-filling) for FDTD simulation of optoelectronic and photonic devices,” 2005 Technical Digest of the Annual Conference on Lasers and Electro-Optics, Paper QTuD7, Baltimore, MD, May 2005.

**13. **S. T. Ho, P. Kumar, and J. H. Shapiro, “Quantum theory of nondegenerate multiwave mixing (I) - General formulation,” Phys. Rev. A **37**, 2017–2032 (1988). [CrossRef]

**14. **S. T. Ho and P. Kumar, “Quantum optics in a dielectric: Macroscopic electromagnetic-field and medium operators for a linear dispersive Lossy medium-A microscopic derivation of the operators and their commutation relations,” J. Opt. Soc. Am. B **10**, 1620–1636 (1993). [CrossRef]

**15. **S. T. Ho, P. Kumar, and J. H. Shapiro, “Vector-field quantum model of degenerate four-wave mixing,” Phys. Rev. A **34**, 293–303 (July 1986). [CrossRef]

**16. **J. J. Sakurai, *Advanced Quantum Mechanics*, (Addison Wesley,1967).

**17. ***ĉ*_{sj}(*t*) in semiconductor corresponds to the spatially localized operator${\hat{c}}_{\mathit{sj}}\left(t\right)\equiv {\hat{c}}_{s}({\mathbf{r}}_{\mathit{nj}},t)={\displaystyle \sum _{k}}{a}_{k}{\hat{c}}_{k}{e}^{\mathit{i}\mathbf{k}\xb7\mathbf{r}}.$

**18. **W. H. Louisell, *Quantum Statistical Properties of Radiation*, (Wiley-Interscience, New York, 1990).

**19. **
For example, if three upper levels can decay to a single ground level, then each upper level will be associated with a transition dipole so that the total number of dipoles involved will be three, which is equal to the number of the upper levels.

**20. **R. F. Kazarinov, C. H. Henry, and R. A. Logan, “Longitudinal mode self-stabilization in semicondcutor lasers,” J. Appl. Phys. **53**, 4631–4644 (1982). [CrossRef]

**21. **S. Marrin, B. Deveaud, F. Clerot, K. Fuliwara, and K. Mitsunaga, “Capture of photoexcited carriers in a single quantum well with different confinement structures,” IEEE J. Quantum Electron. **27**, 1669–1675 (1991). [CrossRef]

**22. **L. A. Coldren and S. W. Corzine, *Diode lasers and photonic integrated circuits*, (Wiley, John & Sons.1995).

**23. **J. L. Oudar, D. Hulin, A. Migus, A. Antonetti, and F. Alexandre, “Subpicosecond spectral hole burning due to nonthermalized photoexcited carriers in GaAs,” Phys. Rev. Lett. **55**, 2074–2077 (1985). [CrossRef]

**24. **D. Y. Chu, M. K. Chin, S. Z. Xu, T. Y. Chang, and S. T. Ho, “1.5 μm InGaAs/InAlGaAs Quantum-well microdisk lasers,” IEEE Photonics Technol. Lett. **5**, 1353–1355 (1993). [CrossRef]

**25. **W. Fang, J. Y. Xu, A. Yamilov, H. Cao, Y. Ma, S. T. Ho, and G. S. Solomon, “Large nhancement of spontaneous emission rates of InAs quantum dots in GaAs microdisks,” Opt. Lett. **27**, 948–950 (2002). [CrossRef]

**26. **J. P. Zhang, D. Y. Chu, S. L. Wu, W. G. Bi, R. C. Tiberio, C. W. Tu, and S. T. Ho, “Photonic-wire laser,” Phys. Rev. Lett. **75**, 2678–2681 (1995). [CrossRef]