Quantitative Verification of Ab Initio Self-consistent Laser Theory

We generalize and test the recent"ab initio"self-consistent (AISC) time-independent semiclassical laser theory. This self-consistent formalism generates all the stationary lasing properties in the multimode regime (frequencies, thresholds, internal and external fields, output power and emission pattern) from simple inputs: the dielectric function of the passive cavity, the atomic transition frequency, and the transverse relaxation time of the lasing transition. We find that the theory gives excellent quantitative agreement with full time-dependent simulations of the Maxwell-Bloch equations after it has been generalized to drop the slowly-varying envelope approximation. The theory is infinite order in the non-linear hole-burning interaction; the widely used third order approximation is shown to fail badly.

state. However time-independent methods to find these stationary properties in the multi-mode regime for an open laser cavity did not exist until the recent work of Tureci et al. [2,3,4] presented an "ab initio" self-consistent (AISC) formalism which generates all of the lasing properties including the output power and emission pattern from a few simple inputs. The laser cavity can be of arbitrary complexity and openness, including, e.g., chaotic dielectric disk lasers [5], photonic lattice defect mode lasers [6] and random lasers [4,7]. Here we generalize this infinite order non-linear theory by extending it beyond the slowly-varying envelope approximation. With this improvement it gives remarkably good agreement with time-dependent simulations of the Maxwell-Bloch (MB) equations, while the standard third order approximation to the non-linear hole-burning interaction fails badly.
The AISC theory builds on the original ideas of Haken and coworkers [8,9] that the inversion of the lasing medium will be approximately time-independent when γ ≪ γ ⊥ , ∆ (where γ ⊥ is the transverse (polarization) relaxation rate, ∆ is the typical mode spacing, and γ is the longitudinal (inversion) relaxation rate). The only significant approximation in the theory presented below is this approximation of stationary inversion (SIA) (well-satisfied in many lasers of interest). In addition to the excellent agreement we find between the AISC theory and the MB simulations when the ratios γ /γ ⊥ , γ /∆ are very small, we develop below a perturbative treatment of the beating terms which are neglected in SIA to extend the theory to larger values of these ratios.
The key improvements contained in the AISC theory are: 1) Treatment of the openness of the cavity exactly. 2) Inclusion of the space-dependent non-linear modal interactions (spatial holeburning) to all orders, in contrast to standard third-order treatments [1]. We show below that the third order treatment fails quantitatively and qualitatively even for the simple laser resonator we study here.
To perform a well-controlled comparison of MB and AISC results we chose to study the simple one-dimensional microcavity edge emitting laser [2,3] consisting of a perfect mirror at the origin and a dielectric region of uniform index n 0 and length L terminating abruptly on air (see inset, Fig. 1). The MB equations are simulated in time and space using a FDTD approach for the Maxwell equations, while the Bloch equations are discretized using a Crank-Nicholson scheme. To avoid solving a nonlinear system of equations at each spatial location and time step, we adopt the method proposed by Bidégaray [10], in which the polarization and inversion are spatially aligned with the electric field, but are computed at staggered times, along with the magnetic field. Modal intensities are computed by a Fourier transform of the electric field at the boundary after the simulation has reached the steady state.
The AISC theory consists of a set of coupled non-linear time-independent integral equations of size equal to the number of lasing modes at a given pump. The AISC theory presented in Refs. [2,3] is a solution to the MB equations [1] after two standard simplifications, the rotatingwave approximation (RWA) and slowly-varying envelope approximation (SVEA). The SVEA involves factoring out the rapid time-dependence of the electric field and the atomic polarization field at the atomic frequency, k a , (here and below we set c = 1 and use frequency and wavevector interchangeably), and neglecting the second time derivatives in the Maxwell wave equation of the remaining envelope function of the fields. The resulting non-linear system contains only first time-derivatives of the fields and the inversion and is sometimes referred to as the Schrödinger-Bloch (SB) equations [11]. The SVEA works well when the cavity frequencies are negligibly shifted from the atomic frequency. For microcavities these shifts need not be negligible; our simulated cavities have n 0 k a L = 60 corresponding to roughly ten wavelengths of radiation inside the cavity, approaching the microcavity limit. For the case studied here, we find noticeable discrepancies between MB solutions and the AISC theory of Refs. [2][3][4] which incorporates the SVEA (see Fig. 1). This motivated us to generalize the AISC theory to drop the SVEA (the RWA is found to be well-satisfied in all cases). (1). Excellent agreement is found with no fitting parameter. Colored lines represent individual modal output intensities; the black lines the total output intensity. Dashed lines are results of AISC calculations when the slowly-varying envelope approximation is made as in Ref. [3] showing significant quantitative discrepancies. For example, in the n = 3 case the differences of the third/fourth thresholds between the MB and AISC approaches are 46% / 63%, respectively, but are reduced to 3% and 15% once the SVEA is removed. The spectra at D 0 = 10 and the gain curve are shown as insets in (a) and (b) with the solid lines representing the predictions of the AISC approach (Eq. (1)) and with the diamonds illustrating the height and frequency of each lasing peak. The schematic in (a) shows a uniform dielectric cavity with a perfect mirror on the left and a dielectric-air interface on the right.
The generalization was as follows. Again stationary periodic solutions are assumed for the electric field E( and for the polarization fields, which oscillate at unknown lasing frequencies, k µ . The spatial variation of the field amplitude Ψ µ (x x x) is also unknown, and not assumed to be a cavity resonance, but is determined self-consistently. For a high finesse cavity and near the first threshold it was shown [2] that Ψ µ (x x x) is well-approximated by a single cavity resonance, but above threshold and for lower finesse this is not at all the case [3,4]. The treatment of the matter equations does not involve the SVEA and is exactly the same as in Ref. [2] , where the key assumption is stationary inversion, which allows the non-linear polarization in the Maxwell equation to be replaced by a non-linear function of the electric field itself. The new element is that we keep the second time-derivative of the polarization in the Maxwell equation and evaluate it by differentiating the polarization equation. The resulting improved AI/MB equations for the mode functions Ψ µ and the frequencies k µ are: Here is the dielectric function of the cavity (for the microcavity edge emitting laser n(x x x) = n 0 and the pump is assumed uniform (d 0 (x x x) = 0)). Electric field and pump strength are measured in units where g is the dipole matrix element of the gain medium. With a slight change in notation this equation differs from that derived in Ref. [2] only by the additional factor k 2 µ /k 2 a multiplying the integral. This is consistent with the expectation that the SVEA is good when the lasing frequency is very close to the atomic frequency; incorporating this change into the iterative algorithm for solving the system is trivial, but crucial for quantitative agreement with the current MB simulations. However we do not find qualitative changes from dropping the SVEA, either for this simple laser or for the more complicated two-dimensional random laser studied elsewhere [4].
In Figs. 1a, 1b we show results for n 0 = 1.5, 3, finding remarkable agreement between MB and AI approaches with no fitting parameters for the case of three mode and four mode lasing respectively. Both thresholds, modal intensities and frequencies are found correctly by the AISC approach. Note that all of these quantities are direct outputs of the AISC theory, whereas they must be found by numerical fourier analysis of the MB outputs, which can introduce some numerical error. If the earlier AI approach is used with the SVEA then significant discrepancies are found, for example for n 0 = 3 the threshold of the third mode is found to be higher than from MB while the fourth threshold is too low. Note however that the theory with the SVEA does get the right number of modes and the correct linear behavior for large pumps. We believe that this original AISC approach does solve very accurately the Schrödinger-Bloch equations and that the same discrepancies would arise between MB and SB simulations, although we haven't confirmed this.
Almost all studies of the MB equations in the multimode regime have used the approximation of treating the non-linear modal interactions to third-order (near threshold approximation) and in fact this approximation is used quite generally and uncritically throughout laser theory. From examination of the form of Eq. (1) it is clear that it treats the non-linear interactions to all orders, while the third order treatment would arise from expanding the denominator to the leading order in |Ψ ν | 2 . This third-order version of the AISC theory then becomes similar to standard treatments of Haken [8,9], with the improvement of correctly treating the openness of the cavity and the self-consistency of the lasing modes in space [2,3]. An early version of this improved third order theory was found to have major deficiencies: it predicted too many lasing modes and the intensities did not scale linearly at large pump, but exhibited a spurious saturation [2,12]. In Fig. 2 we present comparisons of the third order approximation to Eq. (1), which is improved over Ref. [2] because it drops the SVEA and the "single-pole approximation" used there. We find that this improved third order theory still does a very poor job of reproducing the multimode MB results: it still predicts too many modes (in this case seven, when there should only be four at D 0 = 10 in Fig. 2), and shows the same spurious saturation as found earlier [2] because the third-order approximation cannot give the correct linear behavior for large pumps. The infinite order treatment of Eq. (1) is both qualitatively and quantitatively essential.  Fig. 1b. The dashed lines are the results of the third order approximation to Eq. (1). The frequently used third order approximation is seen to fail badly at a pump level roughly twice the first threshold value, exhibiting a spurious saturation not present in the actual MB solutions or the AI theory. In addition, the third order approximation predicts too many lasing modes at larger pump strength. For example, it predicts seven lasing modes at D 0 = 10, while both the MB and AISC show only four. Right inset just shows the same data on a larger vertical scale.
The central approximation required for Eq. (1) is that of stationary inversion. Previous work by Haken [9] argued that SIA holds for the MB equations when γ ≪ γ ⊥ , ∆ , where γ is the relaxation rate of the inversion and ∆ is the frequency difference of lasing modes.
For a typical semiconductor laser γ ⊥ ≃ 10 12 − 10 13 s −1 and γ ≃ 10 8 − 10 9 s −1 , or equivalently γ /γ ⊥ ≃ 10 −3 −10 −5 . For the microcavity edge-emitting laser we are modeling we took γ ⊥ , ∆ ∼ 1 and γ = 10 −3 , and we found the excellent agreement shown in Figs. 1a, 1b. In addition, direct analysis of the inversion vs. time obtained from the MB simulations confirm very weak timedependence in the steady state, justifying the use of SIA. The previous work did not develop a systematic theory in which γ /γ ⊥ , γ /∆ appear as small parameters in the lasing equations, allowing perturbative treatment of corrections to the SIA; we are now able to do this within the AI formalism.
Note first that in Eq. (1) the electric field is measured in units e c ∼ √ γ but, unlike γ ⊥ , γ does not appear explicitly. Hence the solutions of Eq. (1) depend on γ only through this scale factor, and Eq. (1) makes the strong prediction of a universal overall scaling of the field intensities: |E(x)| 2 ∼ γ when dimensions are restored. The perturbative corrections to Eq. (1) are obtained by including the leading effects of the beating terms between the different lasing modes which lead to time-dependence of the inversion at multiples of the beat frequencies. These population oscillations non-linearly mix with the electric field and polarization to generate all harmonics of the beat frequencies in principle, but the multimode approximation assumes all the newly generated fourier components of the fields are negligible. The leading correction to this approximation is to evaluate the effect of the lowest sidebands of population oscillation on the polarization at the lasing frequencies and on the static part of the inversion, both of which will enter Eq. (1). For simplicity we present a sketch of the correction to Eq. (1) in the two-mode regime; details and the straightforward generalization to more modes will be given elsewhere. We start by writing the electric field and the polarization in the multiperiodic forms e(x x x,t) = , and allow for the first two side-bands at the beat frequency ∆ = k 1 − k 2 in the inversion, so that the total inver- where the dimensionless function f (k 1 , k 2 ) = −(i + ∆/(2γ ⊥ ))/(1 +k 1k2 − i∆/γ ⊥ ),k ν = (k ν − k a )/γ ⊥ and the fields are not yet measured in units of e c . The component d + will mix with the field Ψ 2 to yield a contribution to the polarization p 1 at frequency k 1 , ( and similarly d − and Ψ 1 mix to contribute to p 2 ), where p (1) 2 (x x x) is obtained by interchanging subscripts 1, 2. The correction to the AI formalism is the term in the numerator explicitly proportional to the small parameter γ /∆. However, having found a correction to the polarizations p 1 , p 2 we must then evaluate their contribution to the static inversion. We find where the dimensionless function g(k 1 , k 2 ) = (2 +k 2 1 +k 2 2 )(1 −k 1k2 )/[(1 +k 2 1 )(1 +k 2 2 )] 2 and now and below electric fields have been scaled by e c . Note that the correction term in Eq. (4) is explicitly proportional to the second small parameter, γ /γ ⊥ . For our simulations γ ⊥ ≈ ∆ and the functions f (k 1 , k 2 ), g(k 1 , k 2 ) are order unity. The full correction to the non-linear polarization in Eq.
. Eq. (5) predicts corrections to the universal behavior, |E(x x x)| 2 ∼ γ , found in Fig. 1. There is no correction to the first mode below the second threshold as the correction terms all vanish (there needs to be two modes to have beats). However, the theory predicts a non-trivial correction to the threshold of the second mode. Note that the correction to the numerator in Eq. (5) does not vanish below the second threshold but contributes self-consistently to that threshold. This correction can be regarded as modifying the dielectric function of the microcavity to take ; the effective dielectric function then becomes complex and varying in space according to the intensity of the first mode. This in turn changes the threshold for the second mode. If the modes are on opposite sides of the atomic frequency and k 2 < k 1 , the imaginary part of the effective index is always amplifying and tends to decrease the second threshold; we find this effect dominates over the change in the real part and increasing γ uniformly decreases the thresholds. The opposite effect is possible and observed in other cases we have studied (not shown). In Fig. 3 we show the results of MB simulations as γ is varied from 0.001 to 0.1 (with γ ⊥ = 4). Note that the universal behavior (in units scaled by γ ) is well obeyed until γ = 0.1, encompassing most lasers of interest. The qualitative effect predicted by the perturbation theory is clearly seen, the higher thresholds are reduced as γ is increased. The effect is small for the 2nd threshold but large for the third as we expect as the corrections scale with the product of the intensities of lower modes. The inset to Fig. 3 shows that the perturbation theory for the third threshold (a suitable generalization of the two-mode Eq. (5)) yields semiquantitative agreement with the threshold shifts found from the simulations. Detailed comparisons between multimode perturbation theory and simulations above threshold will be given elsewhere. Note that the simulations also find additional modes turning on for γ = 0.1 but their intensities are very small and they are not shown in Fig. 3.  Fig. 1a. The inset shows the shifts of the third threshold as a function of γ . The perturbation theory (squares, with the line to guide the eyes) predicts semiquantitatively the decrease of the threshold as γ increases found in the MB simulations. The MB threshold is not sharp and we add an error bar to denote the size of the transition region.
In conclusion, for the case studied, the recently developed ab initio self-consistent laser theory without the slowly varying envelope approximation presents a very accurate solution of the steady-state Maxwell-Bloch semiclassical lasing equations without solving the time-dependent problem. Third order treatments fail badly and our infinite order treatment is essential. The theory is a well-controlled expansion in the small parameters γ /∆, γ /γ ⊥ and leading corrections in these small parameters can be evaluated and understood qualitatively. The AI method is flexible and more convenient than MB simulations for complex two-dimensional and threedimensional microcavities and provides more physical insight into the lasing properties. The accuracy of the method suggests it can be useful in the analysis and design of novel laser systems.