Steady-state ab initio laser theory for complex gain media

We derive and test a generalization of Steady-State Ab Initio Laser Theory (SALT) to treat complex gain media. The generalized theory (C-SALT) is able to treat atomic and molecular gain media with diffusion and multiple lasing transitions, and semiconductor gain media in the free carrier approximation including fully the effect of Pauli blocking. The key assumption of the theory is stationarity of the level populations, which leads to coupled self-consistent equations for the populations and the lasing modes that fully include the effects of openness and non-linear spatial hole-burning. These equations can be solved efficiently for the steady-state lasing properties by a similar iteration procedure as in SALT, where a static gain medium with a single transition is assumed. The theory is tested by comparison to much less efficient Finite Difference Time Domain (FDTD) methods and excellent agreement is found. Using C-SALT to analyze the effects of varying gain diffusion constant we demonstrate a cross-over between the regime of strong spatial hole burning with multimode lasing to a regime of negligible spatial hole burning, leading to gain-clamping, and single mode lasing. The effect of spatially inhomogeneous pumping combined with diffusion is also studied and a relevant length scale for spatial inhomogeneity to persist under these conditions is determined. For the semiconductor gain model, we demonstrate the frequency shift due to Pauli blocking as the pumping strength changes.


Introduction
Semiclassical laser theory describes an inverted gain medium coupled to classical electromagnetic radiation described by Maxwell's equation. It is sufficient to describe all properties of interest in lasers (e.g. the number, spatial distribution, frequencies, output power and turn-on thresholds for the lasing modes) except for those which arise from quantum fluctuations (such as the linewidth). These equations are the main tool for modeling and design of lasers. The full semiclassical laser equations are coupled non-linear partial differential equations in space and time, describing the coupling of the polarization and level populations in the gain medium to light, and are not easily solved even numerically. In most treatments the spatial variation of the quantities in these equations is either completely neglected or addressed through an approximate expansion in closed cavity states which are taken to be orthogonal, hence turning the resulting equations into coupled ordinary non-linear differential equations in time. Furthermore, the gain medium is often described as a simple two-level atomic system and modeled using the Bloch equations.
The computational complexity of the laser equations is substantially increased when it becomes necessary to include many atomic levels and multiple lasing transitions in the mathematical description, as is necessary for many modern laser systems with complex gain media. Examples of such systems are standard semiconductor lasers [1,2], quantum cascade lasers [3,4], and rotationally excited gas lasers [5]. For example, though the band structure of the semiconductor gain medium can be approximated as a series of two level atomic transitions, multiple transitions are required to represent the effects of Pauli blocking [6][7][8]. Cascaded-transition quantum cascade lasers are specifically designed with two lasing transitions to operate at longer wavelengths [9,10]. Finally, in these lasers, and a number of others, the carriers or the gain atoms are able to diffuse after injection or excitation, which tends to smooth out the effects of spatial hole burning, and strongly affect the number of lasing modes for a given pump. This critical aspect of laser modeling has received relatively little attention in the literature [4,[11][12][13] and has been treated only with great simplifications. One exception to this is the work of Böhringer and Hess, who provide a more accurate account of the diffusion of gain carriers in semiconductor media, but within the context of constructing a time-domain simulation algorithm [7]. It is possible to include gain diffusion in brute force (FDTD) simulations, but this adds significant computational overhead that makes it difficult to simulate cavities in multiple dimensions [14,15]. While such simulations are done in space and time they often are only aimed at finding the steady-state lasing properties. In the current work we demonstrate how the effects of many levels, multiple lasing transitions, and gain diffusion can all be included in a theoretical framework which leads to a computationally tractable algorithm that directly calculates the steady-state properties. The algorithm is efficient enough to treat both complex gain media and complicated resonator geometries. Below it is tested against FDTD results for simple one-dimensional cavities without diffusion and excellent agreement is found. It is then applied to two-dimensional deformed disk lasers, giving results which explain certain experimental observations. Calculations of this type are beyond the scope of previous approaches to the semiclassical laser equations.
Our new framework is a generalization of the theoretical approach introduced by Türeci, Stone, et al. now known as Steady-state Ab initio Laser Theory (SALT) [16][17][18][19]. This theory employs only a single important approximation, the stationary inversion approximation (SIA), to treat multimode lasing, and is essentially exact for treating single-mode lasing; and it includes fully the spatial degrees of freedom of the fields and the openness of the laser system. The resulting SALT equations are frequency domain wave equations for the lasing modes, coupled non-linearly through the spatially varying cross gain saturation. These equations can be solved efficiently for the steady-state properties of laser cavities of arbitrary complexity, in any number of dimensions. (The vector character of the lasing modes in three dimensions can also be treated [20], but in this work we will analyze only scalar limits of the electromagnetic field equations). This approach allows a clearer and partially analytic understanding of the lasing solutions, and has already led to a number of new discoveries, such as mode frustration in partially pumped cavities [21], control of emission properties of random lasers through selective pumping [22][23][24][25][26], and, through a quantum generalization, a more general form of the Schawlow-Townes linewidth formula [27][28][29]. For single transition gain media without diffusion, SALT has also been shown to give excellent agreement with FDTD simulations at a substantially reduced computational cost [30,31]. In the case of a single transition the gain susceptibility for N-level lasers can be explicitly calculated in terms of the stationary inversion profile and the lasing fields and inserted into the wave equation so as to require only a single set of coupled electric field equations [31]. It has been shown that in this case the N-level model can be reduced to the two-level Maxwell-Bloch model with renormalized relaxation and pump parameters [31]. Such a simplification is not possible with multiple lasing transitions, because it is no longer possible to associate a particular lasing mode with a particular pair of inverted levels, as we will show below. Thus a significant generalization of the theory is required, and it will no longer be possible to map the full solution to an effective Maxwell-Bloch model. However we will show below that a more general stationarity assumption (The Stationary Population Approximation, SPA) leads instead to coupled sets of field and population equations that can be solved iteratively almost as efficiently as the SALT equations. The generalized theory, which we call complex-SALT (C-SALT), allows us to treat steady-state lasing with (i) an arbitrary number of atomic levels and lasing transitions, and (ii) gain diffusion, and, just as in SALT, the full effects of non-linear spatial hole-burning are treated for an arbitrary cavity geometry.
The ability of C-SALT to treat multimode lasing in the presence of diffusion allows us to demonstrate explicitly a fundamental aspect of laser physics: the competition between spatial hole-burning and gain clamping/saturation. For a completely spatially uniform lasing mode (e.g. in a ring laser), the mode which has the highest gain, after it reaches threshold, clamps the gain allowing no further modes to turn on [32]. For the typical spatially non-uniform (e.g. standing wave) mode, the effective gain varies in space (even with uniform pumping) because where the existing lasing modes have high intensity the gain is highly saturated and reduced, whereas near nodal regions of the lasing modes other modes can experience higher gain and turn on. This is the standard spatial hole-burning effect, which allows multimode lasing. However, if there is gain diffusion, then excited gain atoms/carriers tend to diffuse into the regions of low field intensity, smoothing out the effective gain profile. Thus as a function of the strength of the diffusion the laser will transition from highly multimode lasing to a more uniform single or few mode lasing state. This picture is well known, but it has not been possible to study it quantitatively in previous approaches, as we do below using the full C-SALT equations.
We identify two distinct physical situations in which carrier diffusion can have a substantial effect. The first,as just discussed, is when the carrier diffusion competes with spatial holeburning in a uniformly pumped cavity, leading to a transition between multimode lasing to a gain-clamping regime in which only a single lasing mode reaches threshold. We define a crossover scale in the diffusion rate for this clamping to occur and test it in simulations. The second situation is when carrier diffusion competes with non-uniform pumping so as to mitigate the spatial selection effects which occur in the absence of diffusion. A different crossover scale is found to control this effect. We are not aware of prior quantitative theoretical studies of either of these effects.
The remainder of this paper is organized as follows. In section II we introduce the semiclassical lasing equations and analyze the most general conditions under which they lead to steady-state multimode lasing. Section III derives the equations for the steady-state (A) for a single transition, leading to the SALT equations (B),(C) for multiple lasing transitions and gain diffusion, leading to the C-SALT equations, which are the main result of this work. Section IV (A) demonstrates quantitative agreement between C-SALT and FDTD simulations for the case of multiple lasing transitions without diffusion. (B) demonstrates the transition between spatial hole-burning and gain clamping in the presence of of variable diffusion. In section V a discussion of the Stationary Population Approximation is presented. Section VI demonstrates how C-SALT can be used to simulate lasers described by a free-carrier semiconductor model. Finally, a conclusion is given in section VII.

Overview of the semiclassical lasing equations
In this section we present the derivation of the equations of motion for the density matrix elements for an atomic gain medium with an arbitrary number of levels and lasing transitions. The basic steps here are well known but we want to take the opportunity to emphasize certain key approximations which set the limits of the standard semiclassical approximation and also to define our notation. As noted above, the semiclassical lasing equations neglect the operator nature of the electromagnetic fields, so that the lasing fields are determined by Maxwell's wave equation (coupled to a quantum gain medium), where the effect of the passive cavity is described by the linear cavity dielectric function, ε c (x), which in general varies in space and frequency (although typically we will neglect the frequency variation). The amplifying response of the gain medium is described by the non-linear gain polarization, P g , which acts as a source for this equation and includes contributions from all of the lasing transitions of every gain atom in the cavity. At this stage we will assume that the gain medium arises from independent identical atomic, molecular or defect centers; the case of semiconductor gain media with band excitations acting as gain centers will described in Section 6 below. The gain polarization can be expressed as in which the index α runs over all of the gain atoms in the cavity at locations x (α) , e is the charge of an electron, and the M × M density matrix of atom α is denoted byρ (α) , where M is the number of atomic levels involved in the lasing process. Among the M levels a subset of them will contribute to lasing and support lasing transitions; the others will simply be part of the downward cascade of electronic excitations involved in the pumping and emission in steady-state. Lasing transitions will arise from pairs of level which have sufficiently large polarization due to their population inversion to contribute substantially to a lasing line at a nearby frequency. These pairs have a dipole moment, The dipole moment is always zero for n = m due to spatial symmetry. By assuming that the full polarization can be expressed in terms of the density matrix for individual atoms, we are ignoring interatomic coherence effects which are quite small for conventional lasers (but not e.g. for polariton lasers).
To complete the semiclassical lasing equations, one must consider the quantum equations of motion for the average polarization of the gain medium, where we have now assumed identical "atoms" and written the density matrix elements as ρ nm = n|ρ|m , and we initially assume a fixed density of gain atoms N(x). The evolution of the density matrix can be found from the Heisenberg equation of motion, where the atomic Hamiltonian H 0 |n = E n |n , and the interaction Hamiltonian can be written as H I = ex · E(x,t). Upon evaluating the commutator and simplifying, the evolution of the density matrix elements can be re-written as where From this equation, we can see that off-diagonal density matrix elements, which determine the gain polarization, can couple to one another within manifolds of atomic transitions. If we include all terms in the equation of motion, then the time evolution of a specific off-diagonal element, ρ nm , will depend not only on the level populations (ρ nn , ρ mm ), but also on other offdiagonal elements, i.e. on the polarization of other transitions. If this is the case one cannot arrive at lasing equations of the standard form and one cannot define the polarization associated with a specific transition. However, physically these off-diagonal terms correspond to coherent multiple excitation, leading to effects such as electrically induced transparency, and inversionless lasing. Conventional lasers do not typically operate in this regime. To reach this regime, the non-radiative relaxation rates between non-lasing transitions must be of similar order to those between lasing levels, which makes it difficult for the gain medium to build up the necessary inversion to lase [33]. As such, we will assume that we are in the weakly coupled polarization regime and thus that off-diagonal density matrix elements depend only on the level populations of that specific pair of levels. The equation of motion for the off-diagonal elements for such a pair is where we have now added the effect of environmental dephasing on the gain atoms in the standard manner in terms of a transverse relaxation/dephasing rate γ ⊥,nm . With the above assumption, the total gain polarization can now be broken up into N T constituent polarizations of each lasing transition, where each transition, previously labeled by the pair of levels, n, m, has been relabeled with a transition index j, where p + j is the positive frequency part of p j = p + j + p − j , and by definition ω nm > 0.
From the form of the gain polarization as a transition sum, Eq. (8), it is clear that every constituent polarization in the gain medium contributes to the total source term in the wave equation, Eq. (1). Strictly speaking, it is thus impossible to say that one portion of the electric field is driven by only a single transition if multiple transitions are present in the gain medium, though practically there are many cases where these transitions are so well separated in frequency that this is effectively the case. Each constituent polarization, from Eq. (7), obeys its own equation of motion, in which the properties of the jth constituent polarization are given in terms of the dephasing rate, γ ⊥, j , the atomic transition frequency, ω a, j , the constituent inversion, , and the dipole matrix element, θ θ θ j . The equation of motion for the density matrix, Eq. (6), also determines the evolution of the populations in each atomic level, given by the diagonal elements of the density matrix, As can be seen here, atomic level populations couple only to the constituent polarizations with which they share a transition, and atomic levels which are not a part of any lasing transition do not appear at this stage on the right hand side of the equation for the populations. However typically all levels are coupled non-radiatively through other degrees of freedom ("the bath") and these effects need to be included phenomenologically in the standard manner, leading to Here γ nm represents either the non-radiative decay rate between a higher level m and a lower level n or pumping rate from lower level m to higher level n, ρ n ≡ N(x)ρ nn is the total number of electrons in level n of atoms at location x. ξ n, j represents the relationship between the population of level, n and a given lasing transition, j; ξ n, j = 1 if n is the upper level of the transition, ξ n, j = −1 if n is the lower level of the transition, and ξ n, j = 0 if n is not involved in that lasing transition.
Thus the full set of semiclassical lasing equations are Eqs. (1), (8), (10), and (12), which define the wave equation for the electric field, the total polarization in terms of the constituent polarizations, the equation of motion for each constituent polarization, and the equation of motion for the populations in each atomic level respectively. Since now the time evolution of the off-diagonal elements of the density matrix are contained in the polarization equations, we will henceforth represent the populations (diagonal elements) in terms of a density vector with components ρ n (x).

Two-level gain medium
In the limit of a two-level gain medium with only a single lasing transition, Eqs. (8), (10), and (12) are simplified to where γ = (γ 12 + γ 21 ) is the relaxation rate of the inversion and is the inversion in the absence of an electric field. These equations, Eqs. (13) and (14), along with the wave equation, Eq. (1), comprise the Maxwell-Bloch equations, the most basic version of semiclassical laser theory. As noted above, it is these equations which are solved directly for the steady-state properties using SALT. For single mode lasing the SALT solution is essentially exact; for multimode lasing it requires the stationary inversion approximation (SIA), ∂ t d = 0. This approximation holds if the multimode beating terms which drive time-dependence of the populations are negligible. The stationary inversion approximation (SIA), requires two conditions to be valid. First, that the relaxation rate of the inversion, γ , be small compared to the modal spacing, Δ, and the dephasing rate of the polarization, γ Δ, γ ⊥ . This condition is usually met for microlasers, for example a Fabry-Perot laser cavity typically requires L ≤ 100μm for Δ γ [34]. The second condition is that the mode spacing must be well separated from the relaxation oscillation frequency, so as to avoid resonantly driving relaxation fluctuations which could destabilize the multimode solution. The relaxation oscillation frequency is ω r ∼ √ κγ [35], where κ is the field decay rate of the cavity. For microcavities, κ ≤ Δ, so ω r ≤ Δγ < Δ. As such, the SIA will hold and the multimode solution found by SALT will be correct when γ Δ, γ ⊥ [36].
Essentially the same conditions will be required for C-SALT to describe multimode lasing, although, as for SALT, C-SALT solutions should be accurate in general for the single-mode case.
Continuing our brief review of SALT, we now simplify our analysis by treating slab or twodimensional geometries for which the electric fields in the transverse magnetic (TM) modes, can be taken to be a scalar, E → E), noting that the treatment discussed here is still completely applicable in geometries for which the fields must be treated as vectors [20]. Having done so, we make a multimode ansatz, stating that the electric field and polarization can be broken up into N L components with distinct frequencies representing each lasing mode, where the plus superscript denotes the positive frequency component of the field, E = 2Re[E + ], and Ψ μ (x) and p μ (x) are the spatial profiles of the electric field and corresponding polarization of the lasing mode with frequency ω μ . The multimode ansatz allows us to match frequency components of the electric and polarization fields through Eq. (13), Upon inserting this definition into Eqs. (14) and (1) and making the rotating wave approximation (RWA), which is well satisfied for all lasers of interest, one recovers the coupled set of SALT equations. As noted above, these take the form of a set of wave equations, one for each lasing mode of the cavity, which are coupled through the non-linear hole-burning interaction in the gain medium [16,17,19], where is the Lorentzian gain curve evaluated at ω ν . To solve the SALT equations, Eqs. (19) and (20), simultaneously, each of the lasing modes is expanded over a basis, where the basis functions f n are fixed, (but possibly dependent upon the lasing frequency ω μ ), and the expansion coefficients a n are found via a non-linear solution algorithm. There are at least two useful choices for the basis set. The first, and most developed to this point, is to use a basis set of purely outgoing non-hermitian states which are termed the Constant Flux (CF) states [16,19,37]. More recently it has been shown that one can simply use a position space basis [20] combined with a perfectly matched layer (PML) to implement the outgoing boundary condition. We will not go into the details of these computational approaches here as either can be used with the C-SALT equations which are our main focus. Instead we refer the interested reader to the cited references; the computations presented here are based on the CF approach. SALT assumes only a single lasing transition, and the above equations are then typically presented in natural, scaled units for electric field and inversion, , that involve the relaxation rates associated with that transition. With multiple transitions with differing relaxation rates, all contributing to the lasing lines, there will be no such natural scaling for C-SALT. As noted above, previous work has shown that the steadystate equations for N-level lasing with a single transition can be exactly mapped onto the SALT equations with renormalized relaxation parameters, so that computationally the N-level case is negligibly different from the two-level Maxwell-Bloch case [31].

Multiple lasing transitions
When multiple lasing transitions are present, the multimode ansatz is expanded to include constituent polarizations, which still allows one to match frequency components for each constituent polarization using Eq. (10) to find To derive the C-SALT equations, this solution for the constituent polarization is inserted into the equation of motion for the atomic levels, Eq. (12), and the RWA is again made, resulting in To proceed, we first rewrite this equation in terms of a M × M population density vector at each position in the cavity, ρ ρ ρ(x), whose components are the atomic level populations ρ n (x). Additionally, we make a slightly stronger version of the SIA, the stationary population approximation (SPA), which states that the beating between different lasing modes does not lead to significant time-dependence in the level populations, so ∂ t ρ ρ ρ ≈ 0. While a more detailed discussion of SPA will be given in section 5, it is worth noting here that it is valid for most laser systems of interest, i.e. those with fast decays into and out of the lasing levels, while the upper lasing level of each transition is metastable and thus long lived. These fast relaxation rates need not be small compared to Δ and γ ⊥ , only the relaxation rates of the metastable upper lasing transitions need be slow in this sense, a condition which is numerically tested and confirmed by the FDTD simulations below. Thus, the above equation can be rewritten as where R is a matrix containing information about the pump and decay rates, and Ξ j is a matrix containing information about which level populations are coupled to the partial polarizations and constitute the jth lasing transition. The full forms of these matrices are given in Appendix A.
Equation (25) is a homogeneous equation satisfied by the "vector" of atomic populations at each point in space. It requires knowledge of the lasing modes to be solved and hence will need to be solved simultaneously with the electric field equations. In addition, in the absence of gain diffusion, the total number of gain atoms at each point in space, N(x), is externally fixed, and is a given of the problem. Hence the homogeneous Eq. (25) is to be solved subject to the normalization condition which uniquely determines the level population vector. It is convenient to incorporate this normalization condition directly into Eq. (25) by defining a matrix B and an M-component total number vector N(x) such that where neither the matrix B, nor the vector N(x) are uniquely defined, but must be chosen to represent Eq. (26). The normalization can then be inserted into Eq. (25), resulting in To incorporate multiple lasing transitions into the SALT formalism to recover the C-SALT equations, all that must be altered is the equation for the electric susceptibility, Eq. (20), not the wave equation itself. Using Eq. (23), the electric susceptibility can be written as where m is determined from Eq. (28). The main difference here is now the atomic population densities cannot be directly inserted into the wave equation through the use of a scalar inversion equation. Instead we must simultaneously solve the wave equation, Eq. (19), and the population equation, Eq. (28). Using a non-linear iteration algorithm to solve the problem numerically, one inserts an initial guess for the field profiles, {Ψ μ (x)}, and uses these to solve for the full spatial profile of the atomic population densities. This result for the population densities generates a guess for the susceptibility, which can be inserted into the wave equation and iterated back and forth to self-consistency.
There is one other important difference between Eqs. (20) and (29) in the case of a "partially pumped" laser, in which pumping is applied to only a subset of the regions containing gain media [22][23][24][25][26]. One needs to distinguish physically between two situations: 1) Having a (given) spatial density of gain atoms, which can be non-uniform, and hence will lead to non-uniform but always positive gain under conditions of uniform spatial pumping. This is taken into account from the specification of N(x) under the assumed conditions of spatially uniform pumping. 2) Spatial non-uniform pumping "partial pumping" in which there is a uniform distribution of gain atoms, which are not equally pumped. This would appear in our formalism as a variation with spatial position of the pumping rate parameters in the R matrix of Eq. (28). In this case non-pumped regions would act as absorbers for laser light, which would automatically be taken into account in terms of a spatial variation in the susceptibility function, χ g (x), which would now have an absorbing form in those unpumped regions.

Gain diffusion
The formalism developed in Section 3.2 can be extended to include gain diffusion, a phenomenon found in many types of gain media. A term representing this effect can be added to Eq. (24), resulting in where D n is the longitudinal diffusion coefficient for the atomic level |n . Despite the similarities of Eqs. (24) and (30), there is one important difference, namely that the atomic populations at each spatial location are now coupled together. Thus, when diffusion is present, the population density vector ρ ρ ρ is an M × P dimensional vector whose components are atomic populations at each of P discretized spatial locations. The stationary population approximation can still be made, resulting in a generalized homogeneous equation for the population density vector, where D is the matrix of longitudinal diffusion coefficients at each spatial location, and the other two matrices R, and Ξ have also now been similarly expanded over the position basis as well.
To correctly normalize Eq. (31), we note that one consequence of SPA is ∂ t ∑ n ρ n = 0, and when performing this sum on Eq. (30), one finds that the total population density is homogeneous in the steady state in the presence of diffusion, Furthermore, the walls of the cavity prevent any flux of gain atoms across its borders, which is represented by the Neumann boundary condition ∂ x ρ n | x=0,L = 0.
This, together with Eq. (32), yields the normalization This is a simple restatement of the fact that a diffusive gain medium in a passive cavity will be evenly distributed in the steady state, unlike a non-diffusive gain medium where the density of gain atoms can fluctuate depending on their distribution. Note however, that while the diffusion condition implies that the number of "atoms" will be uniform in space, their excitation distribution will not be. Spatial hole-burning due to the spatial variation of the lasing modes above threshold will lead to a non-uniform gain even in the presence of uniform pumping and diffusion, and this effect is still captured by the theory. We will see below that these two effects compete; as the diffusion coefficient of the medium is increased, the excitation distribution will become more uniform in steady-state, counteracting the effects of spatial hole-burning.
Again, the normalization requirement can now be expressed as where the matrix B and vector N are likewise expanded over the spatial basis. Inserting Eq. (35) into Eq. (31) yields the generalized level population equation The solution now proceeds in the same way as in Section 3.2. Knowing the atomic level populations ρ ρ ρ, we can use Eq. (29) to solve for the susceptibility, which is then inserted into the wave equation, Eq. (19). The problem thus reduces to a set of differential equations, one per lasing mode, coupled through the generalized level population equation, Eq. (36). The latter is now a much bigger matrix equation, MP × MP, as opposed to M × M, increasing the computational cost, but still keeping it within a manageable range, even for two-dimensional lasers (see below). These two generalizations, Eqs. (28) and (36), along with their coupling to SALT, Eqs. (29) and (19), are the main results of this paper. They extend the capabilities of SALT to cover many types of gain media, including newly-developed ones [5,10]. C-SALT can thus be used as an efficient tool for the design and study of devices in which gain diffusion is present alongside spatial hole-burning, a regime which is challenging for traditional numerical methods to handle [7,8,14,15].

Results
To perform a well controlled test of the stationary population approximation in C-SALT, for the case of multiple transitions without diffusion, we studied 1D microcavity lasers for which FDTD simulations (without diffusion) are tractable. We used a FDTD scheme similar to the one proposed by Bidégaray [38], altered to include multiple lasing transitions and additional atomic populations. The simulations were run for a total duration at least 40 times greater than the longest time scale in the model, to ensure a steady state was reached. FDTD simulations of Maxwell's equations coupled to atomic populations and polarizations do not make the steadystate assumption used in C-SALT, thus observing a stable output in FDTD serves as a test of the SPA used in C-SALT. The simulated laser cavity consists of a dielectric slab of background refractive index n = 1.5, with a perfectly-reflecting mirror on one side and an interface with air on the other. Distributed uniformly within the slab is a six-level gain medium, with two atomic transitions of slightly different frequencies and widths. For these set of simulations, the diffusion constant was set to zero in the C-SALT equations. The gain linewidths of the transitions overlap, as shown in the inset of Fig. 1, so each lasing mode receives significant gain contributions from both transitions. As shown in Fig. 1, excellent agreement is seen between C-SALT and FDTD simulations, both in predictions of modal intensity and frequency, thus quantitatively verifying the use of SPA. In these simulations, only the relaxation rates of the metastable upper lasing levels are small compared to Δ and γ ⊥, j ; the relaxation rates of other atomic levels are of the same order. This provides evidence for our earlier claim that the only rates which must be small when compared to Δ and γ ⊥, j are those of the metastable upper lasing state.
We now move to the case of lasers with multimode lasing and gain diffusion, for which FDTD is a very challenging computational effort, but which are relatively tractable using the C-SALT approach. In Fig. 2, we demonstrate how gain diffusion affects the transition from gainclamped single-mode lasing (which can be described by the simplest form of laser theory [32]) to multimode lasing (which is possible as a result of spatial hole-burning). The left panel of Fig. 2 shows C-SALT simulations for a two level atomic medium with three different values for the diffusion coefficient. The solid lines show the modal intensities of the medium without diffusion, and dotted and dot-dashed lines of the same color show the evolution of the modal intensities as the diffusion coefficient is increased. We observe, as expected, that increasing the diffusion coefficient postpones the transition from single-mode to multimode operation, by increasing the threshold of the second and higher-order lasing modes. In fact, for the largest diffusion coefficient the third lasing mode does not reach threshold within the pump values that were simulated. The right panel of Fig. 2 shows the inversion of the gain medium inside the cavity as a function of position within the cavity for a given pump value, d 0 = 0.345, a point The results in Fig. 2 make intuitive sense: an increased diffusion coefficient spatially "smooths out" the population inversion of the first lasing mode; this supplies the first mode with more gain, acting against the effects of spatial hole-burning. Specifically, the right panel of Fig. 2 shows that increasing the diffusion coefficient flattens the inversion close to the cavity mirror, where the lasing modes spatially overlap. Near the end facet of the cavity, the gain is already being used fairly uniformly, so the effects of diffusion are substantially reduced. These results can also be understood quantitatively by writing down a steady-state, effective inversion equation for an atomic gain medium, d, which can be done as there is only a single lasing transition [31], where d 0 (x) is the equilibrium value of the inversion density in the absence of both fields inside the cavity and diffusion. The final term on the right hand side can be identified as the rate of stimulated emission, which is spatially dependent and proportional to the intensity of the local electric field. Using the SPA, Eq. (37) can be solved for the inversion density, as Thus, for diffusion to be germane to the system, where k D is the wavevector associated with the scale of the inhomogeneity in the inversion. If the inversion is less than this, the gain atoms are unable to move very far before they either nonradiatively decay or undergo stimulated emission, preventing the diffusion from washing out the effects of the spatial inhomogeneity in the inversion. For spatial hole-burning, this variation in the inversion is on the order of the atomic transition wavelength, k D = 2k a , as the inversion oscillates twice as fast as the electric field [4]. This prediction is consistent with the numerical results shown in Fig. 2, where the onset of strong diffusion suppresses multimode operation and leads to gain clamped behavior. For partially pumped cavities [22][23][24][25][26], there is, in addition to the scale associated with smoothing out spatial hole-burning, another relevant scale for measuring the strength of diffusion. This is the scale at which the diffusion begins to overcome the spatially inhomogeneous pumping, an inhomogeneity on the scale of the entire cavity length. Quantitatively, this can still be understood from Eq. (40), except that now k D = 2π/L, where L is the length of the cavity, and we are assuming that the variation of the pumping is on the scale of the entire cavity. This leads to the criterion that when 4π 2 D/L 2 > γ + γ SE (I) diffusion will strongly reduce the effects of non-uniform pumping by allowing the inversion to penetrate into the non-pumped regions.
Both types of diffusion-induced transitions are demonstrated in Fig. 3, which shows the results of C-SALT simulations of a four-level gain medium with a single transition, which is only pumped in half the cavity. A two-level medium would not be suitable for this test as it would have strong absorption in the unpumped region. As the diffusion coefficient is increased, the first transition from the spatial hole-burning regime to the spatially-averaged gain saturation regime is observed, as well as the effects of gain clamping which increases the second lasing threshold. In this regime, the inversion does not penetrate into the unpumped region. However, as the diffusion coefficient is further increased, the inverted atoms are able to penetrate further into the unpumped region. For this example, the effective spatial scale of the partial pumping is the length of the cavity, L, by choice. As expected, we observe this second transition when 4π 2 D/L 2 γ ∼ 1 + γ SE (I)/γ . Once the first transition has occurred the gain is sufficiently clamped, even though it is not uniform over the entire cavity, and we only find single mode lasing. If one intentionally pumps non-uniformly on a smaller scale than L, that would decrease the diffusion coefficient needed for the second transition until in the case of wavelength scale non-uniform pumping the two transitions would coincide roughly.
To demonstrate the scalability of C-SALT to multiple dimensions we ran 2D TM simulations of quadrupole-shaped dielectric cavities, whose boundary is defined by the equation r(θ ) = r 0 (1 + ε cos(2θ )), where ε represents the degree of deformation from a circular cavity. Such cavities have been extensively studied in the context of wave chaos theory and experiments [37,39,40]. It was found that the spatial profile of the first lasing mode depends upon both the deformation and the pumping profile [37,39,40]. Here we show in Fig. 4 that for a quadrupole cavity with ε = 0.16 which is being partially pumped in the middle of the cavity, the spatial profile of the threshold lasing mode changes with the strength of the carrier diffusion in the system. When the diffusion strength is too weak to overcome the partial pumping of the cavity, the first threshold mode is found to have strong angular dependence in its far-field intensity output and heavily overlap the center of the cavity where the gain medium is inverted. As the diffusion coefficient is increased, the threshold lasing mode changes to become one that lives closer to the edge of the cavity, increasing its lifetime. Finally, when the gain becomes nearly uniform due to diffusion, we find that the threshold lasing mode is a whispering gallery mode.

Discussion of SPA
In the derivation of the population equation given above, Eq. (28), we took as an assumption that the beating terms from the coupling of different mode amplitudes averaged to zero and thus could be neglected. In this section we will make explicit this approximation. Similar to the discussion for a two-level gain medium in section 3.1, there are two criteria that go into the SPA: the populations do not acquire beating terms in the presence of multiple lasing modes, and that relaxation oscillations are not resonantly enhanced by being driven at the beat frequency. As discussed above, for a two level medium, the former criterion requires that γ Δ, γ ⊥ , and the later requirement that ω r ∼ √ κγ < Δ. The main difficulty encountered when generalizing these equations in the presence of multiple transitions is the absence of an explicit formula for an effective γ parameter entering these inequalities. However, the example of the N-level single transition case [31], for which we have such an explicit formula, strongly suggests if all of the {γ lu } are small in the relevant sense, (where {γ lu } is the set of non-radiative decay rates from the upper lasing level to the lower lasing level of each lasing transition) then the SPA will hold. Thus, since κ ≤ Δ (they are comparable for most of the lasers studied here), as long as √ Δ > γ j lu , ∀ j, we expect that the SPA will be satisfied. This is consistent with our FDTD simulation results, where a steady-state output was observed without any assumption of stationary behavior, as noted above in the discussion of Fig. 1.
In the limit of small violations of the SPA, it is possible to correct perturbatively for the effects of the beating populations within a generalized SALT framework. This calculation is discussed in Appendix B.

Semi-SALT: a free-carrier model
A natural extension of the discussion on atomic gain media with multiple transitions is to a treatment of bulk semiconductor gain media, in which there is a continuum of available transitions for the electrons between the conduction and valence bands. In doing so, we must also consider Pauli exclusion and Fermi-Dirac statistics. The polarization-Bloch equation for semiconductor media was originally derived by Lindberg and Koch, who took into account many-body effects [41], in which ρ cv,q is the off-diagonal density matrix element between the conduction and valence bands at electron momentum q, θ q is the dipole matrix element for the transition at q, Δε q is the re-normalized energy difference between the conduction and valence states, γ q is the dephasing rate, and V s is the Coulomb interaction. Note that the inversion term is the simplification of f c ), the probability that a conduction state is filled and the relevant valence state is open minus the probability that a valence state is filled and the conduction state is open. The macroscopic polarization field in Maxwell's equations is then given by [1] To proceed we will make the free-carrier approximation, setting V s = 0, with full understanding that Coulomb repulsion is an important effect in semiconductor lasers, and thus that the results found here are just the first step towards a more complete theory. Doing so allows one to write down the free-carrier susceptibility, where the factor of two appearing in the numerator accounts for spin degeneracy. The factors of f c,q and f v,q are simply the occupation probabilities for finding an electron at momentum q in the conduction and valence bands respectively and depend both upon the applied electric potential, φ , as well as the magnitude of the electric field within the cavity. The inversion equations, taking into account Fermi-Dirac statistics, take the form: where γ ,q is the non-radiative interband relaxation rate between the conduction and valence bands at momentum state q, and f (0) c,q is the Fermi function for the conduction band at momentum q in which we have introduced the electro-chemical potential μ + eφ . In this theory the applied voltage φ , related to the injected current, will play the role of the pump. In writing such a simple inversion equation we are neglecting intraband transitions, which is why a population equation similar to Eq. (28) is not needed here.
Next, we expand E(r,t) and ρ cv,q (r,t) as a summation of distinct lasing modes with different spatial profiles in the same manner as done previously in Eqs. (16) and (17). Using the above equations, we again assume stationary inversion of all of the constituent transitions ∂ t D q (r) = 0, assuming that the modal beating terms in the product Eρ * cv,q are negligible, so that we can solve for the susceptibility and insert it into Maxwell's wave equation, resulting in the Semi-SALT equations: is the Lorentzian linewidth. In general, Semi-SALT is able to treat both direct and indirect band gap materials, so long as the energy difference between the conduction and valence bands, Δε q , is specified. To continue to analytically treat the Semi-SALT equations here we will now assume that we are modeling a direct band gap semiconductor laser where the renormalized energy gap can be written as in which m r is the reduced mass and E g is the energy gap when q = 0. If we take both γ q and θ q to be independent of q, we find that the integral defining the real part of χ is divergent, as both the numerator and the denominator ∝ q 4 , and for large q, D q → −1. This is a known problem [2], which stems from the fact that the Lorentzian line-shape of the dipole transition is too broad and that other many-body effects truncate the line-shape. However, within the assumptions already used in this treatment the integral can be regularized by incorporating the correct q dependence in to θ q [2], θ q = θ 0 1 +¯h in which and where |λ represents a lattice periodic function,p is the momentum operator, and as such θ 0 has no further dependence upon q. Note that even when many-body effects are considered, γ q is still relatively constant as a function of q.
Numerically solving the Semi-SALT equations is substantially more computationally demanding than the usual SALT [31], or C-SALT equations discussed above. The reason for this is two-fold. First, the electric susceptibility no longer depends linearly upon the pump variable. Instead the applied electric potential, φ , appears in the Fermi-Dirac functions defining the equilibrium inversion of the semiconductor, yielding a nonlinear dependence in the susceptibility on the applied potential. This has the effect of making the lasing threshold problem non-linear in the pump variable and thus increasing the computational difficulty, both to find the first lasing threshold, and to find subsequent laser thresholds for additional modes. Second, the integral over q must be performed at each spatial location, another computationally expensive task. Both of these difficulties can be surmounted; the details can be found in Appendix C.
Even under the limiting assumptions made above, Semi-SALT is able to predict two features unique to semiconductor lasers, as seen in Fig. 5. Unlike the atomic gain media discussed above, the gain curve of a semiconductor gain media is in homogeneously broadened and asymmetric due to the continuum of transitions available above the energy gap. As such, we expect that the thresholds of lasing modes with frequencies above the energy gap will experience more gain and thus have lower thresholds, an effect clearly seen in the left panel of Fig. 5. The frequencies of semiconductor laser modes are also expected to shift away from the energy gap as the pump is increased due to Pauli blocking, additional electrons excited across the band gap from the increasing the applied electric potential must find higher energy states to transition to, yielding more energetic stimulated emission transitions. This effect is clearly predicted by Semi-SALT, as shown in the right panel of Fig. 5.
In general, quantitative computational study of semiconductor lasers is very challenging and is only feasible with supercomputers and specialized codes. We hope that incorporating the steady-state ansatz and the SIA into the electromagnetic part of the calculation, as is done in Semi-SALT, could eventually lead to more efficient computational approaches to parts of this problem.

Summary
We believe that the C-SALT equations derived here are the most general and accurate form of the steady-state semiclassical lasing equations, assuming only stationary populations and neglect of atomic and multi-level coherences, as is typically valid for most lasers of interest. As we have shown, these equations include both the effects of multiple lasing transitions and of gain diffusion and are computationally tractable, even in more than one dimension. This will allow efficient simulation and modeling of a number of new laser systems with these features. The ability to model gain diffusion revealed two relevant scales for the diffusion constant. Spatial hole burning disappears and single mode lasing (gain clamping) sets in when 4k 2 a D ∼ γ + γ SE (I). This could be demonstrated because C-SALT treats the above threshold gain competition exactly. Non-uniform pumping effects disappear when 4π 2 D/L 2 ∼ γ + γ SE (I), where L is the typical scale of the non-uniformity in the pump profile. Comparisons with experiment can be used to extract information about the relevant diffusion constant(s). Finally, a free-carrier version of SALT, semi-SALT, was derived which fully incorporates Fermi-Dirac statistics, showing the effect of Pauli blocking. It is hoped that further work along these lines will incorporate many-body effects leading to advances in semiconductor laser modeling.

A. Matrix definitions
The rate matrix, R(x), is an M × M matrix where M is the number of atomic levels, and can be spatially dependent when partial pumping is applied to the cavity. In Eq. (25), R can be calculated explicitly as where in most cases the non-radiative decay rates γ nm are a property of the atomic medium, and hence are not spatially dependent. The coupling matrices Ξ j also have size M × M in Eq. (25), and are defined as where the indices u and l refer to the upper and lower lasing levels of the constituent polarization j respectively. In Eqs. (31) and (36), these matrices have size PM ×PM where P is the number of discretized spatial points in the cavity and M is the number of atomic levels. The definitions of these matrices are then adjusted to be zero everywhere except at the location x ∈ P, and have only a potentially non-zero block of size M × M at this location in the matrix.

B. Treatment of beating populations
We will now investigate the breakdown of the SPA when the beat frequencies are picked up in the atomic population dynamics. Following Ge et al. [30], we assume that there are only two lasing modes in the cavity and will only be concerned with the lowest frequency components such as ω 1 − ω 2 , neglecting the other side band terms that are generated such as 2ω 1 − ω 2 as these will oscillate much faster. We first calculate corrections to the atomic population densities defining ρ ρ ρ(x,t) = ρ ρ ρ s (x) + (ρ ρ ρ b (x,t) + c.c.) ,