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.

 

Fig. 1. four-level two-electron model

Download Full Size | PDF

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{ħ}{\omega }_a{\hat{N}}_u;H_{\mathit{Field}}={\displaystyle \sum _{k,\sigma }}\hslash {\omega }_k({\hat{a}}_{k\sigma }^{+}{\hat{a}}_{k\sigma }+1⁄2);H_{\mathit{AF}}=-\hat{\mathit{\mu }}·\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|er ̂|g〉. The first-order differential equations of the dipole operator with the empirical dephasing term γV̂ are given by:

We define the direction of the dipole along the z-axis such that µ=µe _𝑧. Therefore, the dipole operators are

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 ω ₁₂ 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 ω₃₀ is the pumping frequency.

We neglect the Rabi oscillation term 2${\omega }_{30}^2$(µb·A)²/ħ ² in Eq. (8), which is important only when the external electric field is very high. This yields the governing equations

where P_a and P_b have the resonant frequencies ω_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πε ₀ c ³/(${\omega }_{21}^2$ τ ₂₁). 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^-13s^-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

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 µ=µ ₀ (1-N). Here, µ ₀ 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:

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·(dP/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 τ₃₂ and τ₃₀ respectively. Electrons in Level 2 spontaneously decay to Levels 1 with decay time constants τ₂₁.

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 Pb. Here, the use of explicit second-order finite-differences centered at time-step n requires only knowledge of E at n. Next, we apply Eq. (12) to update E to time-step n+1. Next, we apply Eq. (11) to update N ₃ to time-step n+1. This involves the Pauli exclusion term (1-N ₂), and N₂ is taken from time-step n. Next, we apply in sequence Eq. (11) to update N₂ and N₁ . This allows N₀ to be calculated by using the conservation of electron populations. Finally we update H to time-step 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 τ₃₂ =τ₁₀ =100 fs and τ₂₁ =τ₃₀ =300 ps (which takes into account the dephasing time). The initial population density is N₁ =N₀ =5×10²³/m³. 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.

 

Fig. 2. Electron population density probability. Left shows the inversion between levels 1 and 2.

Download Full Size | PDF

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].

 

Fig. 3. Intensity output of pumping and lasing signals.

Download Full Size | PDF

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.

 

Fig. 4. Output intensity vs. pump intensity.

Download Full Size | PDF

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.