The Maxwell-Bloch (MB) equations provide the foundation of semiclassical laser theory [1] and are the simplest description which captures the full space-dependent non-linear behavior of a laser. These time-dependent equations can be simulated to determine the stationary lasing 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–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 hole-burning) 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 micro-cavity edge emitting laser [2, 3] consisting of a perfect mirror at the origin and a dielectric region of uniform index n ₀ 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 rotating-wave 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 ₀ k _a L=60 corresponding to roughly ten wavelengths of radiation inside the cavity, approaching the micro-cavity limit. For the case studied here, we find noticeable discrepancies between MB solutions and the AISC theory of Refs. [2–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).

 

Fig. 1. Modal intensities as functions of the pump strength D ₀ in a one-dimensional micro-cavity edge emitting laser of γ _⊥=4.0 and γ _‖=0.001. (a) n=1.5, k _a L=40. (b) n=3, k _a L=20. Square data points are the result of MB time-dependent simulations; solid lines are the result of time-independent ab initio calculations (AISC) of Eq. (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 ₀=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.

Download Full Size | PDF

The generalization was as follows. Again stationary periodic solutions are assumed for the electric field E(x,t)=2Re[e(x,t)]=2Re[∑_µ Ψ _µ(x)exp(-ik _µ t)] and for the polarization fields, which oscillate at unknown lasing frequencies, k _µ. The spatial variation of the field amplitude Ψ_µ(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) 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 AISC/MB equations for the mode functions Ψ_µ and the frequencies k _µ are:

Here G(x,x ^′;k _µ) is the Green function of the open cavity, Γ_ν=Γ(k _ν) is the gain profile evaluated at k _ν, D ₀(x)=D ₀(1+d ₀(x)) is the spatial pump, and ε(x)=n ²(x) is the dielectric function of the cavity (for the micro-cavity edge emitting laser n(x)=n ₀ and the pump is assumed uniform (d ₀(x)=0)). Electric field and pump strength are measured in units $e_c=\frac{\overline{h}\sqrt{{\gamma }_{\perp }{\gamma }_{\parallel }}}{2g,D_{0c}}=\frac{4\pi k_a^2{g}^2}{\overline{h}{\gamma }_{\perp }}$ , 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 ² _µ/k ² _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 ₀=1.5, 3, finding remarkable agreement between MB and AISC 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 AISC approach is used with the SVEA then significant discrepancies are found, for example for n ₀=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 theMB 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 |Ψ_ν|². 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 ₀=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. 2. Modal intensities as functions of the pump strength D ₀ in a one-dimensional micro-cavity edge emitting laser of n=3, k _a L=20, γ _⊥=4.0 and γ‖=0.001; the solid lines and data points are the same as in 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 AISC 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 ₀=10, while both the MB and AISC show only four. Right inset just shows the same data on a larger vertical scale.

Download Full Size | PDF

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¹²-10¹³ s ^-1 and γ _‖≃10⁸-10⁹ s ^-1, or equivalently γ _‖/γ _⊥≃10^-3-10^-5. For the micro-cavity 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 time-dependence 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 AISC formalism.

Note first that in Eq. (1) the electric field is measured in units $e_c~\sqrt{r_{\left|\right|}}$ 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)|²~γ _‖ 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,t)= 2Re[∑² _µ=1Ψ_µ (x)exp(-ik _µ t)] and p(x,t)=2Re[∑² _µ=1pµ(x)exp(-ik _µ t)], and allow for the first two side-bands at the beat frequency Δ=k ₁-k ₂ in the inversion, so that the total inversion is D(x,t)=D _e(x)+d+(x)exp(-iΔt)+d-(x)exp(+iΔt) where the real quantity D _s(x) is the time-independent part of the inversion and d+(x)=d-(x)*. This ansatz is inserted into the Bloch equation for the inversion [2]; solving for the component of the inversion equation which oscillates at exp(-iΔt) relates d+(x) to the product of the field and polarization and then substitution of the zeroth order result for the polarization in terms of the zeroth order static inversion D ⁽⁰⁾ _s gives

where the dimensionless function f (k ₁,k ₂)=-(i+Δ/(2γ _⊥))/(1+k̃₁ k̃₂ -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 Ψ₂ to yield a contribution to the polarization p ₁ at frequency k ₁, (and similarly d- and Ψ₁ mix to contribute to p ₂),

where p ⁽¹⁾ ₂(x) is obtained by interchanging subscripts 1, 2. The correction to the AISC formalism is the term in the numerator explicitly proportional to the small parameter γ _‖/Δ. However, having found a correction to the polarizations p ₁, p ₂ we must then evaluate their contribution to the static inversion. We find

where the dimensionless function g(k ₁,k ₂)=(2+k̃² ₁+k̃² ₂)(1-k̃₁ k̃₂)/[(1+k̃² ₁)(1+k̃² ₂)]² 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 ₁,k ₂),g(k ₁,k ₂) are order unity. The full correction to the non-linear polarization in Eq. (3) is obtained by replacing D ⁽⁰⁾ _s with D _s of Eq. (4). The corrected polarization leads to corrected version of Eq. (1) of the AISC theory:

Ψ₁(x) satisfies the same equation with the subscripts 1 and 2 interchanged.

Equation (5) predicts corrections to the universal behavior, |E(x)|²~γ _‖, 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 micro-cavity to take the form $\epsilon \text{'}(\mathit{x})=\frac{\epsilon (\mathit{x})}{[1+\frac{{\gamma }_{\left|\right|}}{\Delta }f(k_2,k_1)|{\psi }_1(\mathit{x}){|}^2]}$ ; 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 ₂<k ₁, 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. 3. Modal intensities for the micro-cavity edge emitting laser of Fig. 1a as γ _‖ is varied (n=1.5, k _a L=20, γ _⊨=4.0 and mode-spacing Δ≈1.8. Solid lines are AISC results from Fig. 1a (with γ _‖=0.001); dashed lines are γ _‖=0.01 and dot-dashed are γ _‖=0.1. The color scheme is the same as in 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.

Download Full Size | PDF

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 AISC method depends on an iterative solution of Eq. (1). The numerical method involves three modules: a complex eigenvalue solver to find the basis set of biorthogonal states in which to expand the solutions, Ψ_µ(x), construction of threshold matrices to find the modal thresholds from a linear analysis, and solution of the non-linear eigenvalue problem posed by Eq. (1) itself. Some details are given in ref. [4], and particularly in the on-line supporting materials. We have not found standard numerical packages to implement this algorithm, nor has our method been systematically optimized. In contrast the MB code is constructed out of standard numerical components. In our current versions, generating the modal intensities for one value of D ₀ as in Fig. 1a takes roughly 60 times longer with the MB code as compared to the AISC code, and the AISC code works roughly 130 times more efficiently to capture the entire range of D ₀.

However the qualitative advantage of the AISC approach from the numerical (as opposed to the conceptual) point of view arises for more realistic two-dimensional and three-dimensional laser structures. Note that the reason we can employ the stationary inversion approximation in deriving Eq. (1) is that there is a separation of time scales in typical lasers. The approach to steady-state is controlled by the long time scale, γ _‖ ^-1, whereas the field oscillates at the rapid rate, ~k _a. Full wave MB simulations must bridge these two time scales and hence are quite time-consuming. That is why the SVEA is employed, to transform to the Schrodinger-Bloch equations [11]. However our results demonstrate that the SVEA is not very accurate. Full wave MB simulations in two and three dimensions over the full range of time scales necessary to describe multi-mode lasing are probably not yet feasible, and have not yet been done. In contrast, two and three dimensional AISC calculations are quite feasible; two-dimensional random laser simulations have already been performed [4], and two-dimensional micro-cavity laser simulations are in progress [13]. Therefore the AISC method can give access to accurate laser simulations which are currently not possible by any other method. This method also provides more physical insight into lasing properties. The accuracy of the method suggests it can be useful in the analysis and design of novel laser systems.