Table Of Content

Diagrammatic semiclassical laser theory Oleg Zaitsev1,2,∗ and Lev Deych3,† 1Physikalisches Institut der Universita¨t Bonn, Nußallee 12, 53115 Bonn, Germany 2Fachbereich Physik, Universita¨t Duisburg-Essen, Lotharstraße 1, 47048 Duisburg, Germany 3Physics Department, Queens College of City University of New York, Flushing, NY 11367, U.S.A. We derive semiclassical laser equations valid in all orders of nonlinearity. With the help of a diagrammatic representation, the perturbation series in powers of electric field can be resummed in terms of a certain class of diagrams. The resummation makes it possible to take into account a weakeffectofpopulationpulsationsinacontrolledway,whiletreatingthenonlinearityexactly. The proposedlasertheoryreproducestheall-ordernonlinearequationsintheapproximationofconstant 0 populationinversionandthethird-orderequationswithpopulation-pulsationterms,asspecialcases. 1 The theory can beapplied to arbitrarily open and irregular lasers, such as random lasers. 0 2 PACSnumbers: 42.55.Ah,42.55.Zz n a J I. INTRODUCTION introduced for outside of the resonator. The openness of 7 the cavity in this approach is reproduced by introduc- 2 ing coupling between discrete set of inside modes and Interestin randomlaserswith coherentfeedback [1, 2] the continuous spectrum of outside modes. This way, ] andlasersbasedonchaoticmicroresonatorswithoutmir- s the modal expansion of the field becomes possible both rors [3] revealed a number of shortcomings of conven- c inside and outside of the cavity. i tionallaser theory [4, 5] that complicatedits application pt tosuchsystems. Amongthepropertiesthatcharacterize In the present work we use a different approach based o randomand chaotic lasersare strong openness, irregular onkeepingthespectralparameteroftheoutsideoutgoing . spatial dependence of the refractive index, and, possi- field real, while making inside field to satisfy continuity s c bly, nontrivial shape of the resonator. Lasing in these conditions at the boundary of the cavity. This approach si systemscanbe accompaniedbystrongcouplingbetween also results in the discrete spectrum of the cavity field y the modes, which requires a more careful treatment of with complex-valued frequencies, but, since the outside h nonlinear effects than is necessary for regular lasers. field is forced to depend on real spectral parameter, it p does not diverge and is characterizedby a constant flux. An essential part of a laser description is the choice [ ThisapproachwasusedinRef.[14]forspecialkindofvi- of an appropriate basis of normal modes in which elec- 1 tromagnetic field and other system functions can be ex- brationproblems andwasadaptedfor lasersinRef. [15], v panded. In an open system it is not possible to de- whererespectivemodesweredubbedthe“constant-flux” 2 (CF) modes. One can introduce two adjoint biorthogo- fine a Hermitian eigenvalue problem whose eigenfunc- 5 nalsystemsofCFmodes,whichcanbeusedtorepresent tionswouldformanorthogonalbasis. Instead,onehasto 9 field inside the cavity. introduce a biorthogonalsystem of so called quasimodes 4 In the standard semiclassical laser theory [4, 5] las- . as left and right eigenfunctions of a non-Hermitian op- 1 ingmodesareusuallytakentocoincidewithquasimodes erator. A number of methods to construct quasimodes 0 ofthe respectivecavities,while their amplitudesandfre- were discussed in the literature. Among the earliest are 0 quenciesarefoundfromequationsbasedonperturbation 1 the Fox-Li modes [6–8] which are useful in resonators expansionscontainingtermslinearandcubicinthefield. : with a preferred propagation direction and clearly de- v In random lasers this picture needs to be revised. First, fined transverse plane. In a more general setting, there i X were attempts to use solutions of the wave equation sat- it was shown [16, 17] that normal modes in the presence of gain differ from the passive modes even in the lin- r isfyingoutgoingboundaryconditionsatinfinity(Siegert- a Gamow boundary conditions), with complex eigenfre- ear approximation if the refractive index and/or unsat- urated population inversion are nonuniform. Second, it quencies corresponding to scattering resonances [9–11]. wasrealized[18–20]thatself-andcross-saturationcoeffi- However,these modes divergeatinfinity, which makesit cientsbeforethecubictermscanhavedifferentstatistical problematic to use them as a basis [9], while some ways propertiesin differentsystems, leadingto differentmode aroundthisproblemhavebeendiscussedinRefs.[10,11]. statistics. Finally, it was pointed out [15, 21] that non- Another possibility is to use the so called system-and- linear effects can significantly contribute to modification bath type approaches [12, 13], where cavities are de- ofthe lasingmodes comparedto those of the empty cav- scribed by an orthogonal system of wavefunctions of a ity. A theory suggested in Refs. [15, 21] allowed for self- Hermitian problem with an independent sets of modes consistent calculations of not only lasing frequencies but also of the spatial distributions of the respective modes. Neglectingpulsationsofthepopulationinversionthe au- ∗E-mail: [email protected] thors of Refs. [15, 21] were able to derive equations for †E-mail: [email protected] field amplitudes and frequencies beyond the usual third- 2 order approximation. an electric field E(r,t) is governedby the wave equation Inthepresentworkwealsogeneralizetheconventional ∂2 lasertheory,but,unliketheapproachofRefs.[15,21],we ǫ(r) E+ ( E)=0, (1) do not begin by introducing a special approximation for ∂t2 ∇× ∇× populationinversion. Instead,wecarryouttheperturba- whereGaussianunitswiththevelocityoflightinvacuum tionexpansionuptotheinfiniteorderinthefieldkeeping c=1 are used. all the terms which do not have fast temporal oscilla- In the absence of gain and absorption the field will tions. This also includes a part of population-pulsation decay in time due to the openness. It is convenient to contributionswhichisconsistentwiththeslowly-varying- representthe decaying field as a superposition of certain envelope approximation. The classificationof the result- normal modes, the quasimodes, that have only outgoing ingtermsbecomespossibleduetoaspecialdiagramtech- components outside the system. These modes can be nique developedto representthe terms ofthe expansion. constructed as families of constant-flux (CF) states [15] This technique, however, differs significantly from usual ψ (r,ω) depending on a real continuous parameter ω. Feynman diagrams widely used in solid state and high- k Explicitly,theCFmodessatisfythe differentialequation energy physics, because we have to deal with terms of everincreasingdegreeofnonlinearity. Thus,diagramsin [ ψ (ω)]=ω2ψ (ω) (2) ourtechniquearenotusedtoliterallyrepresenteachterm ∇× ∇× k k of the expansion, but play a more limited role, as a tool in the exterior of the cavity with the outgoing-wave assisting in the classification of the terms. Nevertheless, this technique allows for standard division of diagrams boundary conditions at infinity. Inside the system, the into connected and disconnected, with disconnected dia- same state satisfies a different equation: grams submitting to easy resummation in terms of only connected ones. The connected diagrams can be classi- 1 ψ (ω) k =Ω2(ω)ψ (ω). (3) fied according to the order of magnitude of the popula- ǫ(r)∇×"∇× ǫ(r)# k k tion pulsations. The resulting laser equations generalize the nonlinear equations of Refs. [15, 21] in two respects. p p For each ω, the complex eigenfrequency Ω (ω) is quan- k First, equations derived in this paper are dynamic, al- tized, as the eigenfunctions are required to match lowing for studying time dependence of the amplitudes, smoothly at the interface. whiletheequationsofRefs.[15,21]canonlydescribesta- Conjugate wavefunctions φ (r,ω) obey Eq. (2) out- tionarylasingoutput. Second,ourequationsincorporate k side with the incoming-wave boundary conditions and terms responsible for the population-pulsation contribu- the equation tion in all orders of the perturbation theory; the equations of Refs. [15, 21] are reproduced in our approach if 1 φ (ω) only the lowest-order diagram is taken into account. k =[Ω2(ω)]∗φ (ω). (4) The structure of our paper is as follows. In Sec. II we ǫ(r)∇×"∇× ǫ(r)# k k recall the definition and properties of the constant-flux p p quasimodesof opensystem. Standardsemiclassicallaser inside the system. The CF functions and their conju- equationsarewritteninSec.III infrequencyrepresenta- gates are biorthogonal and can be chosen to satisfy the tionforlaterconvenience. Coupledequationsforelectric condition field, polarization, and population inversion are reduced to equations for the field alone in Sec. IV using infinite- drφ∗(r,ω) ψ (r,ω)=δ , (5) k · k′ kk′ order perturbation theory. In Sec. V we formulate the ZI diagrammatic technique and resum the perturbation ex- where the integration is over the interior . pansion in terms of connected diagrams. In Sec. VI we AFouriercomponentE (r)oftheinternIalfieldcanbe reproduce the results of linear theory, third-ordertheory ω expanded in the CF modes as with population-pulsationterms,and all-ordernonlinear theory in the constant-inversionapproximation. Finally, we write corrections to the all-order theory that are of Eω(r)=ǫ−1/2(r) ak(ω)ψk(r,ω), (6) the first order in the population pulsations. k X a (ω)= dr ǫ(r)φ∗(r,ω) E (r). (7) k k · ω ZI p When continued to the exterior, this expansion yields a II. CONSTANT-FLUX QUASIMODES wave at the frequency ω propagating in the free space away from the system. In a stationary lasing regime ω We consider an open system defined by real dielec- becomes subjected to an equation, which has a discrete tric constant ǫ(r), with ǫ = 1 outside of the system’s setofsolutionscorrespondingtothefrequenciesoflasing boundary. In the Coulomb gauge [ǫ(r)E(r,t)] = 0, modes. ∇· 3 III. SEMICLASSICAL LASER EQUATIONS where the pump ∆n0(r) is assumed constant in time. Inthe semiclassicaltheory oflasers[4,5]the fields are described classically at the level of Maxwell equations and the active medium is treated by quantum mechan- ics. To this end, the wave equation (1) is written with a sourceterm,thepolarizationP(r,t)ofthegainmedium: IV. ALL-ORDER PERTURBATION THEORY ∂2 ∂2 ǫ(r) E+ ( E)= 4π P(r,t). (8) ∂t2 ∇× ∇× − ∂t2 A. Expansions for polarization and population Inthesimplestmodel,theactivemediumisacollectionof inversion homogeneously broadened two-level atoms. Their state is fully describedbyP(r,t) andthe population-inversion density ∆n(r,t) (difference between populations of the Equations (14)-(16)canbe reduced to an equation for upperandlowerlevelsperunitvolume). Thesefunctions the electric field alone using perturbation theory in the satisfy the equations of motion [5] field amplitude. In particular, one needs to construct an ∂2 +2γ ∂ +ν2 P= 2νd2E(r,t)∆n(r,t), (9) expansionofPω(r)inthe (odd)powersofthe fieldusing ∂t2 ⊥∂t − ~ Eqs.(15)and(16). Thenthisexpansionissubstitutedin (cid:18) (cid:19) Eq.(14)producingtherequiredequationforthefield. In ∂ 2 ∂ ∂t∆n−γk[∆n0(r,t)−∆n]= ~νE(r,t)· ∂tP(r,t). ttohethcoentvheinrdtioonradlelrasinerEthe(orr)y, w[4h,i5ch]Pyωie(lrd)sistheexpsaatnudreadtiuopn (10) ω terms in the rate equations. In this paper we will carry where d is the magnitude of the atomic dipole matrix out the expansion up to an arbitrary order in the field’s element, ν is the atomic transition frequency, and γ amplitude and use diagrammaticmethod to sortout the ⊥ (γ )isthepolarization(population-inversion)decayrate. respective terms. k If the right-hand side of Eq. (10) vanishes, ∆n relaxes tothe unsaturatedpopulationinversion∆n (r,t), which WebeginbyneglectingthequadratictermsinEq.(16) 0 which gives the zero-order expression for the popula- is a measure of the pump strength. The coupled equations (8)-(10)determine, inprinciple,electric fieldinthe tion inversion: ∆n(ω0) = 2π∆n0(r)δ(ω). Substitut- system, if ∆n (r,t) is given. ing this expression in the polarization Eq. (15), we 0 It is convenient,at this stage,to rewrite the equations obtain the first-order term for polarization P(1) = ω of motion in the frequency representation. The Fourier- i(d2/~γ )D(ω)∆n (r)E , where ⊥ 0 ω − transformed variables are given in the form 1 ∞ E(r,t)= Re dωE (r)e−iωt, (11) π ω ω ν −1 Z0 D(ω) 1 i − . (17) 1 ∞ ≡ − γ P(r,t)= Re dωP (r)e−iωt, (12) (cid:20) ⊥ (cid:21) ω π Z0 1 ∞ ∆n(r,t)= dω∆n (r)e−iωt (13) 2π ω This expression is substituted back in Eq. (16), from Z−∞ whichthe second-ordercorrectionto the inversion∆n(2) ω that facilitates application of the rotating-wave approx- isdetermined. Thisiterationprocedureyieldsthepertur- imation. Additionally, we assume that E and P are ω ω bation series for polarization and population inversion, negligible outside of a small vicinity of ν. Then the frequency squared can be approximated as ω2 = (ν +ω ν)2 ν2 +2ν(ω ν) and the lasing equations becom−e effec≈tively the firs−t-order differential equations in time, P (r)= P(q)(r), (18) ω ω ǫ(r)( ν2+2νω)E + ( E )=4πν2P , qXodd ω ω ω − − ∇× ∇× ∆n (r)= ∆n(q−1)(r), (19) (14) ω ω qodd X d2 ∞ [ i(ω ν)+γ ]P = i dω′E ∆n , − − ⊥ ω − 2π~ ω′ ω−ω′ Z0 (15) in odd and even powersof the electric field, respectively. ( iω+γ )∆n =2πγ ∆n (r)δ(ω) Henceforth we restrict the calculations to the case of k ω k 0 − scalar field. The general terms of these expansions are i ∞ dω′(E∗ P E P∗ ), (16) derived by induction in Appendix A. The resulting ex- − π~ ω′−ω· ω′ − ω′+ω · ω′ Z0 pressions are 4 A(q+1)/2 P(q) =2i~γ ∆n (r)D(ω) dω dω E E∗ E E∗ E ω k πq−1 0 1··· q−1 ω1 ω2··· ωq−2 ωq−1 ω−ω1+ω2−···−ωq−2+ωq−1 Z D (ω ω )D (ω ω +ω ω ) D (ω ω + +ω ω ) k 1 2 k 1 2 3 4 k 1 2 q−2 q−1 × − − − ··· − ··· − [D(ω )+D∗(ω )][D(ω ω +ω )+D∗(ω ω +ω )] 1 2 1 2 3 2 1 4 × − − ··· [D(ω ω +ω ω +ω )+D∗(ω ω +ω ω +ω )], (20) 1 2 3 q−3 q−2 2 1 4 q−4 q−1 × − −···− − −···− A(q+1)/2 ∆n(q+1) =2 ∆n (r)D (ω) dω dω E E∗ E∗ E E∗ ω πq 0 k 1··· q ω1 ω2··· ωq−1 ωq −ω+ω1−ω2+···−ωq−1+ωq Z D (ω ω )D (ω ω +ω ω ) D (ω ω + +ω ω ) k 1 2 k 1 2 3 4 k 1 2 q−2 q−1 × − − − ··· − ··· − [D(ω )+D∗(ω )][D(ω ω +ω )+D∗(ω ω +ω )] 1 2 1 2 3 2 1 4 × − − ··· [D(ω ω +ω ω +ω )+D∗(ω ω +ω ω +ω )] 1 2 3 q−3 q−2 2 1 4 q−4 q−1 × − −···− − −···− D(ω ω +ω ω +ω )+c.c.(ω ω), (21) 1 2 3 q−1 q × − −···− →− whereqisodd. Thenotationc.c.(ω ω)standsforthe to the equation →− first part of the equation to which complex conjugation [Ω (ω ) ω ]a δ(ω ω ) accompaniedby changeofthe signofω wasapplied. We k k k k k − − also introduced the following definitions = W(q) a a∗ a∗ a k,k1,...,kq k1 k2··· kq−1 kq qXodd k1X,...,kq d2 A≡−2~2γ γ , (22) ×δ(ω−ωk1 +ωk2 −···+ωkq−1 −ωkq), (25) ⊥ k −1 where the coefficients W(q) arefunctions ofthe fre- Dk(ω)≡(cid:20)1−iγωk(cid:21) . (23) qthueenrcigiehst-ωhka,nωdk1s,id..e.c,oωnktqka.,iknC1,.lt.eh.,akerqldy,elmtaosfutnocftitohnesttehrmatsaorne differentfromδ(ω ω )ontheleft-handside. Inthetime k The aboveexpressionsdescribe nonlinear (q 3)correc- − ≥ domain, these terms would introduce oscillations with tionsto polarizationandpopulationinversion,whichare frequencies different from ω . This means that, in gen- k used below to obtain equations for the field amplitudes. eral, no stationary lasing solutions exist unless, for some reasons, the terms oscillating at the “beat” frequencies are small and can be neglected. Usually the selection of the slowly changing contributions to Eq. (25) is done by B. Equations for electric field leaving only terms with pairwise coinciding indices k , i so that the respective frequency differences cancel out Equation for the amplitudes (7) of normal modes, (more on this procedure can be found below). However, since the number of ω in this equation is odd, one will ki alwaysremainwiththeexpressionδ(ω ω ),wherek′ is [Ωk(ω)−ω]ak(ω)=2πν drǫ−1/2(r)φ∗k(r,ω)Pω(r), theindexofoneoftheremaininguncan−celek′dfrequencies. ZI (24) This problem can only be resolved by requiring that, follows from Eq. (14) and biorthogonality condition (5). in a given lasing mode, l, each contribution a of the lk The right-hand side is a sum of the infinite number of quasimode k oscillates at the same frequency, and the nonlinearcorrectionstothepolarizationandis,therefore, general solution for a takes the form of k afunctionalofall amplitudesa (ω). Bydefinition,alas- k ingmodeissuchasolutionwhich,inthelimitt ,ap- ak(ω)=π alkδ(ω ωl). (26) →∞ − proachesapurelyharmonicform, exp( iω t). Hereω l k k X ∝ − is some frequency, which is determined self-consistently Amplitudes a can be shown to obey the equation lk from Eq. (24). In the frequency domain these solutions are characterized by delta-functional frequency depen- [Ω (ω ) ω ]a = V(l)a , (27) k l − l lk kk′ lk′ dence δ(ω ω ). ∝ − k Xk′ Generallyspeaking,lasingmodesaredifferentfromthe (l) quasimodesofthepassivesystem. Indeed,anattemptto where the coefficients Vkk′ depend on allfrequencies and (l) search for the solutions of Eq. (24) in the form a (ω) = amplitudes. An explicit expression for V is given be- k kk′ a δ(ω ω ) leads, after the frequency integration (20), low. One can see from this equation that the lasing k k − 5 modes are those combinations of the quasimodes that amounts to neglect of time derivatives of the nonlinear diagonalize the matrix V(l), while lasing frequencies are corrections. The resulting equation is kk′ real eigenvalues of this matrix [16, 17]. One has to realize, though, that Eq. (27) is not an or- d dinary eigenvalue problem, since the matrix Vk(kl)′ itself (cid:26)−i[1−Ω′k(ωl)]dt +[Ωk(ωl)−ωl](cid:27)alk(t) depends on the amplitudes a . Unlike linear eigenvalue lk = V(l)(t)a (t), (28) problems, that determine frequencies for which nonzero kk′ lk′ solutions for the amplitudes can exist, solving Eq. (27) Xk′ one must be able to find the frequencies, as well as which, in the stationary limit, coincides with Eq. (27). the field amplitudes. This is possible, because the requirement that the respective eigenfrequencies must be Note that V(l) is now a slowly varying function of time kk′ real provides an additional constraint on solutions of via its dependence on the amplitudes al′′k′′(t). The cor- Eq. (27) [21, 22]. rectiondue to Ω′ is a new term, whichhas not been dis- k If matrix V(l) is calculated in the constant-inversion cussedinanyoftheprevioustreatmentsoflasingdynam- kk′ ics. While itdoesnotaffectthesteady-statesolutions,it approximation(seebelow),Eq.(27)reproducesthemain mightchange their stability, andis, therefore, important result of Ref. [23]. However, while the derivation of this forstronglyopencavities. Moredetailedstudyofitsrole equation in Ref. [23] is only valid in the strictly station- isoutsideofthescopeofthispaperandwillbepresented ary limit, the argumentspresentedhere canbe extended elsewhere. toanonstationarycase. Indeed,wecanrepeattheabove (l) arguments for a weakly nonstationary situation, requir- ThepolarizationmatrixV intheslowlyvaryingam- kk′ ing that the amplitudes a be slowlychangingfunctions plitude approximation can be presented as lk of time. Formally, we replace the mode expansion (26) wsuimthedakt(oωb)e=sharlpalylk(pωea−keωdl)a,twωhe=reωa.lkI(nωt−hiωslc)asiseawse- Vk(kl)′(t)=2πν drǫ−1(r)φ∗k(r,ωl)ψk′(r,ωl)ηl(r,t), P l ZI cantransformEq.(24)tothetime domainbyexpanding (29) Ω (ω) as Ω Ω (ω )+(ω ω )Ω′(ω ), where Ω′(ω) is where we introduced the nonlinear susceptibility thke derivatikve≈of Ωk lwith res−peclt tok thle spectral pkaram- η(r,t) P (r,t)/E (r,t) defined as the ratio of the k l l l eterω. ItwasfoundinRef.[24]that,atleastinthe case slowly v≡arying polarization amplitude P (r,t) and the l of modes of a disk resonator, this derivative is not small field E(r,t) in the mode l. The expression for the sus- l and must be taken into account. In nonlinear terms we ceptibility is found from perturbation expansion for po- simply replace a (ω ω ) πa (t)δ(ω ω ), which larization P(r,t), Eqs. (18) and (20), and is given by lk l lk l − → − ηl(r,t)=2i~γkD(ωl)∆n0(r) Aq+21 r |El2(r,t)|2|El4(r,t)|2···|Elq−1(r,t)|2 qXodd l1X,...,lq D (ω ω )D (ω ω +ω ω ) D (ω ω + +ω ω ) × k l1 − l2 k l1 − l2 l3 − l4 ··· k l1 − l2 ··· lq−2 − lq−1 [D(ω )+D∗(ω )][D(ω ω +ω )+D∗(ω ω +ω )] × l1 l2 l1 − l2 l3 l2 − l1 l4 ··· [D(ω ω +ω ω +ω )+D∗(ω ω +ω ω +ω )]. (30) × l1 − l2 l3 −···− lq−3 lq−2 l2 − l1 l4 −···− lq−4 lq−1 Here the order of nonlinearity q, introduced in Eq. (18), drops out of the equation, since the amplitude E must lq determines the number of different indices l , which take be equal to some other amplitude E . It is assumed i li values from 1 to N , where N is the number of lasing that the slowly varying field amplitudes are expressedin m m modes. The superscript “r” at the sum symbol specifies terms of quasimode components a (t) of the respective lk that the possible values of the indices are restricted by lth lasing mode using the resonance condition E (r,t)=ǫ−1/2(r) a (t)ψ (r,ω ). (32) l lk k l ω ω +ω ω +ω ω =0 (31) l1 − l2 l3 −···− lq−1 lq − l Xk whichensurescancelationoffastoscillatingterms. Inthe Substituting Eqs. (30) and (29) in Eq. (28) one obtains absenceofaccidentaldegeneracies,thisconditionimplies a closed system of dynamic equations for a valid to all lk thateachofthe indicesl ,l ,...,l mustbe equaltoone orders in the field amplitude. 1 3 q of the indices l ,l ,...,l ,l. This leads, in particular, One of the fundamental difficulties of the theory of 2 4 q−1 to the appearance of absolute squares of the field in the lasers is that the number N of lasing modes is a priori m first line of Eq. (30). Moreover, the index l effectively unknown and depends on the strength of the pumping q 6 and the spatial distribution of the electric field in the l l 1 2 cavity. In Ref. [23] this value is determined by the num- l l ber ofpossible solutionsofEq.(27)withrealfrequencies 3 4 at a given pumping strength. This approach, however, does not take into account stability of the found solu- l l tions, which can only be determined by considering the q− 2 q− 1 l l time-dependent Eq. (28). Using this equation one could q start by assuming that N is equal to the size of the m basis of quasimodes, N , andstudy their time evolution. b FIG.1: Labelingofverticesinadiagramoforderq=1,3,... ThoseE whichdonotcorrespondtoreallasingsolutions l at given pumping would decay to zero, and the number (a) oflasingmodeswouldbedeterminedaposterioriwithout the need for a prior knowledge of N . This approach is l l m 1 ~0 not free of difficulties either, because of possible multi- X 11 stable behavior and hysteresis. Analysis of these issues, however, is beyond the scope of this paper. (b) (c) l l l l 1 2 1 2 V. DIAGRAMMATIC TECHNIQUE l l l l 3 ~0 3 ~0 X X 31 32 A. Diagrammatic representation of the perturbation series FIG. 2: First-order diagram (a) and third-order diagrams (b, c). Inthissubsectionwepresentadiagrammatictechnique developedtoclassifydifferentnonlineartermsinEq.(30). It should be noted, however, that our diagrams, unlike are shown in Figs. 2 and 3. diagramsofthefieldormany-particletheory,donotpro- Eachdiagramspecifiesaparticularcontributiontothe vide one-to-one correspondence between different terms series (30). The latter will be written in the form of Eq. (30) and elements of the diagrams. The role of the diagrams here is more limited: we use them to clas- η(r,t)=2i~γ ∆n (r)D(ω ) , (33) l k 0 l l sify different pairing possibilities for the lasing mode in- X e dices l ,l ,...,l ,l in the perturbation series (30). Nev- Nq erthele1ss,2as it iqs shown below, this technique allows for Xl ≡ Aq+21 Xq0j. (34) classificationandpartialsummation of the classes of the qXodd Xj=1 terms in a manner very similar to traditional diagram- e r matic methods. Unlike pairing of vertices in traditional where each X0 represents a partial sum in ( ), qj ··· diagrammatic techniques, which reflects Vick’s theorem Eq.(30),in which pairsof indices are chosento be equal X for creation-annihilationoperatorsor Gaussianstatistics to each othereaccording to the links connecting respec- of respective random processes, the pairing procedure in tive vertices in the diagram. For example, the first three the situation under consideration hinges upon the condition expressed by Eq. (31). The resonance condition guaranteesthe absence ofthe fastoscillatingterms, and, l l 1 2 hence, the validity of the slowly changing amplitude approximation. l l 3 4 To construct a diagram X0 of order q = 1,3,..., we qj l l place q+1 vertices in two columns as shown in Fig. 1. 5 ~ ~ ~ 0 0 0 X X X The left vertices are labeleed l1,l3,...,lq and the right 51 52 53 vertices are labeled l ,l ,...,l ,l. The vertex l is dif- 2 4 q−1 ferent from the other vertices, because there is no summation over the index l in Eq. (30). After that, each vertex on the left is connected with exactly one vertex on the right. The index j = 1,...,N labels all distinct q connection possibilities in an arbitrary order. To obtain ~0 ~0 ~0 X X X alldiagramsoforderq,wecan,first,econnectthevertices 54 55 56 by (q+1)/2 horizontal links and then reshuffle the vertices,say,ontheleftwithoutcuttingthelinks. Thus,the FIG. 3: Fifth-order diagrams. The dash-dotted lines are number of possible diagrams of order q is the number of thecutsthatsplitdisconnecteddiagramsintoconnecteddia- permutations N =[(q+1)/2]!. Diagrams for q =1,3,5 grams. q e 7 diagrams correspond to the following expressions: Our diagrammatictechnique possessesthe basic prop- ertythatdisconnecteddiagramsaregivenbyproductsof X0 =1, (35) their connected parts. Thus, we can write for our exam- 11 ples X0 = E (r,t)2D (ω ω )[D(ω )+D∗(ω )], e31 | l2 | k l− l2 l l2 lX26=l X505 =X101X31, (41) e (36) X0 =X0 (X )2. (42) X0 = E (r,t)22Re[D(ω )]. (37) e56 11 11 32 | l2 | l2 Themultiplicativityisduetothefactthattheresonance Xl2 condition of the typee (31) is fulfilled for each connected e subdiagram (see also Appendix B). The multiplicativity The restriction l = l in the diagram X0 excludes the 2 6 31 property allows to express the series (34) in terms of the term with l =l =l =l, which enters X0 . In general, 1 2 3 32 connected diagrams as the terms with more than two indiceseequal belong to m the diagram in which the links connectineg these indices Nq ∞ Nq do not cross each other. Another example are the fifth- = X0 X Xl  qj  qj order terms with l2 = l3 = l4 = l5 6= l, which enter X502, qXodd Xj=1 mX=0 qXodd Xj=1 but not X501. Expression for arbitrary Xq0j is given in  NqX0   Appendix B. e = qodd j=1 qj . (43) e e 1−P qodPd Nj=q1Xqj This resummatPion formPula is the main result of our pa- B. Resummation of the diagrams per. A diagramis called connected if it cannot be cut by a horizontal line without cutting a link. For instance, the VI. LIMITING CASES AND DISCUSSION diagrams X0 , X0 , X0 , and X0 are connected, while 31 51 52 53 the diagrams X0 , X0 , X0 , and X0 are disconnected. To make the meaning of Eq. (43) more transparent 32 54 55 56 To simplifyetheenotateion, we oredered all connected dia- we apply it in several well-known special cases. We grams before tehe diseconneected diagerams for given q. We start with the linear approximation in Sec. VIA and will label connected diagrams as X0 with consider the effect of gain-induced coupling of passive qj modes [16, 17]. In Sec. VIB we reproduce semiclassical X0 =X0 , j =1,...,N , (38) equations of the standard third-order laser theory [4, 5]. qj qj q In Sec. VIC we discuss all-order nonlinear theory in the approximation of constant population inversion [15, 21– where N (< N ) is thee number of connected diagrams. q q 23] and derive, using our theory, the first diagrammatic The horizontalcuts separate disconnected diagrams into correction to it. one connectedediagram containing the vertex l (denoted byanunfilleddotinthegraphicrepresentation)andsev- eral connected subdiagrams without such vertex. The A. Linear gain-induced mode coupling lattersubdiagramswillbedenotedasX ,j =1,...,N , qj q where q+1 is the number of vertices in the subdiagram. In the linear approximation to the polarization P = Inplaceofthevertexlthesediagramshaveavertexwith l ηE [Eq. (30)] only the lowest diagram X0 contributes the index lq+1 which runs over all lasing modes, as the l l 11 other indices lj. For example, the diagram X505 consists to Xl. Equation (28), where the matrix Vk(kl)′ is calcu- latedusingEqs.(29),(33),and(35),yields the following ofX0 andX andthediagramX0 consistsofX0 and 11 31 56 11 equations for the slowly-varying amplitudes: two diagrams X , where e 11 e d δ +i Ω (ω ) δ ω a =0, (44) X11 = |El2(r,t)|22Re[D(ωl2)], (39) k′ (cid:26) kk′dt kk′ l − kk′ l (cid:27) lk′ Xl2 X (cid:2) (cid:3) where the Ω′ term is henceforth neglected. The fre- X = E (r,t)2 E (r,t)2 k 31 | l2 | | l4 | quency matrix lX2,l4 l26=l4 Ωkk′(ω)=Ωk(ω)δkk′ Vkk′(ω), (45) − D (ω ω )[D(ω )+D∗(ω )]2. (40) × k l4 − l2 l4 l2 is modified by the linear gain term, A general expression for Xqj is given in Appendix B. Vkk′(ω)= Note that X is of the order q+1 in the electric field. qj d2 ∆n (r) Connected diagrams contain (q 1)/2 nontrivial factors 2iπν D(ω) dr 0 φ∗(r,ω)ψ (r,ω), (46) Dk =1. − − ~γ⊥ ZI ǫ(r) k k 6 8 proportional to the overlap integral. Clearly, the ma- strengthbeloworatthethreshold. However,thebasisof trices V and Ω are nondiagonal if the pump or di- normalmodescanbeusedasastartingpointinnonlinear kk′ kk′ electric constant are not homogeneous. In this case the theories. biorthogonalquasimodesofthesystemψ andφ areno k k longer ψ and φ , but are determined by the right and k k left eigenvectors a(r,l) of Ω according to lk kk′ ψ (r,ω)= a(r)ψ (r,ω), (47) B. Third-order theory l lk′ k′ k′ X φl(r,ω)= a(lkl)′φk′(r,ω). (48) inTthoeofibetaldin, wane kaepepprotxhiemdatiaiognratmosXXl 0ofatnhde tXhi0rdi-nortdheer k′ 11 31 X numerator and the diagram X in the denominator of 11 The righteigenvectorsa(r) are normalmodes of Eq.(44) Eq. (43) and expand the latter: lk′ whose amplitudes a (t) obey the equation l a˙l+i Ωl(ωl) ωl al =0, (49) Xl ≈X101+X101X11+X301. (52) − (cid:2) (cid:3) where Ωl(ω) are eigenvalues of Ωkk′(ω). We recall It is convenient to write lasing equations in the basis of that electric field in the mode l has a time dependence quasimodes ψ and φ that diagonalize the linear part. k k al(t)exp( iωlt). Thus, the lasing frequency ωl is deter- Due to nonlinear effects, the lasing modes above the − mined from the requirement threshold, Re[Ω (ω )]=ω (50) l l l E (r,t)=ǫ−1/2(r) a (t)ψ (r,ω ), (53) l lk k l and the threshold condition for this mode is k X Im[Ω (ω )]=0. (51) l l are, in general, linear combinations of individual quasi- As follows from Eq. (49), the mode amplitudes diverge modes. Equation for the amplitudes a (t) follows from lk exponentially above the threshold. Hence, applicabil- Eq. (28), after taking into account the results of linear ity of the linear approximation is limited to the pump theory (Sec. VIA), and has the form: a˙ +i Ω (ω ) ω a = πν d2 2D(ω ) dr∆n0(r)φ∗(r,ω )E (r,t) lk k l − l lk −~γk (cid:18)~γ⊥(cid:19) l ZI ǫ(r) k l l (cid:2) (cid:3) Nm p E (r,t)2 2Re[D(ω )]+(1 δ )D (ω ω )[D(ω )+D∗(ω )] , k =1,...,N , (54) l′ l′ ll′ k l l′ l l′ b × | | − − l′=1 X (cid:8) (cid:9) where N is the size of the basis of quasimodes and for the intensities I and phases ϕ , b l l N is the number of lasing modes. The lasing frequen- m cies ω and the mode thresholds need to be determined l from these N N equations in the stationary regime I˙ 2Im[Ω (ω )]I b m l l l l a˙ = 0 using, ×e.g., a selfconsistent procedure described − lk 2πν d2 2 in Ref. [21]. = I Re D(ω ) B I , (55) −~γk (cid:18)~γ⊥(cid:19) l " l l′ ll′ l′{···}# X ϕ˙ +Re[Ω (ω )] ω l l l − l Insomecasesthestandardassumptionofatraditional πν d2 2 lasing theory that the lasing modes coincide with the = Im D(ω ) B I , (56) quasimodes of the cavity remain valid. In this case the −~γk (cid:18)~γ⊥(cid:19) " l l′ ll′ l′{···}# totalnumberofequations(54)isreducedtoN sincethe X b amplitudesareapproximatedasa (t)=a (t)δ . Repre- lk l lk senting a (t) =√I exp(iϕ ) and separating the real and The terms enclosed in the braces are the same as those l l l imaginarypartsinEq.(54),weobtain2N realequations in Eq. (54). We defined overlap integrals for the quasi- b 9 modes as where Bll′ =ZIdr∆[ǫn(r0)(]r2)φ∗l(r,ωl)ψl(r,ωl)|ψl′|2. (57) Υn = 2~2dγ2⊥γk l′6=l|El′(r,t)|2Dk(ωl−ωl′) X The lasing frequencies ω are determined, together with l [D(ω )+D∗(ω )], (61) the stationary intensities, from Eqs.(55) and (56) in the × l l′ stitoantison(5a5ry) arengdimfreeqI˙ule=ncy0 eaqnudatϕ˙ioln=s (05.6)Tgheenerraatleizeeqtuhae- Υd = 2~2dγ2⊥γk l′′6=l′|El′′(r,t)|2Dk(ωl′′ −ωl′) standard third-order semiclassical theory with the self- X [D(ω )+D∗(ω )]2 and cross-saturation terms [4, 5] to systems with strong l′′ l′ . (62) openness and arbitrary distribution of refractive index. × 2Re[D(ωl′)] Taking into account that the electric field in the lasing C. Constant-inversion approximation and modeshasatypicalmagnitudeof~√γ⊥γk/d,onecansee corrections that deviations from the constant-inversion approximation is of the order of D (ω ω ) [the term D(ω ) is of k l l′ l − the order of unity since frequencies of lasing modes are In many physical situations the population inversion ∆n(r,t) is approximately constant in time. Oscillations concentrated within the width of the gain profile]. The differencebetweenlasingfrequenciescanbeestimatedas ofthepopulation(populationpulsations)areresponsible ω ω γ /N . Then the condition D 1 can be for the terms proportional to Dk(ωl − ωl′), l 6= l′, in |exlp−ressle′d|≈asN⊥ γm/γ 1. Giventhatγ iksu≪suallysev- the expansion of the nonlinear susceptibility (30) and, m k ⊥ ≪ k eral orders of magnitude smaller than γ , this condition hence, in the expressions for the diagrams. Comparing ⊥ isinmostsituationsfulfilled. ItwasreportedinRef.[22], thetermsinthebracesinEq.(54),wecanconcludethat however,thatnonlinearinteractionbetweenlasingmodes the population pulsations can be neglected if D (ω k l − canpushtheir frequencies towardeachother makingthe ω ) 1, i.e., ω ω γ for all lasing frequencies l′ l l′ k ≪ | − | ≫ intermodespectralintervalmuchsmallerthanthetypical ω =ω . l l′ 6 value givenabove. Such pairs of modes can resultin sig- In the approximation of constant inversion only the nificant corrections to the constant-population approxi- diagrams X0 and X , which do not contain the D 11 11 k mation. Adding to the expansion additional connected functions, contribute to (43). Equation (28) for the l X diagrams with up to q+1 vertices, one can improve the coefficientsofexpansion(32)oftheelectricfieldbecomes constant-population approximation by constructing las- d2 ∆n (r) ing equations valid in the order (q 1)/2 in Dk. a˙ +i[Ω (ω ) ω ]a =2πν dr 0 − lk k l − l lk ~γ⊥ ZI ǫ(r) φ∗k(r,ωl)El(r,t) p . (58) VII. CONCLUSIONS × 1+ d2 E (r,t)22Re[D(ω )] 2~2γ⊥γk l′| l′ | l′ We presented a diagrammatic semiclassical laser the- In contrast to Eq. (54)P, the linear contribution here is oryvalidinallordersofelectricfield. Theoriginalpertur- not diagonalizedand is containedin the right-handside. bation series in the powers of the field can be resummed Equation (58) is valid in all orders in nonlinearity if in terms of a certain class of diagrams, the connected the population inversion is constant. With the help of diagrams. The resummation allows to construct a con- Eq. (43) it is straightforward to write out the correc- trolled expansion in the small parameter γ /γ , which k ⊥ tions due to the populationpulsations. The terms ofthe is a measure of population pulsations, while treating first order in Dk ≪1 are contained in the diagrams X301 the nonlinearity exactly. Our lasing equations gener- and X31 [Eqs. (36) and (40)], so that l can be approxi- alize the all-order nonlinear equations in the constant- X mated as inversion approximation and the third-order equations X0 +X0 with population-pulsation terms. The use of constant- l 11 31 . (59) flux quasimodes as basis functions makes it possible to X ≈ 1 X X − 11− 31 apply the theory to stronglyopen andirregularsystems, These corrections modify both the numerator and de- suchasrandomlasersandlaserswithchaoticresonators. nominator of Eq. (58), which can now be presented as d2 ∆n (r) a˙ +i[Ω (ω ) ω ]a =2πν dr 0 Acknowledgments lk k l − l lk ~γ⊥ ZI ǫ(r) φ∗(r,ω )E(r,t)(1 Υ ) Financial support was provided by the Deutsche k l l − n p . × 1+ d2 E (r,t)22Re[D(ω )](1 Υ ) Forschungsgemeinschaft via the SFB/TR12 and 2~2γ⊥γk l′| l′ | l′ − d FOR557(O.Z.)andbyPSC-CUNYviagrantsNo.62680- (60) P 00 39 and No. 61788-0039 (L.D.). 10 Appendix A: Derivation of Eqs. (20) and (21) we first exchange the labels ω1 ω2, ω3 ω4, ..., ↔ ↔ ω ω and then define the variable q−2 q−1 → For a few lowest values of q the validity of Eqs. (20) and (21) can be checked directly. To prove these re- ωq−2 =ω ω′ ω1+ω2 +ωq−3+ωq−1. (A4) − − −··· lations by induction, we assume that ∆n(q−1) is given ω Transforming the integral over ω′ as by Eq. (21), then derive P(q) (20), and, finally, ob- ω tain ∆n(q+1). ω dω′D (ω ω′)E E According to Eq. (15), k − ω′ ω−ω′−ω1+ω2−···+ωq−3+ωq−1 Z Pω(q) =i~γkAπD(ω) dω′Eω′∆n(ωq−−ω1′). (A1) =Z dωq−2Dk(ω1−ω2+···+ωq−2−ωq−1) Z E E (A5) Substituting ∆n(q−1) from Eq. (21), we immediately see × ω−ω1+ω2−···−ωq−2+ωq−1 ωq−2 ω that the factor before the integral in Eq. (20) is reproduced. To calculate the integral, we consider sep- we obtain the part of the Pω(q) integrand proportionalto arately the∗ two contributions: ∆n(ωq−−ω1′) = ∆n(ωq−−ω1′) + D∗N(ωex2t−wωe1d+er·iv·e·+∆ωnq(q−+11)).using Eq. (16), (q−1) (q−1) ω ∆n ,where∆n isgivenexplicitlybyEq.(21) ω′−ω ω f a(cid:20)nfd ∆n(cid:21)(q−1) ∗ isfc.c.(ω ω). When integrating ∆nω(q+1) =π−~γi Dk(ω) dω′Eω∗′−ωPω(q′) −ω → − k Z (q−(cid:20)1) (cid:21) +c.c.(ω ω). (A6) ∆nω−ωf′ we introduce a new variable Clearly, the factor before the→in−tegral in Eq. (21) follows ω =ω′ ω+ω ω + ω +ω . (A2) aftersubstitutionofP(q) (20). Introducingthenewvari- f q−1 − 1− 2 ···− q−3 q−2 ω′ able Then the integral over ω′ becomes ω =ω′ ω +ω ω +ω (A7) q 1 2 q−2 q−1 − −···− dω′D (ω ω′)E E∗ k − ω′ ω′−ω+ω1−ω2+···−ωq−3+ωq−2 we rewrite the ω′ integral Z = dω D (ω ω + +ω ω ) q−1 k 1 2 q−2 q−1 Z − ··· − dω′D(ω′)Eω∗′−ωEω′−ω1+ω2−···−ωq−2+ωq−1 E E∗ . (A3) Z × ω−ω1+ω2−···−ωq−2+ωq−1 ωq−1 = dω D(ω ω + ω +ω ) q 1 2 q−1 q − ···− Z (q−1) Comparisonwith Eq. (20) shows that ∆nω−ω′ contribu- E∗ E . (A8) tionyieldsthe partofthe P(q) integrandproportionalto × −ω+ω1−ω2+···−ωq−1+ωq ωq ω f ∗ (q−1) D(ω ω + +ω ). When integrating ∆n , This completes the proof of Eqs. (20) and (21). 1− 2 ··· q−2 ω′−ω (cid:20) (cid:21) f Appendix B: General expressions for diagrams The diagrams X0 (q =3,5,...) are given by the following analytical expression: qj Xq0j =e |El2(r,t)|2|El4(r,t)|2···|Elq−1(r,t)|2 XXeq0j e D (ω ω )D (ω ω +ω ω ) D (ω ω + +ω ω ) × k l1 − l2 k l1 − l2 l3 − l4 ··· k l1 − l2 ··· lq−2 − lq−1 [D(ω )+D∗(ω )][D(ω ω +ω )+D∗(ω ω +ω )] × l1 l2 l1 − l2 l3 l2 − l1 l4 ··· [D(ω ω +ω ω +ω )+D∗(ω ω +ω ω +ω )]. (B1) × l1 − l2 l3 −···− lq−3 lq−2 l2 − l1 l4 −···− lq−4 lq−1 The symbol denotes a summation over the lasing-mode indices l ,l ,...,l =1,...,N according to the rules: 2 4 q−1 m XXeq0j (i) if the indices l and l are connected in the diagram X0 , set l = l and (ii) the terms with four or more indices i j qj i j l ,l ,...,l ,l equal are excluded unless the links connecting the affected vertices do not intersect. 1 2 q WedenotebyX0 connecteddiagramsandsubdiagramsecontainingthevertexl. WiththeorderingofX0 suchthat qj qj e