Steady-state ab initio laser theory for complex gain media
aa r X i v : . [ phy s i c s . op ti c s ] J un Steady-state ab initio laser theory for complex gain media
Alexander Cerjan, Y. D. Chong, and A. Douglas Stone ∗ Department of Applied Physics, Yale University, New Haven, CT 06520, USA School of Physical and Mathematical Sciences, Nanyang Technological University, Singapore 637371, Singapore (Dated: July 10, 2018)We derive and test a generalization of Steady-State Ab Initio Laser Theory (SALT) to treatcomplex gain media. The generalized theory (C-SALT) is able to treat atomic and molecular gainmedia with diffusion and multiple lasing transitions, and semiconductor gain media in the freecarrier approximation including fully the effect of Pauli blocking. The key assumption of the theoryis stationarity of the level populations, which leads to coupled self-consistent equations for thepopulations and the lasing modes that fully include the effects of openness and non-linear spatialhole-burning. These equations can be solved efficiently for the steady-state lasing properties bya similar iteration procedure as in SALT, where a static gain medium with a single transition isassumed. 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 varyinggain diffusion constant we demonstrate a cross-over between the regime of strong spatial hole burningwith multimode lasing to a regime of negligible spatial hole burning, leading to gain-clamping, andsingle mode lasing. The effect of spatially inhomogeneous pumping combined with diffusion is alsostudied and a relevant length scale for spatial inhomogeneity to persist under these conditions isdetermined. For the semiconductor gain model, we demonstrate the frequency shift due to Pauliblocking as the pumping strength changes.
PACS numbers: 42.55.Ah, 42.55.Px, 42.70.Hj, 42.60.Pk
I. INTRODUCTION
Laser theories describe the properties of self-organizedelectromagnetic oscillators based on stimulated emissionin an inverted gain medium. In semiclassical laser the-ory the electromagnetic equations are simply Maxwell’sequations, typically in the form of a Helmholtz equation,coupled to matter through a polarization or susceptibilitywhich describes the pumped, inverted gain medium. Thesusceptibility is sometimes modeled as a classical chargedoscillator, but more typically is calculated quantum me-chanically. The use of Maxwell’s equations implies thatquantum fluctuations of the electromagnetic fields are ne-glected at this level of the theory. In the full quantumtheory of the laser the fields are represented by operatorsand calculation of their higher moments yields the fluc-tuation properties. However most of the basic propertiesof interest for a laser are determined by the semiclassicallaser equations, which in addition to Maxwell’s equationsinclude equations for the polarization and atomic popu-lations of the gain medium (the diagonal and off-diagonalmatrix elements of the density matrix). These propertiesinclude the lasing thresholds, modal intensities, modalfrequencies, and spatial variation of the internal fieldsand output intensities, as well as the spatial variation ofthe inversion of the gain medium and some other dynami-cal properties such as the relaxation oscillation frequencyand decay rate of each lasing mode. The full semiclassicallaser equations are coupled non-linear partial differentialequations in space and time and are not easily solved even ∗ [email protected] numerically. In typical textbook treatments the spatialvariation of the quantities in these equations is eithercompletely neglected or treated through an approximateexpansion in closed cavity states which are taken to be or-thogonal, hence turning the resulting equations into cou-pled ordinary non-linear differential equations in time.Furthermore, the gain medium is often described as asimple two-level atomic system and modeled using theBloch equations.The past two decades have seen an explosion in thetypes and configurations of novel laser systems. Thisdevelopment has been fueled by advances in microfabri-cation techniques and motivated by applications to in-tegrated on-chip optics, for more efficient optical com-munications and sensing, as well as basic scientific inter-est. Many of these advances involve lasers with complexgain media, such as semiconductor lasers [1, 2], quan-tum cascade lasers [3, 4], and rotationally excited gases[5], for which the two level atomic system approxima-tion is poorly suited. For example, though the bandstructure of the semiconductor gain medium can be ap-proximated as a series of two level atomic transitions,multiple transitions are required to represent the effectsof Pauli blocking [6]. Cascaded-transition quantum cas-cade lasers are designed with two lasing transitions tooperate at longer wavelengths [7, 8]. Additionally, forboth rotationally excited gases and semiconductor mediathe carriers are allowed to diffuse through the cavity, aneffect that is not addressed in the Maxwell-Bloch equa-tions. Traditional methods for simulating lasers, suchas Finite-Difference Time-Domain (FDTD) calculations,are too computationally demanding to study a large pa-rameter space of the system, particularly if it is necessaryo include carrier diffusion effects [9–11],Recently, a theoretical approach has been proposed byT¨ureci, Stone, et al. to obtain directly the steady statesolutions for complex laser systems without time inte-gration [12–15]. This theory, the Steady-state Ab initioLaser Theory (SALT), employs only a single importantapproximation, the stationary inversion approximation(SIA), to treat multimode lasing, and no significant ap-proximation in treating single-mode lasing, and it treatsthe spatial degrees of freedom exactly. The SALT equa-tions are frequency domain wave equations for the lasingmodes, coupled non-linearly through the spatially vary-ing cross gain saturation. These equations can be solvedefficiently for the steady-state properties of laser cavitiesof arbitrary complexity, in any number of dimensions, us-ing a specialized non-hermitian basis set which will be in-troduced below. This approach allows a clearer and par-tially analytic understanding of the lasing solutions, andhas already led to a number of new discoveries, such asmode frustration in partially pumped cavities [16], con-trol of emission properties of random lasers through se-lective pumping [17–21], and a more general form of theSchawlow-Townes linewidth formula [22, 23]. For singletransition gain media without diffusion, SALT has alsobeen shown to give excellent agreement with FDTD sim-ulations at a substantially reduced computational cost[24, 25]. In the case of a single transition the gain sus-ceptibility for N-level lasers can be explicitly calculatedin terms of the stationary inversion profile and the lasingfields and inserted into the wave equation so as to requireonly a single set of coupled electric field equations [25]. Itcan then be shown that the N-level model can be reducedto the two-level Maxwell-Bloch model with renormalizedrelaxation and pump parameters [25].Such a simplification is not possible with multiple las-ing transitions, but, as we will show, a stationarity as-sumption (The Stationary Population Approximation)which is only slightly more general than that used inSALT, leads instead to coupled sets of field and popu-lation equations that can be solved iteratively almost asefficiently as the SALT equations. The generalized the-ory, which we call complex-SALT (C-SALT), allows usto treat steady-state lasing with (i) an arbitrary numberof atomic levels and lasing transitions, and (ii) gain dif-fusion. As already noted, direct integration in space andtime of the lasing equations with diffusing gain centers[9–11] is very challenging. C-SALT can obtain the samesteady-state solution with relative ease, allowing one touse the method for exploration of basic laser physics, and,potentially, for device design.Below we first solve the C-SALT equations for mul-timode lasing with multiple transitions but withoutgain diffusion and demonstrate excellent agreement withFDTD simulations. Including gain diffusion in the C-SALT equations and studying both 1D and 2D laserscavities, we identify two distinct physical situations inwhich carrier diffusion can have a substantial effect. Thefirst is when the carrier diffusion competes with spatial hole-burning in a uniformly pumped cavity, leading to atransition between multimode lasing to a gain-clampingregime in which only a single lasing mode reaches thresh-old (as one finds when spatial hole-burning is absent [26]).The second situation is when carrier diffusion competeswith non-uniform pumping so as to mitigate the spatialselection effects which occur in the absence of diffusion.We are not aware of prior quantitative theoretical studiesof these effects.The remainder of this paper is organized as follows. Insection II we introduce the semiclassical lasing equationsand analyze the most general conditions under whichthey lead to steady-state multimode lasing. Section IIIderives the equations for the steady-state (A) for a sin-gle transition, leading to the SALT equations (B),(C) formultiple lasing transitions and gain diffusion, leading tothe C-SALT equations, which are the main result of thiswork. Section IV demonstrates quantitative agreementbetween SALT and FDTD simulations as well as show-ing the transition between spatial hole-burning and gainclamping as the diffusion constant is increased. In sectionV a discussion of the Stationary Population Approxima-tion is presented. Section VI demonstrates how SALTcan be used to simulate lasers described by a free-carriersemiconductor model. Finally, a conclusion is given insection VII. II. OVERVIEW OF THE SEMICLASSICALLASING EQUATIONS
As noted above, the semiclassical lasing equations ne-glect the operator nature of the electromagnetic fields, sothat the lasing fields are determined by Maxwell’s waveequation (coupled to a quantum gain medium), h ∇ × ∇ × − ε c c ( x ) ∂ t i E ( x , t ) = 4 πc ∂ t P g ( x , t ) , (1)where the effect of the passive cavity is described bythe linear cavity dielectric function, ε c ( x ), which in gen-eral varies in space and frequency (although typically wewill neglect the frequency variation). The amplifying re-sponse of the gain medium is described by the non-lineargain polarization, P g , which acts as a source for thisequation and includes contributions from all of the las-ing transitions of every gain atom in the cavity. At thisstage we will assume that the gain medium arises from in-dependent identical atomic, molecular or defect centers;the case of semiconductor gain media with band excita-tions acting as gain centers will described in Section VIbelow.The gain polarization can be expressed as P g ( x , t ) = − X α δ ( x − x ( α ) )Tr[ˆ ρ ( α ) e ˆ x α ] , (2)in which the index α runs over all of the gain atoms inthe cavity at locations x ( α ) , e is the charge of an electron,nd the M × M density matrix of atom α is denoted byˆ ρ ( α ) , where M is the number of atomic levels involved inthe lasing process. Among the M levels a subset of themwill contribute to lasing and support lasing transitions;the others will simply be part of the downward cascade ofelectronic excitations involved in the pumping and emis-sion in steady-state. Lasing transitions will arise frompairs of level which have sufficiently large polarizationdue to their population inversion to contribute substan-tially to a lasing line at a nearby frequency. These pairshave a dipole moment, θ ( α ) nm = e h n | ˆ x ( α ) | m i . (3)The dipole moment is always zero for n = m due tospatial symmetry. By assuming that the full polariza-tion can be expressed in term of the density matrix forindividual atoms, we are ignoring interatomic coherenceeffects which are quite small for conventional lasers (butnot e.g. for polariton lasers).To complete the semiclassical lasing equations, onemust consider the quantum equations of motion for theaverage polarization of the gain medium, ∂ t P g ( x , t ) = − N ( x ) M X n M X m ∂ t ( ρ nm ) θ mn , (4)where we have now assumed identical “atoms” and writ-ten the density matrix elements as ρ nm = h n | ˆ ρ | m i , andwe initially assume a fixed density of gain atoms N ( x ).The evolution of the density matrix can be found fromthe Heisenberg equation of motion, ∂ t ρ nm = − i ~ h n | [ H + H I , ˆ ρ ] | m i , (5)where the atomic Hamiltonian H | n i = E n | n i , andthe interaction Hamiltonian can be written as H I = e ˆ x · E ( x , t ). Upon evaluating the commutator and sim-plifying, the evolution of the density matrix elements canbe re-written as ∂ t ρ nm = − iω nm ρ nm − i ~ M X k ( θ nk ρ km − ρ nk θ km ) · E ( x , t ) , (6)where ω nm = (1 / ~ )( E n − E m ) is the transition frequency.From this equation, we can see that off-diagonal den-sity matrix elements, which determine the gain polar-ization, can couple to one another within manifolds ofatomic transitions. If we include all terms in the equa-tion of motion, then the time evolution of a specific off-diagonal element, ρ nm , will depend not only on the levelpopulations ( ρ nn , ρ mm ), but also on other off-diagonalelements, i.e. on the polarization of other transitions. Ifthis is the case one cannot arrive at lasing equations ofthe standard form and one cannot define the polarizationassociated with a specific transition. However, physicallythese off-diagonal terms correspond to coherent multi-ple excitation, leading to effects such as electrically in-duced transparency, and inversion-less lasing. Conven-tional lasers do not typically operate in this regime. To reach this regime, the non-radiative relaxation rates be-tween non-lasing transitions must be of similar order tothose between lasing levels, which makes it difficult forthe gain medium to build up the necessary inversion tolase [27]. As such, we will assume that we are in theweakly coupled polarization regime and thus that off-diagonal density matrix elements depend only on thelevel populations of that specific pair of levels. The equa-tion of motion for the off-diagonal elements for such a pairis ∂ t ρ nm = − ( γ ⊥ ,nm + iω nm ) ρ nm + i ~ ( ρ nn − ρ mm ) θ nm · E ( x , t ) , ( m = n ) (7)where we have now added the effect of environmentaldephasing on the gain atoms in the standard manner interms of a transverse relaxation/dephasing rate γ ⊥ ,nm .With the above assumption, the total gain polarizationcan now be broken up into N T constituent polarizationsof each lasing transition, P g ( x , t ) = N T X j p j ( x , t ) (8) p + j ( x , t ) = − N ( x ) ρ nm θ mn (9)where each transition, previously labeled by the pair oflevels, 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 > ∂ t p + j ( x , t ) = − ( γ ⊥ ,j + iω a,j ) p + j − id j ~ ( θ j · E ) θ ∗ j , (10)in which the properties of the j th constituent polarizationare given in terms of the dephasing rate, γ ⊥ ,j , the atomictransition frequency, ω a,j , the constituent inversion, d j = N ( x )( ρ ( j ) nn − ρ ( j ) mm ), and the dipole matrix element, θ j .The equation of motion for the density matrix, Eq.(6), also determines the evolution of the populations ineach atomic level, given by the diagonal elements of thedensity matrix, ∂ t ρ nn = − i ~ M X k ( θ nk ρ kn − ρ nk θ kn ) · E ( x , t ) . (11)As can be seen here, atomic level populations couple onlyto the constituent polarizations with which they share aransition, and atomic levels which are not a part of anylasing transition do not appear at this stage on the righthand side of the equation for the populations. Howevertypically all levels are coupled non-radiatively throughother degrees of freedom (“the bath”) and these effectsneed to be included phenomenologically in the standardmanner, leading to ∂ t ρ n = M X m = n γ nm ρ m − M X m = n γ mn ρ n − ~ N T X j ξ n,j (cid:0) p − j − p + j (cid:1) · E . (12)Here γ nm represents either the non-radiative decay ratebetween a higher level m and a lower level n or pumpingrate from lower level m to higher level n , ρ n ≡ N ( x ) ρ nn is the total number of electrons in level n of atoms atlocation x . ξ n,j represents the relationship between thepopulation of level, n and a given lasing transition, j ; ξ n,j = 1 if n is the upper level of the transition, ξ n,j = − 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 areEqs. (1), (8), (10), and (12), which define the wave equa-tion for the electric field, the total polarization in termsof the constituent polarizations, the equation of motionfor each constituent polarization, and the equation of mo-tion for the populations in each atomic level respectively.Since now the time evolution of the off-diagonal elementsof the density matrix are contained in the polarizationequations, we will henceforth represent the populations(diagonal elements) in terms of a density vector with com-ponents ρ n ( x ). III. SOLVING THE SEMICLASSICAL LASINGEQUATIONSA. Two-level gain medium
In the limit of a two-level gain medium with only asingle lasing transition, Eqs. (8), (10), and (12) are sim-plified to ∂ t P + g ( x , t ) = − ( γ ⊥ + iω a ) P + g − iN ( x ) d ~ ( θ · E ) θ ∗ (13) ∂ t d ( x , t ) = − γ k ( d − d ( x )) − ~ (cid:0) P − g − P + g (cid:1) · E (14)where γ k = ( γ + γ ) is the relaxation rate of the inver-sion and d ( x ) = (cid:18) γ − γ γ + γ (cid:19) N ( x ) (15)is the inversion in the absence of an electric field. Theseequations, (13), and (14), along with the wave equation, (1), comprise the Maxwell-Bloch equations, the most ba-sic version of semiclassical laser theory. As noted above,it is these equations which are solved directly for thesteady-state properties using SALT. For single mode las-ing the SALT solution is essentially exact; for multimodelasing it requires the stationary inversion approximation(SIA), ∂ t d = 0. This approximation holds if the multi-mode beating terms which drive time-dependence of thepopulations are negligible.The stationary inversion approximation (SIA), re-quires two conditions to be valid. First, that the re-laxation rate of the inversion, γ k , be small compared tothe modal spacing, ∆, and the dephasing rate of the po-larization, γ k ≪ ∆ , γ ⊥ . This condition is usually metfor microlasers, for example a Fabry-Perot laser cavitytypically requires L ≤ µ m for ∆ ≫ γ k [28]. Thesecond condition is that the mode spacing must be wellseparated from the relaxation oscillation frequency, so asto avoid resonantly driving relaxation fluctuations whichcould destabilize the multimode solution. The relaxationoscillation frequency is ω r ∼ √ κγ k [29], where κ is thefield decay rate of the cavity. For microcavities, κ ≤ ∆,so ω r ≤ p ∆ γ k < ∆. As such, the SIA will hold andthe multimode solution found by SALT will be correctwhen γ k ≪ ∆ , γ ⊥ [30]. Essentially the same conditionswill be required for C-SALT to describe multimode las-ing, although, as for SALT, C-SALT solutions should begenerality accurate for the single-mode case.Continuing our brief review of SALT, we now simplifyour analysis by treating slab or two-dimensional geome-tries for which the electric fields in the transverse mag-netic (TM) modes, can be taken to be a scalar, E → E ),noting that the treatment discussed here is still com-pletely applicable in geometries for which the fields mustbe treated as vectors [31]. Having done so, we make amultimode ansatz, stating that the electric field and po-larization can be broken up into N L components withdistinct frequencies representing each lasing mode, E + ( x, t ) = N L X µ Ψ µ ( x ) e − iω µ t (16) P + g ( x, t ) = N L X µ p µ ( x ) e − iω µ t (17)where the plus superscript denotes the positive frequencycomponent of the field, E = 2Re[ E + ], and Ψ µ ( x ) and p µ ( x ) are the spatial profiles of the electric field andcorresponding polarization of the lasing mode with fre-quency ω µ . The multimode ansatz allows us to match fre-quency components of the electric and polarization fieldsthrough (13), p µ = | θ | ~ dω µ − ω a + iγ ⊥ Ψ µ . (18)Upon inserting this definition into (14) and (1) and mak-ing the rotating wave approximation (RWA), which isell satisfied for all lasers of interest, one recovers thecoupled set of SALT equations. As noted above, thesetake the form of a set of wave equations, one for eachlasing mode of the cavity, which are coupled through thenon-linear hole-burning interaction in the gain medium[12, 13, 15],0 = (cid:2) ∇ + ( ε c ( x ) + 4 πχ g ( x, ω µ )) k µ (cid:3) Ψ µ ( x ) , (19) χ g ( x, ω ) = | θ | ~ d ( x ) ω − ω a + iγ ⊥ ×
11 + | θ | ~ γ k γ ⊥ P N L ν Γ ν | Ψ ν | , (20)where Γ ν = γ ⊥ / ( γ ⊥ + ( ω ν − ω a ) ) is the Lorentzian gaincurve evaluated at ω ν .To solve the SALT equations, (19) and (20), simul-taneously, each of the lasing modes is expanded over abasis, Ψ µ ( x ) = X n a ( µ ) n f n ( x ; ω µ ) (21)where the basis functions f n are fixed, (but possibly de-pendent upon the lasing frequency ω µ ), and the expan-sion coefficients a n are found via a non-linear solutionalgorithm. There are at least two useful choices for thebasis set. The first, and most developed to this point,is to use a basis set of purely outgoing non-hermitianstates which are termed the Constant Flux (CF) states[12, 15, 32]. More recently it has been shown that onecan simply use a position space basis [31] combined witha perfectly matched layer (PML) to implement the out-going boundary condition We will not go into the detailsof these computational approaches here as either can beused with the C-SALT equations which are our main fo-cus. Instead we refer the interested reader to the citedreferences; the computations presented here are based onthe CF approach.SALT assumes only a single lasing transition, and theabove equations are then typically presented in natu-ral, scaled units for electric field and inversion, E c =2 | θ | / ( ~ √ γ ⊥ γ k ), and d c = 4 π | θ | / ( ~ γ ⊥ ), that involve therelaxation rates associated with that transition. Withmultiple transitions with differing relaxation rates, allcontributing to the lasing lines, there will be no such nat-ural scaling for C-SALT. As noted above, previous workhas shown that the steady-state equations for N-level las-ing with a single transition can be exactly mapped ontothe SALT equations with renormalized relaxation param-eters, so that computationally the N-level case is negligi-bly different from the two-level Maxwell-Bloch case [25]. B. Multiple lasing transitions
When multiple lasing transitions are present, the mul-timode ansatz is expanded to include constituent polar- izations, p + j ( x, t ) = N L X µ p j,µ ( x ) e − iω µ t (22)which still allows one to match frequency components foreach constituent polarization using Eq. (10) to find p j,µ = | θ j | ~ d j ω µ − ω a,j + iγ ⊥ ,j Ψ µ . (23)To derive the C-SALT equations, this solution for theconstituent polarization is inserted into the equation ofmotion for the atomic levels, Eq. (12), and the RWA isagain made, resulting in ∂ t ρ n = M X m = n γ nm ρ m − M X m = n γ mn ρ n − N T X j | θ j | ξ n,j d j ~ γ ⊥ ,j X ν,µ Γ ν,j Ψ ∗ ν Ψ µ e − i ( ω µ − ω ν ) t ! . (24)To proceed, we first rewrite this equation in terms of a M × M population density vector at each position in thecavity, ρ ( x ), whose components are the atomic level pop-ulations ρ n ( x ). Additionally, we make a slightly strongerversion of the SIA, the stationary population approxi-mation (SPA), which states that the beating betweendifferent lasing modes does not lead to significant time-dependence in the level populations, so ∂ t ρ ≈
0. While amore detailed discussion of SPA will be given in sectionV, it is worth noting here that it is valid for most lasersystems of interest, i.e. those with fast decays into andout of the lasing levels, while the upper lasing level ofeach transition is metastable and thus long lived. Thesefast relaxation rates need not be small compared to ∆and γ ⊥ , only the relaxation rates of the metastable up-per lasing transitions need be slow in this sense, a condi-tion which is is numerically tested and confirmed by theFDTD simulations below.Thus, the above equation can be rewritten as = R ρ + N T X j | θ j | ~ γ ⊥ ,j X ν Γ ν,j | Ψ ν | ! Ξ j ρ , (25)where R is a matrix containing information about thepump and decay rates, and Ξ j is a matrix containinginformation about which level populations are coupledto the partial polarizations and constitute the j th lasingtransition. The full forms of these matrices are given inAppendix A.Eq. (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 solvedand hence will need to be solved simultaneously with theelectric field equations. In addition, in the absence ofain diffusion, the total number of gain atoms at eachpoint in space, N ( x ), is externally fixed, and is a givenof the problem. Hence the homogeneous eq. (25) is to besolved subject to the normalization condition M X n ρ n ( x ) = N ( x ) , (26)which uniquely determines the level population vector.It is convenient to incorporate this normalization con-dition directly into Eq. (25) by defining a matrix B andan M-component total number vector N ( x ) such that B ρ = N ( x ) , (27)where neither the matrix B , nor the vector N ( x ) areuniquely defined, but must be chosen to represent Eq.(26). The normalization can then be inserted into Eq.(25), resulting in ρ = R + B + N T X j | θ j | ~ γ ⊥ ,j N L X ν Γ ν,j | Ψ ν | ! Ξ j − N ( x ) . (28)To incorporate multiple lasing transitions into theSALT formalism to recover the C-SALT equations, allthat must be altered is the equation for the electric sus-ceptibility (20), not the wave equation itself. Using (23),the electric susceptibility can be written as χ g ( x, ω ) = 1 ~ N T X j | θ j | d j ω − ω a,j + iγ ⊥ ,j , (29)where d j ≡ ρ ( j ) n − ρ ( j ) m is determined from Eq. (28). Themain difference here is now the atomic population den-sities cannot be directly inserted into the wave equationthrough the use of a scalar inversion equation. Insteadwe must simultaneously solve the wave equation (19),and the population equation (28). Using a non-linear it-eration algorithm to solve the problem numerically, oneinserts an initial guess for the field profiles, { Ψ µ ( x ) } , anduses these to solve for the full spatial profile of the atomicpopulation densities. This result for the population den-sities generates a guess for the susceptibility, which canbe inserted into the wave equation and iterated back andforth to self-consistency.There is one other important difference between (20)and (29) in the case of a “partially pumped” laser, inwhich pumping is applied to only a subset of the regionscontaining gain media [17–21]. One needs to distinguishphysically 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 positivegain under conditions of uniform spatial pumping. Thisis taken into account from the specification of N ( x ) un-der the assumed conditions of spatially uniform pump-ing. 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 inour formalism as a variation with spatial position of thepumping rate parameters in the R matrix of (28). In thiscase non-pumped regions would act as absorbers for laserlight, which would automatically be taken into accountin terms of a spatial variation in the susceptibility func-tion, χ g ( x ), which would now have an absorbing form inthose unpumped regions. C. Gain diffusion
The formalism developed in Section III B can be ex-tended to include gain diffusion, a phenomenon found inmany types of gain media. A term representing this effectcan be added to Eq. (24), resulting in ∂ t ρ n = M X m = n γ nm ρ m − M X m = n γ mn ρ n + D n ∇ ρ n − N T X j | θ j | ξ n,j d j ~ γ ⊥ ,j X ν,µ Γ ν,j Ψ ∗ ν Ψ µ e − i ( ω µ − ω ν ) t ! , (30)where D n is the longitudinal diffusion coefficient for theatomic level | n i . Despite the similarities of Eqs. (24) and(30), there is one important difference, namely that theatomic populations at each spatial location are now cou-pled together. Thus, when diffusion is present, the pop-ulation density vector ρ is an M × P dimensional vectorwhose components are atomic populations at each of P discretized spatial locations. The stationary populationapproximation can still be made, resulting in a gener-alized homogeneous equation for the population densityvector, = (cid:0) R + D ∇ (cid:1) ρ + N T X j | θ j | ~ γ ⊥ ,j N L X ν Γ ν,j | Ψ ν | ! Ξ j ρ , (31)where D is the matrix of longitudinal diffusion coeffi-cients at each spatial location, and the other two matrices R , and Ξ have also now been similarly expanded over theposition basis as well. To correctly normalize Eq. (31),we note that one consequence of SPA is ∂ t P n ρ n = 0,and when performing this sum on Eq. (30), one findsthat the total population density is homogeneous in thesteady state in the presence of diffusion,0 = ∇ M X n D n ρ n ! . (32)Furthermore, the walls of the cavity prevent any flux ofgain atoms across its borders, which is represented by theNeumann boundary condition ∂ x ρ n | x =0 ,L = 0 . (33)his, together with Eq. (32), yields the normalization M X n ρ n ( x ) = N. (34)This is a simple restatement of the fact that a diffusivegain medium in a passive cavity will be evenly distributedin the steady state, unlike a non-diffusive gain mediumwhere the density of gain atoms can fluctuate dependingon their distribution.Note however, that while the diffusion condition im-plies 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 modesabove threshold will lead to a non-uniform gain even inthe presence of uniform pumping and diffusion, and thiseffect is still captured by the theory. We will see belowthat these two effects compete; as the diffusion coefficientof the medium is increased, the excitation distributionwill become more uniform in steady-state, counteractingthe effects of spatial hole-burning.Again, the normalization requirement can now be ex-pressed as B ρ = N , (35)where the matrix B and vector N are likewise expandedover the spatial basis. Inserting Eq. (35) into Eq. (31)yields the generalized level population equation ρ = (cid:2) R + D ∇ + B + N T X j | θ j | ~ γ ⊥ ,j N L X ν Γ ν,j | Ψ ν | ! Ξ j − N ( x ) . (36)The solution now proceeds in the same way as in SectionIII B. Knowing the atomic level populations ρ , we canuse Eq. (29) to solve for the susceptibility, which is theninserted into the wave equation, Eq. (19). The problemthus reduces to a set of differential equations, one perlasing mode, coupled through the generalized level popu-lation equation (36). The latter is now a much bigger ma-trix equation, M P × M P , as opposed to M × M , increas-ing the computational cost, but still keeping it within amanageable range, even for two-dimensional lasers (seebelow).These two generalizations, Eqs. (28) and (36), alongwith their coupling to SALT, Eqs. (29) and (19), are themain results of this paper. They extend the capabilitiesof SALT to cover many types of gain media, includingnewly-developed ones [5, 8]. C-SALT can thus be usedas an efficient tool for the design and study of devicesin which gain diffusion is present alongside spatial hole-burning, a regime which is challenging for traditional nu-merical methods to handle [9–11].
36 38 40 420.00.10.20.30.4
FIG. 1. Plot of modal intensities as a function of pumpstrength for a cavity with n = 1 . ω a, = 40, γ ⊥ , = 4, θ = . ω a, = 38, γ ⊥ , = 3, θ = .
1, and 6atomic levels in total, with decay rates as indicated in theschematic. Results from C-SALT using the stationary popu-lation approximation are shown as straight lines, results fromFDTD simulations are shown as triangles. The different col-ors indicate different lasing modes. Inset shows the modalfrequencies and their intensities at P = . c/L . IV. RESULTS
To perform a well controlled test of the stationary pop-ulation approximation in C-SALT, for the case of mul-tiple transitions without diffusion, we studied 1D micro-cavity lasers for which FDTD simulations (without diffu-sion) are tractable. We used a FDTD scheme similarto the one proposed by Bid´egaray [33], altered to in-clude multiple lasing transitions and additional atomicpopulations. The simulations were run for a total du-ration at least 40 times greater than the longest timescale in the model, to ensure a steady state was reached.The simulated laser cavity consists of a dielectric slab ofbackground refractive index n = 1 .
5, with a perfectly-reflecting mirror on one side and an interface with airon the other. Distributed uniformly within the slab isa six-level gain medium, with two atomic transitions ofslightly different frequencies and widths. For these setof simulations, the diffusion constant was set to zero inthe C-SALT equations. The gain linewidths of the tran-sitions overlap, as shown in the inset of Fig. 1, so eachlasing mode receives significant gain contributions fromboth transitions. As shown in Fig. 1, excellent agreementis seen between C-SALT and FDTD simulations, both inpredictions of modal intensity and frequency, thus quan-titatively verifying the use of SPA. In these simulations,only the relaxation rates of the metastable upper lasinglevels are small compared to ∆ and γ ⊥ ,j ; the relaxationrates of other atomic levels are of the same order. Thisprovides evidence for our earlier claim that the only rateswhich must be small when compared to ∆ and γ ⊥ ,j arethose of the metastable upper lasing state. .1 0.2 0.3 0.4 0.5 0.6 0.70.00.51.01.52.02.5 FIG. 2. (Top panel) Plot of the modal intensities calculatedusing C-SALT as a function of pump strength for a dielectricslab cavity, n = 1 .
5, and a two level, single transition atomicgain medium, with ω a = 40 and γ ⊥ = 4, values again in unitsof c/L . Simulations for three different diffusion strengths areshown in solid, dashed, and dot-dashed lines. Different col-ors correspond to different lasing modes within each simula-tion. (Bottom panel) Plot of the inversion in the cavity asa function of position in the cavity at a pump strength of d = 0 . We now move to the case of lasers with multimodelasing and gain diffusion, for which FDTD is a verychallenging computational effort, but which are rela-tively tractable using the C-SALT approach. In Fig.2, we demonstrate how gain diffusion affects the tran-sition from gain-clamped single-mode lasing (which canbe described by the simplest form of laser theory [26])to multimode lasing (which is possible as a result of spa-tial hole-burning). The top panel of Fig. 2 shows C-SALTsimulations for a two level atomic medium with three dif-ferent values for the diffusion coefficient. The solid linesshow the modal intensities of the medium without diffu-sion, and dotted and dot-dashed lines of the same colorshow the evolution of the modal intensities as the diffu-sion coefficient is increased. We observe, as expected,that increasing the diffusion coefficient postpones thetransition from single-mode to multimode operation, byincreasing the threshold of the second and higher-orderlasing modes. In fact, for the largest diffusion coefficient the third lasing mode does not reach threshold withinthe pump values that were simulated. The bottom panelof Fig. 2 shows the inversion of the gain medium insidethe cavity as a function of position within the cavity fora given pump value, d = 0 . d , which can be done as there is only a singlelasing transition [25], ∂ t d = − γ k ( d − d ( x ))+ D ∇ d − | θ | ~ γ ⊥ N L X ν Γ ν | E ν | d, (37)where d ( x ) is the equilibrium value of the inversion den-sity in the absence of both fields inside the cavity anddiffusion. The final term on the right hand side can beidentified as the rate of stimulated emission, γ SE ( I ) = 4 | θ | ~ γ ⊥ N L X ν Γ ν | E ν | , (38)which is spatially dependent and proportional to the in-tensity of the local electric field. Using the SPA, Eq. (37)can be solved for the inversion density, as d ( x ) = (cid:20) γ SE ( I ) γ k − Dγ k ∇ (cid:21) − d ( x ) . (39)Thus, for diffusion to be germane to the system, k D Dγ k & γ SE ( I ) γ k , (40)where k D is the wavevector associated with the scale ofthe inhomogeneity in the inversion. If the inversion isless than this, the gain atoms are unable to move veryfar before they either non-radiatively decay or undergostimulated emission, preventing the diffusion from wash-ing out the effects of the spatial inhomogeneity in theinversion. For spatial hole-burning, this variation in theinversion is on the order of the atomic transition wave-length, k D = 2 k a , as the inversion oscillates twice as fastas the electric field [4]. This prediction is consistent withthe numerical results shown in Fig. 2, where the onset .2 0.3 0.4 0.50.00.20.40.60.8 FIG. 3. (Top panel) Plot of the modal intensities calculatedusing SALT as a function of pump strength for a partiallypumped dielectric slab cavity, n = 1 .
5, containing a fourlevel, single transition atomic gain medium. Simulations offour different values of diffusion are shown as solid, dotted,dot-dashed, and dashed lines. The first lasing mode to turn onin all of the simulations is shown in red, and the second lasingmode, which only turns on in two of the simulations, in blue.(Bottom panel) Plot of the inversion in the cavity as a func-tion of position in the cavity at a pump strength of d = 0 . of strong diffusion suppresses multimode operation andleads to gain clamped behavior.For partially pumped cavities [17–21], there is, in ad-dition to the scale associated with smoothing out spa-tial hole-burning, another relevant scale for measuringthe strength of diffusion. This is the scale at whichthe diffusion begins to overcome the spatially inhomo-geneous pumping, an inhomogeneity on the scale of theentire cavity length. Quantitatively, this can still be un-derstood from Eq. (40), except that now k D = 2 π/L ,where L is the length of the cavity, and we are assum-ing that the variation of the pumping is on the scale ofthe entire cavity. This leads to the criterion that when4 π D/L > γ k + γ SE ( I ) diffusion will strongly reduce theeffects of non-uniform pumping by allowing the inversionto penetrate into the non-pumped regions.Both types of diffusion-induced transitions are demon-strated in Fig. 3, which shows the results of C-SALT ( (cid:1) -4 ) FIG. 4. Plot of the non-interacting modal thresholds as afunction of the diffusion coefficient in a quadrupole cavitywith ǫ = 0 . r = 3 . µm , k = 6 . µm − , γ ⊥ = . µm − , n = 3 + . i . Only the three modes that become the thresh-old lasing mode for different values of the diffusion are plot-ted. Left schematic within the plot shows the boundary ofthe simulated region, black circle, boundary of the cavity,blue quadrupole, and applied pump profile, red circle. Theright schematic within the plot depicts the four level gainmedium that fills the entire quadrupole cavity. Markers in theplot correspond to the threshold lasing mode and correspond-ing inversion profile accounting for the effects of diffusion attheir locations. The inversion profile is shown in a false colorplot, with red corresponding to large inversion, and blue cor-responding to small inversion. The white boundary in theplots denotes the cavity boundary. imulations of a four-level gain medium with a singletransition, which is only pumped in half the cavity. Atwo-level medium would not be suitable for this test asit would have strong absorption in the unpumped re-gion. As the diffusion coefficient is increased, the firsttransition from the spatial hole-burning regime to thespatially-averaged gain saturation regime is observed, aswell as the effects of gain clamping which increases thesecond lasing threshold. In this regime, the inversiondoes not penetrate into the unpumped region. However,as the diffusion coefficient is further increased, the in-verted atoms are able to penetrate further into the un-pumped region. For this example, the effective spatialscale of the partial pumping is the length of the cavity, L ,by choice. As expected, we observe this second transitionwhen 4 π D/L γ k ∼ γ SE ( I ) /γ k . Once the first tran-sition has occurred the gain is sufficiently clamped, eventhough it is not uniform over the entire cavity, and weonly find single mode lasing. If one intentionally pumpsnon-uniformly on a smaller scale than L , that would de-crease the diffusion coefficient needed for the second tran-sition until in the case of wavelength scale non-uniformpumping the two transitions would coincide roughly.To demonstrate the scalability of C-SALT to multipledimensions we ran 2D TM simulations of quadrupole-shaped dielectric cavities, whose boundary is defined bythe equation r ( θ ) = r (1 + ǫ cos(2 θ )), where ǫ representsthe degree of deformation from a circular cavity. Suchcavities have been extensively studied in the context ofwave chaos theory and experiments [32, 34, 35]. It wasfound that the spatial profile of the first lasing modedepends upon both the deformation and the pumpingprofile [32, 34, 35]. Here we show in Fig. 4 that for aquadrupole cavity with ǫ = 0 .
16 which is being partiallypumped in the middle of the cavity, the spatial profile ofthe threshold lasing mode changes with the strength ofthe carrier diffusion in the system. When the diffusionstrength is too weak to overcome the partial pumpingof the cavity, the first threshold mode is found to havestrong angular dependence in its far-field intensity out-put and heavily overlap the center of the cavity where thegain medium is inverted. As the diffusion coefficient isincreased, the threshold lasing mode changes to becomeone that lives closer to the edge of the cavity, increasingits lifetime. Finally, when the gain becomes nearly uni-form due to diffusion, we find that the threshold lasingmode is a whispering gallery mode.
V. DISCUSSION OF SPA
In the derivation of the population equation givenabove, (28), we took as an assumption that the beatingterms from the coupling of different mode amplitudes av-eraged to zero and thus could be neglected. In this sec-tion we will make explicit this approximation. Similarto the discussion for a two-level gain medium in sectionIII A, there are two criteria that go into the SPA: the populations do not acquire beating terms in the presenceof multiple lasing modes, and that relaxation oscillationsare not resonantly enhanced by being driven at the beatfrequency. As discussed above, for a two level medium,the former criterion requires that γ k ≪ ∆ , γ ⊥ , and thelater requirement that ω r ∼ √ κγ k < ∆. The main dif-ficulty encountered when generalizing these equations inthe presence of multiple transitions is the absence of anexplicit formula for an effective γ k parameter enteringthese inequalities. However, the example of the N-levelsingle transition case [25], for which we have such an ex-plicit formula, strongly suggests if all of the { γ lu } aresmall in the relevant sense, (where { γ lu } is the set ofnon-radiative decay rates from the upper lasing level tothe lower lasing level of each lasing transition) then theSPA will hold. Thus, since κ ≤ ∆ (they are compa-rable for most of the lasers studied here), as long as √ ∆ > q γ jlu , ∀ j , we expect that the SPA will be satisfied.This is consistent with our FDTD simulation results, asnoted above.In the limit of small violations of the SPA, it is possibleto correct perturbatively for the effects of the beatingpopulations within a generalized SALT framework. Thiscalculation is discussed in Appendix B. VI. SEMI-SALT: A FREE-CARRIER MODEL
A natural extension of the discussion on atomic gainmedia with multiple transitions is to a treatment of bulksemiconductor gain media, in which there is a contin-uum of available transitions for the electrons betweenthe conduction and valence bands. In doing so, we mustalso consider Pauli exclusion and Fermi-Dirac statistics.The polarization-Bloch equation for semiconductor me-dia was originally derived by Lindberg and Koch, whotook into account many-body effects [36],( ~ ω − ∆ ε q + i ~ γ q ) ρ cv,q ( r, ω )= ( f c,q ( r ) − f v,q ( r )) h g q E ( r, ω )+ 1 V X q ′ V s ( q − q ′ ) ρ cv,q ′ ( r, ω ) i , (41)in which ρ cv,q is the off-diagonal density matrix elementbetween the conduction and valence bands at electronmomentum q , g q is the dipole matrix element for thetransition at q , ∆ ε q is the re-normalized energy dif-ference between the conduction and valence states, γ q is the dephasing rate, and V s is the Coulomb interac-tion. Note that the inversion term is the simplification of f c (1 − f v ) − f v (1 − f c ), the probability that a conductionstate is filled and the relevant valence state is open mi-nus the probability that a valence state is filled and theconduction state is open. The macroscopic polarizationeld in Maxwell’s equations is then given by [1] P ( r, ω ) = 1 V X q g q ρ cv,q ( r, ω ) . (42)To proceed we will make the free-carrier approxi-mation, setting V s = 0, with full understanding thatCoulomb repulsion is an important effect in semiconduc-tor lasers, and thus that the results found here are justthe first step towards a more complete theory. Doing soallows one to write down the free-carrier susceptibility, χ ( r, ω ) = Z d q π ) g q f c,q ( r ; φ, | E | ) − f v,q ( r ; φ, | E | ) ~ ω − ∆ ε q + i ~ γ q , (43)where the factor of two appearing in the numerator ac-counts for spin degeneracy. The factors of f c,q and f v,q are simply the occupation probabilities for finding anelectron at momentum q in the conduction and valencebands respectively and depend both upon the appliedelectric potential, φ , as well as the magnitude of the elec-tric field within the cavity. The inversion equations, tak-ing into account Fermi-Dirac statistics, take the form: ∂ t D q ( r ) = − γ k ,q ( D q ( r ) − D (0) q ) − i ~ (( Eg ∗ q ρ cv,q ) ∗ − c.c. ) (44) D q ( r ) = f c,q ( r ) − f v,q ( r ) (45) D (0) q = f (0) c,q − f (0) v,q = 1 e β ( ε c,q − µ − eφ ) + 1 − e β ( ε v,q − µ + eφ ) + 1 , (46)where γ k ,q is the non-radiative interband relaxation ratebetween the conduction and valence bands at momentumstate q , and f (0) c,q is the Fermi function for the conductionband at momentum q in which we have introduced theelectro-chemical potential µ + eφ . In this theory the ap-plied voltage φ , related to the injected current, will playthe role of the pump. In writing such a simple inversionequation we are neglecting intraband transitions, whichis why a population equation similar to Eq. (28) is notneeded here.Next, we expand E ( r, t ) and ρ cv,q ( r, t ) as a summationof distinct lasing modes with different spatial profiles inthe same manner as done previously (16), (17). Using theabove equations, we again assume stationary inversion ofall of the constituent transitions ∂ t D q ( r ) = 0, assumingthat the modal beating terms in the product Eρ ∗ cv,q arenegligible, so that we can solve for the susceptibility andinsert it into Maxwell’s wave equation, resulting in the Semi-SALT equations:[ ∇ + ( ε c + 4 πχ g ( r, ω )) k µ ]Ψ µ = 0 (47) χ g ( r, ω ) = Z d q π ) g q f (0) c,q − f (0) v,q ~ ω − ∆ ε q + i ~ γ q ×
11 + g q γ k ~ P ν Γ ν,q | Ψ ν ( r ) | (48)where Γ ν,q = ~ γ q ( ~ ω ν − ∆ ε q ) + ~ γ q (49)is the Lorentzian linewidth.To solve the Semi-SALT equations we will assume thatwe are modeling a direct band-gap semiconductor laserwhere the renormalized energy gap can be written as∆ ε q = ~ q m r + E g (50)in which m r is the reduced mass and E g is the energy gapwhen q = 0. If we take both γ q and g q to be independentof q , we find that the integral defining the real part of χ is divergent, as both the numerator and the denominator ∝ q , and for large q , D (0) q → −
1. This is a known prob-lem [2], which stems from the fact that the Lorentzianline-shape of the dipole transition is too broad and thatother many-body effects truncate the line-shape. How-ever, within the assumptions already used in this treat-ment the integral can be regularized by incorporating thecorrect q dependence in to g q [2], g q = g ~ q m r E g (51)in which g = h λ ′ | ˆ p | λ i (52)and where | λ i represents a lattice periodic function, ˆ p isthe momentum operator, and as such g has no furtherdependence upon q . Note that even when many-bodyeffects are considered, γ q is still relatively constant as afunction of q .Numerically solving the Semi-SALT equations is sub-stantially more computationally demanding than theusual SALT [25], or C-SALT equations discussed above.The reason for this is two-fold. First, the electric suscep-tibility no longer depends linearly upon the pump vari-able. Instead the applied electric potential, φ , appearsin the Fermi-Dirac functions defining the equilibrium in-version of the semiconductor, yielding a nonlinear de-pendence in the susceptibility on the applied potential.This has the effect of making the lasing threshold prob-lem non-linear in the pump variable and thus increasingthe computational difficulty, both to find the first las-ing threshold, and to find subsequent laser thresholdsor additional modes. Second, the integral over q mustbe performed at each spatial location, another computa-tionally expensive task. Both of these difficulties can besurmounted; 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 semi-conductor lasers, as seen in Fig. 5. Unlike the atomic gainmedia discussed above, the gain curve of a semiconductorgain media is in homogeneously broadened and asymmet-ric due to the continuum of transitions available abovethe energy gap. As such, we expect that the thresholds oflasing modes with frequencies above the energy gap willexperience more gain and thus have lower thresholds, aneffect clearly seen in the top panel of Fig. 5. The frequen-cies of semiconductor laser modes are also expected toshift away from the energy gap as the pump is increaseddue to Pauli blocking, additional electrons excited acrossthe band gap from the increasing the applied electric po-tential must find higher energy states to transition to,yielding more energetic stimulated emission transitions.This effect is clearly predicted by Semi-SALT, as shownin the bottom panel of Fig. 5.In general, quantitative computational study of semi-conductor lasers is very challenging and is only feasiblewith supercomputers and specialized codes. We hopethat incorporating the steady-state ansatz and the SIAinto the electromagnetic part of the calculation, as is donein Semi-SALT, could eventually lead to more efficientcomputational approaches to parts of this problem. VII. SUMMARY
We believe that the C-SALT equations derived here arethe most general and accurate form of the steady-statesemiclassical lasing equations, assuming only stationarypopulations and neglect of atomic and multi-level coher-ences, as is typically valid for most lasers of interest. Aswe have shown, these equations include both the effectsof multiple lasing transitions and of gain diffusion andare computationally tractable, even in more than one di-mension. This will allow efficient simulation and model-ing of a number of new laser systems with these features.The ability to model gain diffusion revealed two relevantscales for the diffusion constant. Spatial hole burningdisappears and single mode lasing (gain clamping) setsin when 4 k a D ∼ γ k + γ SE ( I ). This could be demon-strated because C-SALT treats the above threshold gaincompetition exactly. Non-uniform pumping effects disap-pear when 4 π D/L ∼ γ k + γ SE ( I ), where L is the typicalscale of the non-uniformity in the pump profile. Compar-isons with experiment can be used to extract informationabout the relevant diffusion constant(s). Finally, a free-carrier version of SALT, semi-SALT, was derived whichfully incorporates Fermi-Dirac statistics, showing the ef-fect of Pauli blocking. It is hoped that further work alongthese lines will incorporate many-body effects leading toadvances in semiconductor laser modeling.
35 40 4524252627 n=1.5
FIG. 5. (Top panel) Plot of the lasing thresholds for eachof the non-interacting modes in the cavity as a function ofthe applied electric potential, φ , for a single sided, slab semi-conductor laser. The energy gap at q = 0 has been set to E g = 40, and the chemical potential set at half that, µ = 20.The non-interacting thresholds of the two modes studied inthe bottom panel are given the same colors as appear in thebottom panel. (Bottom panel) Plot of the frequencies of thetwo lasing modes as a function of the applied electric poten-tial. The frequency for the second lasing mode (red) is onlyshown after it reaches threshold. Appendix A: Matrix definitions
The rate matrix, R ( x ), is an M × M matrix where M is the number of atomic levels, and can be spatiallydependent when partial pumping is applied to the cavity.In Eq. 25, R can be calculated explicitly as R nn ( x ) = − X m γ mn ( x ) (A1) R nm ( x ) = γ nm ( x ) , ( n = m ) (A2)where in most cases the non-radiative decay rates γ nm are a property of the atomic medium, and hence are notspatially dependent. The coupling matrices Ξ j also havesize M × M in Eq. 25, and are defined asΞ ( j ) uu = Ξ ( j ) ll = − ( j ) ul = Ξ ( j ) lu = 1 (A4)where the indices u and l refer to the upper and lowerlasing levels of the constituent polarization j respectively.n Eqs. (31) and (36), these matrices have size P M × P M where P is the number of discretized spatial pointsin the cavity and M is the number of atomic levels. Thedefinitions of these matrices are then adjusted to be zeroeverywhere except at the location x ∈ P , and have only apotentially non-zero block of size M × M at this locationin the matrix. Appendix B: Treatment of beating populations
We will now investigate the breakdown of the SPAwhen the beat frequencies are picked up in the atomicpopulation dynamics. Following Ge et al. [24], we assumethat there are only two lasing modes in the cavity andwill only be concerned with the lowest frequency com-ponents such as ω − ω , neglecting the other side bandterms that are generated such as 2 ω − ω as these willoscillate much faster. We first calculate corrections tothe atomic population densities defining ρ ( x, t ) = ρ s ( x ) + ( ρ b ( x, t ) + c.c. ) , (B1)where the subscript b denotes the portion of the popula-tion density beating at the difference frequency, and s de-notes the steady state component. Using this to rewrite(30) and separating out the time dependent portions wefind = D ∇ ρ s ( x ) + R ρ s ( x )+ X j =1 i ~ (cid:0) Ψ p ∗ ,j + Ψ p ∗ ,j − c.c. (cid:1) ξ j , (B2) ∂ t ρ b ( x, t ) = D ∇ ρ b ( x, t ) + R ρ b ( x, t )+ X j =1 i ~ (cid:0) Ψ p ∗ ,j − Ψ ∗ p ,j (cid:1) e − i ( ω − ω ) t ξ j , (B3)where ξ j is the vector form of the elements ξ n,j .As we are neglecting higher order beating frequencies,there two contributions to each atomic polarization forthe first mode, p ,j , the first from the electric field oscil-lating at ω coupled to the stationary population termsand a second from the electric field oscillating at ω cou-pled to the beating polarization at frequency ω − ω . p ,j = g j γ ,j ~ (Ψ d j,s + Ψ d j, ∆ ) , (B4)where d j,s and d j, ∆ are the stationary and beating inver-sions associated with lasing transition j and γ µ,j = 1 ω µ − ω a,j + iγ ⊥ ,j , (B5)in which the time dependence has been separated fromthe beating population densities, ρ b = ρ ∆ e − i ( ω − ω ) t . Toleading order the polarization of the j th transition canbe written as p (0) µ,j = g j γ µ,j ~ Ψ µ d j,s , (B6) and this is inserted into (B3) to obtain ρ ∆ = 1 − i ( ω − ω ) (cid:0) D ∇ + R (cid:1) ρ ∆ + X j g j ~ (cid:18) γ ∗ ,j − γ ,j ω − ω (cid:19) Ψ Ψ ∗ Ξ j ρ s , (B7)which can be inverted to find ρ ∆ , ρ ∆ = M ( ω , ω )Ψ Ψ ∗ ρ s , (B8)where the matrix M ( ω , ω ) = M ( ω , ω ) ∗ contains allthe information about the diffusion rates, decay rates,and field frequencies. Using this, we are able to writedown the first correction to the atomic polarizations, p (1)1 ,j = g j γ µ,j ~ Ψ (cid:0) M ( ω , ω ) | Ψ | (cid:1) d j,s . (B9)Finally, this can be inserted into (B2) to calculate the cor-rections to the steady state population densities. How-ever, it is difficult to glean any analytic insight fromthese equations as there is no obvious choice of param-eters to compare, with multiple atomic polarization re-laxation rates and interlevel decay rates. Fortunately,as mentioned above, the beating population densities arenegligible in the parameter regimes studied here, namelywhen the upper lasing levels of the atomic transitionsare metastable and have decay rates much less than theatomic polarization decay rates and the free spectralrange. Appendix C: Computation of Semi-SALT
There are two main new difficulties that one encoun-ters upon implementing Semi-SALT that are not seen inSALT computations using atomic gain media. The firstis that the lasing threshold problem is no longer linear inthe pump parameter. To find the non-interacting laserthresholds using the TCF basis [15] one must solve theequation, η n ( ω µ ) = 4 π Z π ) d qg q D (0) q ~ ω µ − ∆ ε q + i ~ γ q ! , (C1)which cannot be reformulated as a linear eigenvalue prob-lem. Here, η n ( ω ) is the eigenvalue being tracked acrossfrequency space of the TCF basis equation for which weare attempting to determine where the threshold lasingstate it corresponds to reaches threshold. The TCF basisis defined by the equations,0 =[ ∇ + ( ε c ( x ) + η n F ( x )) ω ] u n ( x ) (C2) ∂ x u n ( L ) = iku n ( L ) , (C3)where η n is the eigenvalue for the n th solution of thisequation at frequency ω , u n ( x ) is the eigenstate of the th solution of the TCF equation at ω , and F ( x ) is thepump profile. Instead, to find all of the laser thresh-olds one must use both the real and imaginary parts ofthis equation to first solve for φ thr where Re[ η ( ω )] =4 π Re[ χ ( ω, φ thr )], and then solve for the offset in theimaginary part of this equation, δ = Im[ η ( ω )] − π Im[ χ ( ω, φ thr )] . (C4)By monitoring for when δ changes sign as the frequencyis swept through, the lasing threshold can be found.The second main computational problem arises fromattempting to directly solve the Semi-SALT equations,(47),(48), using a non-linear solver above threshold. Suchan algorithm for finding the solution to the Semi-SALTequations requires evaluating the integrals over the elec-tron momentum at every point in space for every guessat the lasing mode amplitudes and frequency until thenon-linear solver converges, and is hopelessly inefficient.Instead, for each iteration of the Semi-SALT equationsin the pump variable above threshold this difficulty canbe sidestepped by performing a Taylor expansion on theelectric susceptibility in terms of the variables solved for at each iteration above threshold, χ ( ω N +1 , ~a N +1 ) = χ ( ω N , ~a N ) + ∂χ∂~a (cid:12)(cid:12)(cid:12)(cid:12) N · ( ~a N +1 − ~a N )+ X ν ∂χ∂ω ν (cid:12)(cid:12)(cid:12)(cid:12) N ( ω ν,N +1 − ω ν,N ) (C5)where χ ( ω N , ~a N ) is the self-consistent solution for theelectric susceptibility for the applied electric potential φ N , and is a function of all of the different lasing fre-quencies present, { ω µ } and the decomposition of the allof the lasing modes { Ψ µ } , where,Ψ µ ( x ) = ~a ( µ ) · ~u ( x ; ω µ ) , (C6)is the spatial decomposition of each lasing mode over theTCF basis at the correct frequency. Using this Taylorexpansion of the susceptibility, now all of the spatiallydependent integrals can be evaluated before invoking thenon-linear solver, dramatically improving performance. ACKNOWLEDGMENTS
We thank Arthur Goestchy and Hui Cao for helpful dis-cussions. We thank the Yale High Performance Comput-ing Center for letting us use their computing resources.This work was supported by NSF grant Nos. DMR-0908437, DMR-1307632. CYD acknowledges support bythe Singapore National Research Foundation under grantNo. NRFF2012-02. [1] H. Haug and S. W. Koch, Phys. Rev. A , 1887 (1989).[2] W. W. Chow and S. W. Koch, Semiconductor-Laser Fun-damentals: Physics of the Gain Materials (Springer, NewYork, 1999).[3] J. Faist, F. Capasso, D. L. Sivco, C. Sirtori, A. L.Hutchinson, and A. Y. Cho, Science , 553 (1994).[4] A. Gordon, C. Y. Wang, L. Diehl, F. X. K¨artner,A. Belyanin, D. Bour, S. Corzine, G. H¨ofler, H. C. Liu,H. Schneider, T. Maier, M. Troccoli, J. Faist, and F. Ca-passo, Phys. Rev. A , 053804 (2008).[5] S. Chua, Y. D. Chong, A. D. Stone, M. Solja˘ci´c, andJ. Bravo-Abad, Opt. Express , 1539 (2011).[6] Y. Huang and S. Ho, Opt. Express , 3569 (2006).[7] C. Sirtori, A. Tredicucci, F. Capasso, J. Faist, D. L.Sivco, A. L. Hutchinson, and A. Y. Cho, Opt. Lett. ,463 (1998).[8] X. Huang, J. L. Zhang, V. Tokranov, S. Oktyabrsky, andC. F. Gmachl, Appl. Phys. Lett. , 051113 (2013).[9] M. Jungo, D. Erni, and W. Baechthold, Opt. QuantumElectron. , 881 (2004).[10] W. H. P. Pernice, F. P. Payne, and D. F. G. Gallagher,J. Lightw. Technol. , 2306 (2007).[11] K. B¨ohringer and O. Hess, Prog. Quantum Electron. ,159 (2008).[12] H. E. T¨ureci, A. D. Stone, and B. Collier, Phys. Rev. A , 043822 (2006).[13] H. E. T¨ureci, A. D. Stone, and L. Ge, Phys. Rev. A ,013813 (2007).[14] H. E. T¨ureci, L. Ge, S. Rotter, and A. D. Stone, Science , 643 (2008).[15] L. Ge, Y. D. Chong, and A. D. Stone, Phys. Rev. A ,063824 (2010).[16] M. Liertzer, L. Ge, A. Cerjan, A. D. Stone, H. E. T¨ureci,and S. Rotter, Phys. Rev. Lett. , 173901 (2012).[17] J. Andreasen and H. Cao, Opt. Lett. , 3586 (2009).[18] J. Andreasen, C. Vanneste, L. Ge, and H. Cao,Phys. Rev. A , 043818 (2010).[19] L. Ge, Y. D. Chong, S. Rotter, H. E. Tuereci, and A. D.Stone, Phys. Rev. A , 023820 (2011).[20] N. Bachelard, J. Andreasen, S. Gigan, and P. Sebbah,Phys. Rev. Lett. , 033903 (2012).[21] T. Hisch, M. Liertzer, D. Pogany, F. Mintert, and S. Rot-ter, Phys. Rev. Lett. , 023902 (2013).[22] Y. D. Chong and A. D. Stone, Phys. Rev. Lett. ,063902 (2012).[23] J. C. Pillay, Y. Natsume, A. D. Stone, and Y. D. Chong,Phys. Rev. A , 033840 (2014).[24] L. Ge, R. J. Tandy, A. D. Stone, and H. E. T¨ureci, Opt.Express , 16895 (2008).[25] A. Cerjan, Y. D. Chong, L. Ge, and A. D. Stone, Opt.xpress , 474 (2012).[26] H. Haken, Light: Laser Dynamics , Vol. 2 (North-HollandPhys. Publishing, New York, 1985).[27] H. Fu and H. Haken, Phys. Rev. A , 4802 (1987).[28] H. Fu and H. Haken, Phys. Rev. A , 2446 (1991).[29] F. T. Arecchi, G. L. Lippi, G. P. Puccioni, and J. R.Tredicce, Opt. Commun. , 308 (1984).[30] A. Cerjan, Y. D. Chong, and A. D. Stone, “Steady-state ab initio theory of lasers with injected signals,” Insubmission.[31] S. Esterhazy, D. Liu, M. Liertzer, J. M. Melenk, A. Cer-jan, L. Ge, A. D. Stone, S. G. Johnson, and S. Rotter, “A scalable approach for the numerical computation ofthe steady-state ab-initio laser theory,” In submission.[32] L. Ge, Steady-state Ab Initio Laser Theory and its Ap-plications in Random and Complex Media , Ph.D. thesis,Yale University (2010).[33] B. Bid´egaray, Numer. Meth. Partial Differential Equa-tions , 284 (2003).[34] J. U. N¨ockel and A. D. Stone, Nature , 45 (1997).[35] C. Gmachl, F. Capasso, E. E. Narimanov, J. U. N¨ockel,A. D. Stone, J. Faist, D. L. Sivco, and A. Y. Cho, Science , 1556 (1998).[36] M. Lindberg and S. W. Koch, Phys. Rev. B38