Nonlinear dispersive three-dimensional finite-difference time-domain analysis for photonic-crystal lasers

The three-dimensional finite-difference time-domain method that can handle dispersive and dynamic nonlinear-gain media is proposed and realized. The effect of carrier diffusion is included through the laser rate equations. Through this three-dimensional nonlinear gain FDTD method, rich laser-dynamics behaviors, such as the lasing threshold, the relaxation oscillation, and the spatial hole burning, are directly observed from a hexapole mode. © 2005 Optical Society of America OCIS codes: (140.5960) Semiconductor lasers; (230.3990) Microstructure devices; (270.2500) Fluctuations, relaxations, and noise. References and links 1. O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, “Two-Dimensional Photonic Band-Gap Defect Mode Laser,” Science 284, 1819 (1999). 2. H. G. Park, J. K. Hwang, J. Huh, H. Y. Ryu, Y. H. Lee, and J. S. Kim, “Nondegenerate monopole-mode twodimensional photonic band gap laser,” Appl. Phys. Lett. 79, 3032 (2001). 3. H. Y. Ryu, S. H. Kim, H. G. Park, J. K. Hwang, Y. H. Lee, and J. S. Kim, “Square-lattice photonic bandgap single-cell laser operating in the lowest-order whispering gallery mode,” Appl. Phys. Lett. 80, 3883 (2002). 4. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, chap. 5, 8, 11 (Artech House, Boston, Mass, 1995). 5. H. G. Park, S. H. Kim, S. H. Kwon, Y. G. Ju, J. K. Yang, J. H. Baek, S. B. Kim, and Y. H. Lee, “Electrically Driven Single-Cell Photonic Crystal Laser,” Science 305, 1444 (2004). 6. K. Nozaki and T. Baba, “Carrier and photon analyses of photonic microlasers by two-dimensional rate equations,” J. Sel. Area. Commun., 23, 1411 (2005). 7. R. J. Luebbers and F. Hunsburger, “FDTD for n-th-order dispersive media,” IEEE Trans. Antennas Propagat., 40, 1297–1301 (1992) 8. J. Schuster and R. Luebbers, “An accurate FDTD algorithm for dispersive media using a piecewise constant recursive convolution technique,” IEEE Antennas and propagation Soc. Internat. Symp. Digest, 4, 2018 (1998). 9. L. A. Coldren and S. W. Corzine, Diode Lasers and Photonic Integrated Circuits, chap. 2 (A Wiley-Interscience Publication, 1995). 10. W. W. Chow, S. W. Koch, and M. Sargent, Semiconductor Laser Physics (Springer-Verlag, Berlin, Germany 1994), Sec. 10-4. 11. G. H. Song, S. Kim, and K. Hwang, “FDTD Simulation of Photonic-Crystal Lasers and Their Relaxation Oscillation,” J. Opt. Soc. Kor. 6, 87 (2002). 12. M. Fujita, A. Sakai, and T. Baba, “Ultrasmall and ultralow threshold GaInAsP-InP microdisk injection lasers: design, fabrication, lasing characteristics, and spontaneous emission factor,” J. Sel. Top. Quantum Electron. 5, 673 (1999). (C) 2005 OSA 28 November 2005 / Vol. 13, No. 24 / OPTICS EXPRESS 9645 #8663 $15.00 USD Received 1 September 2005; revised 3 November 2005; accepted 7 November 2005


Introduction
In recent years, varieties of prototype photonic-crystal (PhC) lasers have been developed [1,2,3].Unlike the case with the conventional semiconductor lasers, analysis approaches employing the simplified rate equations are not adequate to describe dynamics of PhC lasers because of their complicated multi-dimensional mode profiles.The finite-difference time-domain (FDTD) method [4], solving the Maxwell equations linearly, has enabled us to understand various characteristics of passive PhC resonators, such as mode profiles, quality (Q) factors, and so on.However, this method, using time-independent material parameters, is not suitable when dispersive and nonlinear media are involved.So far, there has been less effort to apply the FDTD method in studying dynamical interaction of light and matter in active PhC devices such as the electrically-driven PhC laser reported recently by the present group [5].A PhC microlaser analysis which combined the two-dimensional laser rate equation solver [6] with the FDTD mode solver in a crude fashion, which can hardly be considered self-consistent.
In this study, in order to study the dynamic interaction between the gain medium and electromagnetic fields in the PhC-laser cavity, we have developed a self-consistent scheme which integrates the active FDTD scheme taking account of the Lorentz-dispersive nonlinear medium characteristics, combined with the equation of ambipolar carrier diffusion is directly addressed under current injection.This study has revealed several interesting details in the fast transient dynamics of the PhC single-cell laser such as hole burning and suppressed relaxation oscillation in a micro-size semiconductor injection laser.

Nonlinear dispersive gain FDTD method
The gain medium is approximated by the complex electric permittivity based on the classical model of Lorentzian dipole oscillators, corresponding to the displacement vector D D D(r r r, ω) = ε(r r r, ω)E E E(r r r, ω), and is expressed in the frequency domain as where ω 0 and γ 0 is the resonance frequency and the damping constant respectively, and ε ∞ is the base-level electric permittivity.The numerator Γ(r r r) of the Lorentzian is related to the material gain G in each spatial position: where n(r r r) is the real part of the refractive index of the gain medium.Note that the material gain is dependent on the carrier density N(r r r, t).To convert these relations in the frequencydomain to those in the time domain, the Fourier recursive-convolution scheme [7] is employed for numerical implementation of D D D(r r r, t) ≡ P P P(r r r, t) where ε(r r r, t) = ∞ −∞ ε(r r r, ω) e −i ω t dω/2 π.To construct the formulas suitable for numerical implementation, we have employed the piecewise-constant recursive-convolution (PCRC) technique [8].This technique is considered the best one for computation speed and numerical accuracy in the result from the FDTD simulation of various dispersive media.A realistic medium with inhomogeneous gain spectrum can be represented by a combination of different Lorentzian oscillators spread over the region being simulated.
To obtain the material gain, which is related to the residue of the Lorentz dispersion at each position the following modified rate equation is solved at each step of FDTD computation: where q e and are the elementary charge and the Planck constant, respectively.The ambipolar carrier-diffusion effect is included in the last term in Eq. ( 4) with D being the ambipolar carrierdiffusion coefficient [10].Here, J(r r r, t) is the pump-current density at each calculation position and τ is the carrier life-time.
Note that the term for the loss of the carrier density due to nonradiative recombination is replaced by a well-accepted approximate model of with A, B, and C being the coefficients due to surface recombination, spontaneous radiative recombination, and Auger recombination, respectively [9].The term for the net generation of carriers by light emission, both stimulated and spontaneous, and absorption is contained in the right-hand side of Eq. ( 4), which is to be replaced following the recipe of according to the Poynting theorem [11].Correct handling of the last substitution is important for guaranteeing numerical stability in association with the physical energy balance in the context of FDTD calculation.The values of E E E(r r r, t) and H H H(r r r, t) are to be supplied from the FDTD calculation.When combining with the rate-equations with the FDTD method, the number of the photons created in the FDTD computation is set to be equal to the number of recombining electron-hole pairs in the rate-equation.The simulation confirmed that the afore-mentioned conservation principle is indeed well observerd.The material gain of the nonlinear quantum well (QW) is assumed to have typical logarithmic dependence on the carrier density, as in [9] G(N; r r r, t) = G 0 ln(N(r r r, t)/N tr ), (7) where N tr is the transparent carrier density.If the carrier density is larger than N tr , the material gain G(N) ∝ Γ in Eq. ( 2) becomes positive and the electric field will be amplified.
In order to represent the natural incoherent spontaneous emission, dipole sources with random phases are distributed over the entire gain medium.At each calculation step, the rate of emission from a dipole is set to be equal to that of the radiative recombination, BN 2 (r r r, t), by changing their oscillation amplitudes.In other words, the amplitudes of the dipoles are dependent on the radiative recombination rate.The employed random dipoles are distributed inhomogeneously over the frequency spectrum with the standard deviation corresponding to the half-width of the Lorentzian gain profile used above.For temporal incoherence, the phases and frequencies of the dipoles are rearranged with a tuned and randomized rate in the order of ∼ 1 ps −1 .
In order to combine the FDTD with the nonlinear rate equations properly, electric fields and magnetic fields should be updated in a proper sequence.The temporal development of the nonlinear three-dimensional FDTD method, including the gain medium, is schematically explained in Fig. 1.Remembering that the gain directly amplifies electric field amplitude, the latter should be calculated from the convolution integral, as shown before, which involves the gain dispersion under the FDTD Fourier recursive-convolution technique.The field from spontaneous dipole emission is superposed into the updated electric field.Finally, the material gain in each spatial position is calculated by the laser rate equation.Thus, the electric field and the material gain are fed back to each other during each iteration step.The time increment ∆t of ∆/2 c, smaller than the critical value of ∆/ √ 3 c, was used in our computation to satisfy the condition of numerical stability [4], where c and ∆ are the speed of light and the space increment used for the finite differences in each field component calculation respectively.The scheme of uniaxial perfectly matched layers (UPML) is used for the realization of an appropriate absorbing boundary condition.

Simulation results and discussion
The InGaAsP-InP multiple quantum wells (MQWs) are employed as the gain material.Parameters at room temperature used in the simulation are listed in Table 1 [12].The intrinsic carrier density 1.5 × 10 12 cm −3 is used as the initial carrier density, and thus the QW medium is absorptive in the beginning (t = 0).
A simple PhC single-cell laser is selected as a test structure as shown in Fig. 2. The radii of the outer air holes and the nearest neighbor holes are 0.35 a and 0.20 a, respectively, with a being the lattice constant.The thickness of the PhC slab is 0.5 a.The nonlinear gain medium is buried in the PhC slab as shown in Fig. 2(b).The size of the calculation domain in the FDTD simulation is 10 a × 8 a × 6.5 a.
We choose a hexapole mode [Fig.3(d)], as a test vehicle, artificially.To selectively excite the hexapole mode, the resonance frequency of the hexapole mode is positioned at the center of the gain spectrum, as shown in Fig. 4, so that it prevails over other modes.The Q-factor of this hexapole mode is relatively low at 550, because just two rows of surrounding air holes  cannot contain photons efficiently.Choosing the mode with such a relatively low Q helps the simulation in terms of reducing the excessive computation time in reaching the steady state.The electrical current is set at 10 mA, well above the threshold.Such a strong pump current reduces the turn-on delay of the starting laser, which also saves the computation time substantially.The current is uniformly supplied over the circular region of radius 2.5 a, as shown in Fig. 2(c).
To study the temporal development of the laser mode, the photon number participating in the hexapole laser mode and the modal gain are computed with the result plot given in Fig. 5.The modal gain is evaluated as the material gain weighted by the electric-field intensity profile of the resonance mode in question, as where G mod and G mat are the modal and the material gain, respectively.The integration was taken over the whole FDTD calculation domain.Before lasing, the injected carriers recombine inducing radiative spontaneous emission.However, the gain medium is assumed to be initially absorptive and photons from the spon-   laser mode, the spatial hole-burning effect is clearly observable as shown in Fig. 6.Due to the hole burning, the spatial profile of the carrier density would approximately give an inverted plot from the electric field intensity, which is the inset in Fig. 6(d).The laser mode preferentially consumes the carriers shaping up the region of the strong electric field with a characteristic carrier-density profile.
It is also interesting to observe the relaxation oscillation after the onset of the laser mode.This oscillation gradually decays out.It is expected that both the photon number and the modal gain approach their asymptotic steady-state values in a few ns, since both of the carrier life time and the diffusion time are in the order of ns [9].In fact, the time span of this length is very long in a typical FDTD simulation.After a couple of relaxation oscillations, the modal gain and the photon number approach their steady-state value of 35 cm −1 and 13 000, respectively, in Fig. 5.These values for the modal gain and the photon number agree reasonably well with those obtained by the simple rate equation analyses.Note that the carrier spatial distribution becomes gradually blurred owing to carrier diffusion as shown in Fig. 6.

Summary
In summary, we have developed the three-dimensional dispersive nonlinear-gain FDTD method using Lorentz-gain/dispersion profile which can be applied to the numerical simulation of the semiconductor device with active gain media.Using this method, we have studied lasing dynamics of a PhC single-cell hexapole mode laser.The study includes the creation of the laser mode from spontaneous emission, the spatial hole-burning effect, and the relaxation oscillation.With further refinement, we hope that the method could also be used to study some semiclassical aspects of cavity quantum electrodynamics.

Fig. 1 .
Fig. 1.The scheme of the gain FDTD method.

Fig. 2 .
Fig. 2. The PhC cavity structure in the simulation.The gray area in (a) and (b) is the nonlinear gain region in the PhC slab.The black area in (c) and (d) is the spatially uniform current pumping area with the radius of 2.5 a.

Fig. 4 .
Fig. 4. (a) The Lorentz dispersion of the imaginary part of ε and the locations of the resonant modes of PhC single cell cavity.(b) The H z profiles of the resonant modes.

Fig. 5 .Fig. 6 .
Fig. 5. Dynamics of the hexapole photonic-crystal laser mode.(a) The solid and dashed curves represent the changing number of the participating photons and the temporal variation of the modal gain, respectively.Lasing action begins at A. (b) Temporal behavior of the H z field.

Table 1 .
Parameters of gain medium used in simulation