## Abstract

We report a new finite-difference time-domain (FDTD) computational model of the lasing dynamics of a four-level two-electron atomic system. Transitions between the energy levels are governed by coupled rate equations and the Pauli Exclusion Principle. This approach is an advance relative to earlier FDTD models that did not include the pumping dynamics, or the Pauli Exclusion Principle. Further, the method proposed in this paper is more versatile than the conventional modal expansion of the electromagnetic field for complex inhomogeneous laser geometries constructed in photonic crystals or light-localizing random media. For such complex geometries, the lasing modes are either difficult or impossible to calculate. The present work aims at the self-consistent treatment of the dynamics of the 4-level atomic system and the instantaneous ambient optical electromagnetic field. This permits in principle a much more robust treatment of the overall lasing dynamics of four-level gain systems integrated into virtually arbitrary electromagnetic field confinement geometries.

©2004 Optical Society of America

## Corrections

Allen Taflove, "Finite-difference time-domain model of lasing action in a four-level two-electron atomic system: erratum," Opt. Express**14**, 1702-1702 (2006)

https://opg.optica.org/oe/abstract.cfm?uri=oe-14-4-1702

## 1. Introduction

The finite-difference time-domain (FDTD) computational solution of Maxwell’s equations [1] is attracting increasing interest in full-wave simulations of the dynamics of both passive and active photonic systems. The most basic FDTD technique used in simulating lasing media involves incorporating a classical dispersive Lorentzian gain via the auxiliary differential equation (ADE) method [1, 2]. A second, more advanced ADE technique for this purpose incorporates rate equations that govern the time-domain dynamics of the atomic populations in the lasing medium [3]. Another advanced approach uses the density-matrix method along with FDTD to solve the semi classical Maxwell-Bloch equation via an iterative predictor-corrector method [4]. With increasing interests of researches in quantum optics, such as laser cooling for atoms, quantum computing, enhanced spontaneous emissions in microcavities, there is a need to treat the system either semi classically or fully quantum mechanically to explain these phenomena [6,7].

In this paper we report a new FDTD model of the lasing dynamics of a four-level two electron atomic system. Our model incorporates a simplified quantized electron energies that provide four energy levels for each of two interacting electrons. Transitions between the energy levels are governed by coupled rate equations and the Pauli Exclusion Principle (PEP) [5]. Our new approach is an advance relative to the methods described in [3, 4] which do not include the pumping dynamics for a four-level system, or the PEP. With this semi classical approach, we treat the atom quantum mechanically, and the electromagnetic wave classically. This model could be potentially further extended to a full quantum mechanical treatment of the field-matter interaction.

## 2. Method

Figure 1 illustrates electron transitions in our simplified four-level two-electron model. These transitions are treated as two coupled dipole oscillators. Levels 1 and 2 correspond to dipole *P*_{a} and Levels 0 and 3 correspond to dipole **P**_{b}.

The two-level Bloch equation is used for each oscillator. This is now derived. The atom-photon Hamiltonian for a two-level system can be expressed as:

where

${H}_{\mathit{Atom}}=\mathit{\u0127}{\omega}_{a}{\hat{N}}_{u}\phantom{\rule{.4em}{0ex}};\phantom{\rule{.4em}{0ex}}{H}_{\mathit{Field}}={\displaystyle \sum _{k,\sigma}}\hslash {\omega}_{k}({\hat{a}}_{k\sigma}^{+}{\hat{a}}_{k\sigma}+1\u20442)\phantom{\rule{.4em}{0ex}};\phantom{\rule{.4em}{0ex}}{H}_{\mathit{AF}}=-\hat{\mathit{\mu}}\xb7\hat{\mathit{E}}$

*N̂*
_{u}
=|*u*><*u*|=number operator for the upper-level |*u*>

*N̂*
_{g}
=|*g*><*g*|=number operator for the upper-level |*g*>

*ħω*_{a}
=energy difference between |*u*> and ground level |*g*>

*a*â+_{kσ}=photon-creation operator; *k*=mode number; σ=polarization direction

In Eq. (1), the electric field is given by

and the dipole operator is

where *e* is the electron charge and *V̂*
_{ll}
^{′}=|*l*><*l*
^{′}′| is the atomic transition operator. For a two-level system, the dipole operator can be written as * µ
̂*=

**µ***V̂*†+

***

**µ***V̂*where

*=〈*

**µ***u*|

*e*|

**r**̂*g*〉. The first-order differential equations of the dipole operator with the empirical dephasing term

*γV̂*are given by:

$$\frac{d{\hat{V}}^{+}}{\mathit{dt}}=\frac{i}{\hslash}[\hat{H},\hat{V}]=-i{\omega}_{\mathit{a}}{\hat{V}}^{+}-\gamma {\hat{V}}^{+}+\frac{{\omega}_{\mathit{a}}}{\hslash}\left[{\hat{N}}_{2}-{\hat{N}}_{1}\right]{\mathit{\mu}}^{*}\xb7\hat{\mathit{A}}$$

We define the direction of the dipole along the z-axis such that * µ*=

*µe*

_{𝑧}. Therefore, the dipole operators are

$$=i{\omega}_{a}\left(\mu {\hat{V}}^{\u2020}-{\mu}^{*}\hat{V}\right)-\gamma \hat{\mu}+\frac{2{\omega}_{a}}{\hslash}\left[{\hat{N}}_{2}-{\hat{N}}_{1}\right]{\mid \mu \mid}^{2}\left({\mathit{e}}_{z}\xb7\hat{\mathit{A}}\right)$$

Next we derive the second-order differential equation for the dipole operators in a two-level system. Using *d*/*dt*(*µV̂*†-*µ***V̂*)=*iω*_{a}*µ̂*, Eq. (5) becomes

The polarization density *P*_{a}
(*t*) between |1> and |2> is then

where *ω*
_{12} is the resonance angular frequency; *µ*_{a}
is the dipole matrix element between levels |1> and |2>; *γ*
_{a} is the dephasing rate for *P*
_{a}(*t*); *A*Â is the vector potential; and *N*_{i}
is the atomic population density in level *i*. A similar equation holds for the polarization density *P*_{b}
(*t*) between |0> and |3>:

where ω_{30} is the pumping frequency.

We neglect the Rabi oscillation term 2${\omega}_{30}^{2}$(* µb*·

*A*)

^{2}/

*ħ*

^{2}in Eq. (8), which is important only when the external electric field is very high. This yields the governing equations

$$\frac{{d}^{2}{\mathit{P}}_{b}}{{\mathit{dt}}^{2}}+{\gamma}_{b}\frac{d{\mathit{P}}_{b}}{\mathit{dt}}+{\omega}_{b}^{2}{\mathit{P}}_{b}={\zeta}_{b}({N}_{3}-{N}_{0})\mathit{E}$$

where * P_{a}* and

*have the resonant frequencies*

**P**_{b}*ω*

_{a}and

*ω*

_{b};

*ħω*

_{a}is the energy difference between Levels 1 and 2;

*ħω*

_{b}is the energy difference between Levels 0 and 3; and

*ζ*

_{a}is 6

*πε*

_{0}

*c*

^{3}/(${\omega}_{21}^{2}$

*τ*

_{21}). In Eq. (9), the driving terms are proportional to the population differences and the damping coefficients

*γ*

_{a}and

*γ*

_{b}simulate the non-radiation loss. For the example discussed later,

*γ*

_{a}=

*γ*

_{b}=10

^{-13}s

^{-1}. Further note that the electric field E is an instantaneous value that is composed of contributions from both the pumping and emission signals.

The population operator for the upper level is expressed as

$$=-\gamma (1-{N}_{g}){\hat{N}}_{u}+\frac{1}{\hslash {\omega}_{a}}\frac{\partial \hat{\mathit{\mu}}}{\partial t}\xb7\hat{\mathit{E}}$$

where *N*_{u}
and *N*_{g}
are defined in Eq. (1). Here, we include the spontaneous decay to the lower level by the term -*γ*(1-*N*_{g}
)*N̂*
_{u}
.

Due to the PEP, the presence of electrons within one energy level reduces the efficiency of the pumping or relaxation from other levels since each quantum state can be occupied by only one electron. Similar to semiconductor band structure, for interband interactions (3→0 or 2→1), this results in pumping blocking which takes the form * µ*=

**µ**_{0}(1-

*N*). Here,

**µ**_{0}is the quantum efficiency when there are no electrons in the active region, and

*N*is the electron population density probability. As a consequence, the quantum efficiency drops by a factor of (1-

*N*). Similar to interband relaxations, intraband transitions (3→2 or 1→0) are reduced by a factor (1-

*N*) due to the Fermi distribution of the electron population within the band.

By coupling two dipole oscillators, P_{b} (between level 0 and 3) and P_{a} (between level 1 and 2), The preceding considerations lead to the following rate equations for electron populations within four levels:

$$\frac{{\mathit{dN}}_{2}}{\mathit{dt}}=\frac{{N}_{3}\left(1-{N}_{2}\right)}{{\tau}_{32}}-\frac{{N}_{2}\left(1-{N}_{1}\right)}{{\tau}_{21}}+\frac{1}{\hslash {\omega}_{a}}\mathit{E}\xb7\frac{d{\mathit{P}}_{a}}{\mathit{dt}}\text{}$$

$$\frac{{\mathit{dN}}_{1}}{\mathit{dt}}=\frac{{N}_{2}\left(1-{N}_{1}\right)}{{\tau}_{21}}-\frac{{N}_{1}\left(1-{N}_{0}\right)}{{\tau}_{10}}-\frac{1}{\hslash {\omega}_{a}}\mathit{E}\xb7\frac{d{\mathit{P}}_{a}}{\mathit{dt}}$$

$$\frac{{\mathit{dN}}_{0}}{\mathit{dt}}=\frac{{N}_{3}\left(1-{N}_{0}\right)}{{\tau}_{30}}+\frac{{N}_{1}\left(1-{N}_{0}\right)}{{\tau}_{10}}-\frac{1}{\hslash {\omega}_{b}}\mathit{E}\xb7\frac{d{\mathit{P}}_{b}}{\mathit{dt}}\text{}$$

Here *N*_{i}
is the electron population density probability in Level *i* and *τ*_{ij}
is the decay time constant between levels *i* and *j*. The electron populations vary with pumping * E*·(

*d*/

**P***dt*) and spontaneous emission decay (

*N*

_{i}-

*N*

_{j})/

*τ*

_{ij}. Electrons in Level 3 spontaneously decay to Levels 2 and 0 with decay time constants τ

_{32}and τ

_{30}respectively. Electrons in Level 2 spontaneously decay to Levels 1 with decay time constants τ

_{21}.

Our model couples Eq. (9) and Eq. (11) with the Maxwell-Ampere law:

This permits investigating both transient processes as well as asymptotic behavior at longer time scales.

We now briefly summarize the computational algorithm. At time step *n*+1, we first implement Eq. (9) to update * Pa* and

*. Here, the use of explicit second-order finite-differences centered at time-step*

**P**b*n*requires only knowledge of

*at*

**E***n*. Next, we apply Eq. (12) to update

*to time-step*

**E***n*+1. Next, we apply Eq. (11) to update

*N*

_{3}to time-step

*n*+1. This involves the Pauli exclusion term (1-

*N*

_{2}), and

*N*

_{2}is taken from time-step

*n*. Next, we apply in sequence Eq. (11) to update

*N*

_{2}and

*N*

_{1}. This allows

*N*

_{0}to be calculated by using the conservation of electron populations. Finally we update

*to time-step*

**H***n*+3/2 by applying the Maxwell-Faraday law. This algorithm avoids the need to use the predictor-corrector algorithm of [4].

## 3. Results

To illustrate our model, we investigate a one-dimensional, optically pumped, single-defect, distributed Bragg reflector (DBR) laser cavity. Each DBR has three layers of refractive indices alternating between 1.0 and 2.0 with thickness 375 nm and 187.5 nm, respectively. The gap between the DBRs is 750 nm, corresponding to a cold-cavity defect mode of 1.5 µm. This mode has a cold-cavity Q of 100, which is sufficiently small to achieve lasing with a relatively short FDTD running time. In each DBR, the passband is centered at 0.75 µm, which allows the pumping light at this wavelength to escape the cavity. Further, the stopband is centered at 1.5 µm, which confines the lasing mode at this wavelength within the cavity. To computationally extract the lasing signal, we record the calculated output time-waveform on one side of the DBR cavity, implement a fast Fourier transform, apply a flat-topped Gaussian filter function, and fast-Fourier-transform back to time domain.

In the present example, we choose *τ*_{32}
=*τ*_{10}
=100 fs and *τ*_{21}
=*τ*_{30}
=300 ps (which takes into account the dephasing time). The initial population density is *N*_{1}
=*N*_{0}
=5×10^{23}/m^{3}. With these parameters we estimate the required density probability of population inversion to be approximately 8.4×10^{-4}.

Figure 2 shows the calculated time evolution of the electron populations in all four levels. After the onset of pumping, the population in Level 1 starts to decrease, and the population in Level 2 starts to increase. This is because pumping moves electrons from Level 1 to Level 0 and then to Level 3 and finally to Level 2.

Because of the relatively long decay time from Level 2 to Level 1, the population of Level 2 increases until inversion relative to Level 1 is reached. As shown in Fig. 3, the laser output signal jumps upward at this point, followed by fluctuations that have been predicted in the literature [6].

To obtain the lasing threshold, we plot the output intensity vs. the pump intensity as shown in Fig. 4 for four cases: (a) two-electron model with PEP; (b) two-electron model without PEP; (c) one-electron model with PEP; and (d) one-electron model without PEP. For both cases with the PEP, the lasing threshold is evidenced by a physically plausible sudden jump of the output intensity. On the other hand, omission of the PEP for the one-electron case results in a nonphysical zero-threshold laser operation. Further, omission of the PEP for the two-electron case results in a nonphysical situation where both electrons occupy the ground state and there is no light output for the pump intensity considered. These calculations illustrate the importance of incorporating the Pauli Exclusion Principle.

## 4. Summary and discussion

In summary, we have reported a new FDTD computational model of the lasing dynamics of a four-level two-electron atomic system. Transitions between the energy levels are governed by coupled rate equations and the Pauli Exclusion Principle. Our model combines classical electrodynamics and quantum mechanics for improved understanding of the interaction of electrons and electromagnetic waves.

Specifically, we believe that the new FDTD model is more versatile than previous models employing modal expansions of the electromagnetic field, especially for complex inhomogeneous laser geometries such as those constructed in photonic crystals or light-localizing random media. For such complex geometries, the lasing modes are either difficult or impossible to calculate. Our new model does not require knowledge or calculation of the lasing modes.

In addition, during laser turn-on, turn-off, or pulsing, the generation of broadband transient electromagnetic energy cannot be well described by modal theory, which intrinsically represents sinusoidal steady-state phenomena. The instantaneous nonlinear interactions of such transients with the lasing medium could be important, but are difficult to model using the modal-expansion method. These interactions are naturally treated in our model.

Overall, the present work aims at the self-consistent treatment of the dynamics of the 4-level atomic system and the instantaneous ambient optical electromagnetic field. This permits in principle a more robust treatment of the overall lasing dynamics of four-level gain systems integrated into virtually arbitrary electromagnetic field confinement geometries.

## Acknowledgments

The authors thank H. Cao for valuable discussions, and L. Terry and B. Anderson for assistance in running the FDTD models.

## References and links

**1. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Boston: Artech House, 2000).

**2. **S. C. Hagness, R. M. Joseph, and A. Taflove, “Subpicosecond electrodynamics of distributed Bragg reflector microlasers: Results from finite difference time domain simulation,” Radio Sci. **31**, 931, (1996) [CrossRef]

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

**4. **R. W. Ziolkowski, J. M. Arnold, and D. M. Gogny, “Ultrafast pulse interactions with two-level atoms,” Phys. Rev. A **52**, 3082, (1995). [CrossRef] [PubMed]

**5. **A. E. Siegman, *Lasers* (University Science Books, 1986).

**6. **J. Singh, *Semiconductor Optoelectronics Physics and Technology* (New York: McGraw-Hill, 1995)

**7. **M.O Scully and M.S. Zubairy, *Quantum Optics* (Cambridge, 1997)