Optoelectronic device simulations based on macroscopic Maxwell-Bloch equations
OOptoelectronic Device Simulations based on Macroscopic Maxwell-BlochEquations
Christian Jirauschek , a) Michael Riesch , and Petar Tzenov
Department of Electrical and Computer Engineering, Technical University of Munich, Arcisstr. 21, 80333 Munich,Germany (Dated: 9 June 2020, published as Adv. Theory Simul. 2, 1900018 (2019))
Due to their intuitiveness, flexibility and relative numerical efficiency, the macroscopic Maxwell-Bloch (MB)equations are a widely used semiclassical and semi-phenomenological model to describe optical propagationand coherent light-matter interaction in media consisting of discrete-level quantum systems. This reviewfocuses on the application of this model to advanced optoelectronic devices, such as quantum cascade andquantum dot lasers. The Bloch equations are here treated as a density matrix model for driven quantumsystems with two or multiple discrete energy levels, where dissipation is included by Lindblad terms.Furthermore, the one-dimensional MB equations for semiconductor waveguide structures and optical fibersare rigorously derived. Special analytical solutions and suitable numerical methods are presented. Due to theimportance of the MB equations in computational electrodynamics, an emphasis is placed on the comparisonof different numerical schemes, both with and without the rotating wave approximation. The implementationof additional effects which can become relevant in semiconductor structures, such as spatial hole burning,inhomogeneous broadening and local-field corrections, is discussed. Finally, links to microscopic models andsuitable extensions of the Lindblad formalism are briefly addressed.Keywords: Maxwell-Bloch equations, Lindblad equation, Quantum dots, Quantum cascade laser, WaveguidepropagationThis is the peer reviewed version of the following article: C. Jirauschek, M. Riesch, and P. Tzenov,“Optoelectronic device simulations based on macroscopic Maxwell-Bloch equations”, Adv. Theory Simul. 2,1900018 (2019), which has been published in final form at https://doi.org/10.1002/adts.201900018 .This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions forUse of Self-Archived Versions.
CONTENTS
List of Symbols and Acronyms I. Introduction II. Lindblad equation
III. Optical Bloch Equations a) Electronic mail: [email protected]
D. Non-Redundant Density MatrixRepresentation 13E. Rotating Wave Approximation 14
IV. Maxwell-Bloch Equations
V. Reduction to One-Dimensional Model a r X i v : . [ c ond - m a t . m e s - h a ll ] J un
2. Spatial Hole Burning 28
VI. Analytical Solutions
VII. Numerical Schemes
VIII. Inclusion of Further Effects
IX. Application to Optoelectronic Devices
X. Conclusion and Outlook Conflicts of Interest Supporting Information Acknowledgments Biographies LIST OF SYMBOLS AND ACRONYMS
Symbol Description A Magnetic vector potential A cv Spontaneous emission coefficient A eff Effective mode area A q Quantum active region cross section a Power loss coefficient c Vacuum speed of light D Diffusion coefficient D Displacement field D Dissipation superoperator D ijmn Tensor element of D ˆd Dipole operator d ij Dipole matrix element vector E Electric field E Slowly varying amplitude of E E E = Ee E ± E for forward/backward propagating field E g Bandgap energy E i Eigenenergy of level iE ωp Fourier amplitude of waveguide E-field E t Transverse dependency of E EIT Electromagnetically induced transparency e Elementary charge e Polarization direction of electric field F Modal field distribution F i Wavefunction of semiconductor state i , F i = ϕ i u v i FDTD Finite difference time-domain f ± Forward/backward normalized polarization amplitude g Power gain coefficient g ( ω ) Distribution function of resonance frequencies H Magnetic field H Slowly varying amplitude of H ˆ H Hamiltonian ˆ H Unperturbed system Hamiltonian ˆ H Interaction Hamiltonian H ωp Fourier amplitude of waveguide H-field H t Transverse dependency of H (cid:126) Reduced Planck constant I Optical intensity I s Saturation intensity J f Current density due to free carriers J q Total current density due to quantum systems k
3D bulk/2D in-plane wavevector k ω/ck c Carrier wavenumber with | k c | = n ω c /c L Liouville superoperator ˆ L Lindblad operator ˆ L α → β ˆ L for incoherent transition (cid:96) v g a/ MB Maxwell-Bloch m Mass m ∗ Effective mass m e Electron mass N Number of system levels n Complex refractive index n Background refractive index n QD sheet density n Carrier number density n cv Number density of electron-hole pairs n eff Effective waveguide index β/k n eff Complex effective waveguide index β/k P Optical power P Macroscopic polarization P Slowly varying amplitude of P q P ωp Fourier amplitude of quantum system polarization P q Contribution of P due to quantum systems P s Saturation power ˆp Momentum operator p s Carrier fraction for inhomogeneous broadeningQCL Quantum cascade laserQD Quantum dot q Carrier chargeRNFD Risken-Nummedal finite differencesRWA Rotating wave approximation r Microscopic position vector ˆr Position operator r i Outscattering rate from level ir i → j Transition rate from level i to jS In-plane area of quantum well S Coherence vector, Bloch vector [ u, v, w ] T S z z component of Poynting vectorSHB Spatial hole burningSVAA Slowly varying amplitude approximation s Complex frequency variable s Bloch vector in RWA T Energy relaxation time T Phase relaxation rateTE Transverse electricTM Transverse magnetic t Time variable u Bloch vector component u v Periodic Bloch function of band vV Potential energy V p Probe volume v Bloch vector component ˆv Velocity operator v g Group velocity W γαβ, k Boltzmann rates w Population inversion x Coordinate, e.g., in growth direction x Macroscopic position vector y In-plane coordinate z Propagation coordinate β (cid:60) (cid:8) β (cid:9) β Complex propagation constant β n [d nω β ] ω = ω c Γ Overlap factor Γ ij Pumping rate from level i to jγ Self-phase modulation coefficient γ Energy relaxation rate γ Phase relaxation rate γ ij Dephasing rate between levels i and jγ (cid:48) ij Pure dephasing rate between levels i and j ∆ Frequency detuning ∆ in two-level system ∆ ij Frequency detuning sgn ( ω ij ) ( ω c − | ω ij | )∆ n Transverse refractive index profile ∆ t Time step ∆ z Spatial increment ∆ ω ω − ω c (cid:15) Vacuum permittivity (cid:15) r Dielectric constant (cid:15) r Complex dielectric constant η ij Slowly varying envelope function of ρ ij η ± ij η ij for forward/backward propagating field µ Vacuum permeability ˆ ρ Density operator ρ ± ii Population grating amplitude, ρ + ii = (cid:0) ρ − ii (cid:1) ∗ ρ ii Center value of population grating ρ ij Density matrix element σ Conductivity ϕ i Envelope wavefunction of level i Ψ Time dependent wavefunction ψ i
1D wavefunction of level i in quantum wells Ω Instantaneous Rabi frequency Ω Slowly varying amplitude of ΩΩ g Generalized Rabi frequency ω Frequency variable ω c Carrier frequency ω ij Transition frequency between levels i and j I. INTRODUCTION
Due to advancements in nanotechnology, structuringin the nanometer range is meanwhile routinely exploitedin electronics and photonics. For example, in optoelec-tronic devices such as semiconductor optical amplifiersand lasers, quantum confinement is widely used to con-centrate the carriers in certain energy states, yieldingimproved wall-plug efficiencies and higher output pow-ers. As a further effect, the wavelength can be tunedby changing the size of the confinement structure. Ona commercial basis, mostly one-dimensional confinementis used in form of quantum well structures, which arefabricated based on deposition of nanometer-thin semi-conductor layers of different compositions. In such struc-tures, a quantum well is formed by a layer consisting of alower bandgap material than the adjacent layers, whichrestricts the free electron motion in that layer to the in-plane directions and gives rise to quantized energy statesin growth direction. As a consequence of the further re-striction of the energy spectrum and the even strongercarrier localization, additional improvement can be ex-pected from two- or three-dimensional confinement, re-sulting in quantum wire/dash and quantum dot (QD)structures, respectively. Indeed, QD and quantumdash lasers and laser amplifiers have been shown to ex-hibit excellent characteristics. In Fig. 1, the formationof quantized states in quantum wells, wires and dots isschematically illustrated. The term quantum dash refersto an elongated nanostructure, i.e., some kind of shortquantum wire. By contrast, the term nanowire does notnecessarily indicate strong quantum confinement. Forexample, in nanowire lasers the nanowire geometry typi-cally serves as a single-mode optical waveguide resonator,while the active region is based on a heterostructure orquantum well, as in a conventional laser diode. Semiconductor optoelectronic devices usually rely onelectron-hole recombination, i.e., optical transitions be-tween conduction and valence band states. The associ-ated resonance wavelength is largely determined by thesemiconductor bandgap, which establishes a lower boundon the transition energy. Thus, the coverage of a cer-tain spectral region depends on the existence of suitablesemiconductor materials, which for example restricts theavailability of practical optoelectronic sources and de-
FIG. 1. Schematic illustration of quantum confinement struc-tures: (a) Quantum well, (b) quantum wire, (c) quantum dot. tectors in the mid-infrared and terahertz regions. Analternative concept is based on intersubband devices,which employ so-called intersubband transitions betweenquantized energy states in the conduction (or, in somecases, valence) band of a nanostructure, and thus al-low quantum engineering of the transition wavelengthindependent of the bandgap. Quantum well devicesbased on this concept include quantum cascade lasers, quantum cascade detectors and quantum well infraredphotodetectors. Furthermore, intersubband transitionsare used for QD infrared photodetectors.
Along with quantum confinement, also quantum co-herent effects are found to be increasingly relevant formodern optoelectronic devices. Such effects result fromthe coherent light-matter interaction, which requires thatthe states involved in the optical transition maintain awell-defined phase relationship over a significant time.The coherent interaction manifests itself in so-called Rabiflopping, i.e., carrier population oscillations betweenthe states, which are driven by the optical field. Theresulting carrier dynamics couples back to the opticalfield via the polarization, thus also affecting the propa-gating optical waveform. Besides being an essential pre-requisite for the emerging field of quantum informationtechnology, quantum coherence plays an increasinglyimportant role for modern optoelectronic devices in gen-eral. Due to the strong interaction with the semiconduc-tor environment, e.g., in the form of phonon scatteringand carrier-carrier interactions, this phase relationshiptends to be quickly destroyed, which is commonly re-ferred to as dephasing. However, under favorable condi-tions, signatures of Rabi oscillations have been observedin nanostructured optoelectronic systems and devices.These include quantum well structures, nanowirelasers, quantum cascade lasers, and single QDs atcryogenic temperatures, as well as QD and quantumdash amplifiers at room temperature. Closely relatedis self-induced transparency, where Rabi flopping en-ables a special optical pulse form to propagate withoutbeing attenuated or disturbed. This phenomenon hasmeanwhile also been observed in semiconductor struc-tures such as QD waveguides, with potential appli-cations such as the generation of ultrashort optical pulsesin QD and quantum cascade lasers. Another effectthat relies on quantum coherence is slow light propaga-tion or even complete halting of light, with possi-ble applications such as optical buffers, imaging and quantum memory. This effect has meanwhile alsobeen demonstrated in solid-state media, namely in dopedcrystals.
The use of suitably engineered semiconduc-tor structures would be especially attractive from a prac-tical point of view.
Furthermore, quantum interfer-ence, e.g., in QD or intersubband quantum well systems,is an interesting candidate to realize all-optical switch-ing.
Due to the discrete energy level structure ofQDs, semiconductor devices based thereupon are espe-cially likely to be, at least in part, governed by coher-ence effects, although for QD ensembles the dephas-ing tends to be strong. The same applies to intersub-band quantum well devices, where the levels close to theband edge have parallel dispersion relations, and thusthe quantum dynamics resembles that of discrete-levelsystems. Especially in such devices, coherence effectscan significantly influence the dynamic operation evenfor considerable dephasing.Suitable theoretical models are required for an in-depthunderstanding of the often quite complex interplay of ef-fects determining the dynamic device characteristics, aswell as for quantitative simulation and systematic deviceoptimization. For the optoelectronic devices and struc-tures discussed above, an adequate theoretical descrip-tion must include the coherent carrier dynamics, incoher-ent processes such as scattering or spontaneous emission,as well as the interaction with the optical field. Our fo-cus is here on the very widely used Maxwell-Bloch (MB)equations. The Bloch equations provide a compact modelfor the discrete-level carrier dynamics, which is describedby the density matrix formalism. The coherent single-carrier dynamics is here modeled by the Hamiltonian ofthe quantum system, such as a QD, and also includesthe interaction with a classical optical field. Effects be-yond the single-electron quantum evolution are regardedas interaction with the environment in form of the semi-conductor host, which gives rise to incoherent effects suchas scattering with other carriers and phonons. The re-sulting dissipation in the quantum system is in the Blochequations phenomenologically modeled by relaxation rateterms, which introduce dephasing and incoherent carriertransitions. The Bloch equations were first devised to de-scribe the evolution of the nuclear magnetic moment in amagnetic field, and later on extended to a pair of levelsin resonance with a classical optical field. The modelis closed by coupling the Bloch equations to Maxwell’sequations, which describe the evolution of the clas-sical optical field. This review paper is concerned withthe resulting MB equations, where we go beyond the of-ten applied two-level approximation by consider-ing multiple, albeit discrete, energy levels. Furthermore,we root the phenomenological dissipation terms in theLindblad formalism, which ensures physical behavior ofthe quantum system and allows for the construction ofmore general dissipation terms.
The MB equations offer a generic description of semi-classical light-matter interaction, which can be appliedto different media such as semiconductor structures orgases. The focus of this review lies on semiconduc-tor structures, which is reflected in the treatment ofsome specific issues, such as the concrete embodimentof Maxwell’s equations, or the inclusion of spatial holeburning in linear resonators. Independent of the modeledsystem, the main attractiveness of the MB equations liesin the relatively compact description of the carrier dy-namics, which is helpful for providing intuitive insightinto the device behavior and even allows for closed ana-lytical solutions in some special cases.
From a com-putational point of view, the Bloch equations are widelyused in combination with electromagnetic simulations,e.g., based on the finite-difference time-domain method,as a quantum model of the medium, replacing simplerclassical descriptions such as the Lorentz model. Due tothe relative compactness of the Bloch model, also compu-tationally demanding two- or three-dimensional simula-tions can be carried out. Likewise, the MB equationsenable systematic device optimizations over a large pa-rameter range, as well as long-term simulations, e.g., toinvestigate the steady-state laser dynamics.
An-other important advantage of the MB equations is thatthey can easily be adapted to specific problems by addingfurther effects, such as inhomogeneous broadening orlocal-field corrections.
Clearly, the Bloch equations constitute a compromisebetween accuracy and compactness of the model. A fullmicroscopic treatment of light-matter interaction in asemiconductor, accounting for carrier-phonon and many-body Coulomb interactions as well as for free carrier mo-tion in the unconfined directions, results in the so-calledsemiconductor MB equations.
These do not requirephenomenological input parameters, but the significantlyincreased model complexity usually restricts the model-ing to one spatial dimension and short-term simulations.While the semiconductor MB equations are beyond thescope of this review, they can be used as a basis to derivemacroscopic discrete-level MB equations with Lindbladdissipation and additional correction terms for specificsemiconductor structures.
In detail, our paper is organized as follows: In SectionII, the density matrix formalism and Lindblad model areintroduced, serving as a basis for the Bloch equations.These are treated in Section III, which also includes a dis-cussion of the widely used rotating wave approximation(RWA). In Section IV, the MB equations are introducedin full-wave treatment and invoking the RWA, along withthe slowly varying amplitude approximation (SVAA) forthe field propagation. Section V treats the reduction ofthe MB equations for semiconductor waveguide struc-tures and optical fibers to a spatially one-dimensionalmodel, which is a widely used simplification. Section VIdeals with available analytical solutions for the Bloch andMB equations, while in Section VII, numerical methodsfor the MB equations are covered. Section VIII is dedi-cated to the inclusion of further effects, such as local-fieldcorrections, inhomogeneous broadening and noise. Sec-tion IX deals with the application of the MB model to concrete optoelectronic devices, including bulk as well asinter- and intraband quantum well and QD devices. Thepaper is concluded in Section X, where dissipation mod-els beyond the Lindblad formalism are discussed.
II. LINDBLAD EQUATION
In the following, we consider discrete quantum systemswith N states | i (cid:105) , where i = 1 ..N . We restrict ourselvesto a single-particle description, valid for carrier densi-ties which are sufficiently low to neglect Pauli blocking,but sufficiently high to neglect electron-hole Coulombcorrelation. It has been pointed out that these require-ments are often fulfilled in state-of-the-art semiconductorquantum devices which are the main scope of this paper,and that the Lindblad approach introduced below is thenwell justified. Furthermore, we do not explicitly con-sider spin dependent effects, even though the Lindbladformalism can be extended accordingly. The time evolution of an ideal quantum system isfamously described by the time dependent Schrödingerequation i (cid:126) ∂ t | Ψ (cid:105) = ˆ H | Ψ (cid:105) (1)with the reduced Planck constant (cid:126) , where the systemstate vector | Ψ (cid:105) , and generally also the system Hamil-tonian ˆ H , depend on time t . | Ψ (cid:105) is a pure state, i.e.,a coherent superposition of the basis states | i (cid:105) with | Ψ ( t ) (cid:105) = (cid:80) i c i ( t ) | i (cid:105) , where c i are complex coefficients.In reality, however, no quantum system is perfectly iso-lated, but rather interacts with its environment. Thisinduces decoherence, i.e., loss of quantum coherence inthe system, which must be included into any realistic de-scription. The resulting statistical state of the systemis generally a mixed state which cannot be representedby the system state vector | Ψ (cid:105) , but rather requires anextended description in terms of the density operator ˆ ρ .The corresponding density matrix with respect to thechosen basis states | i (cid:105) has the elements ρ ij = (cid:104) i | ˆ ρ | j (cid:105) ,where the diagonal elements ρ ii give the occupation prob-ability of state | i (cid:105) , while the off-diagonal elements ρ ij represent the coherence between | i (cid:105) and | j (cid:105) . The densityoperator is positive semidefinite which guarantees thatany pure system state | Ψ (cid:105) has a non-negative probabil-ity, i.e., (cid:104) Ψ | ˆ ρ | Ψ (cid:105) ≥ . This also implies Hermiticity,i.e., ˆ ρ = ˆ ρ † and thus ρ ij = ρ ∗ ji , where the dagger andasterisk denote the adjoint and the complex conjugate,respectively. Furthermore, at least for closed systems, the trace must remain constant to ensure particle con-servation and is usually normalized to unity, Tr { ˆ ρ } = 1 .The coherent time evolution of the density operator is inthe Schrödinger picture described by the von Neumannequation i (cid:126) ∂ t ˆ ρ = (cid:104) ˆ H, ˆ ρ (cid:105) . (2) (cid:1) I FIG. 2. Schematic illustration of a quantum system interact-ing with impurities and thermal vibrations in the semicon-ductor lattice.
In realistic scenarios, often many degrees of freedomare relevant for the time evolution and must thus be con-sidered. Usually, only part of these degrees of freedomare of direct interest for the application in mind, andsolving the full Eq. (2) is typically also too demanding.This issue can be addressed by performing a division intoa system containing the degrees of freedom which are ofprimary interest, and a second one with the remainingdegrees of freedom which then constitute the environ-ment. In semiconductor quantum devices, the degrees offreedom of interest may be quantized states in a nanos-tructure such as a quantum well or dot, while decoherencetypically arises from interaction with the semiconductorlattice itself, thus acting as the environment. This sit-uation is schematically illustrated in Fig. 2. There arevarious types of interactions, also referred to as scatter-ing mechanisms, which can induce decoherence in thequantum system of interest. These include the interac-tion with phonons due to (longitudinal- and transverse-optical and -acoustic) thermal lattice vibrations, latticeimperfections in form of impurities (such as dopants), in-terface roughness or atomic disorder in alloys, as well aspiezoelectric fields. Also carrier-carrier interaction canenter the single-particle picture as an additional scatter-ing mechanism.
For a quantum system interactingwith the environment, ˆ ρ and ˆ H in Eq. (2) refer to the fulldynamics of the combined system and environment.The Hamiltonian ˆ H can be written as ˆ H = ˆ H S ⊗ ˆ I E +ˆ I S ⊗ ˆ H E + ˆ H I , where the Hamiltonians ˆ H S , ˆ H E and ˆ H I de-scribe the system S , the environment E and the system-environment interaction, ˆ I S and ˆ I E are the unit opera-tors in the respective Hilbert spaces, and ⊗ denotes thetensor product. The reduced density matrix of the sys-tem of interest is simply obtained by tracing over theenvironmental degrees of freedom, ˆ ρ S = Tr E { ˆ ρ } . Thisstep by itself does obviously not eliminate the depen-dence of Eq. (2) on the environment. Thus, additionalassumptions are necessary to arrive at a model for thenon-unitary time evolution of ˆ ρ S , which is a consequenceof eliminating the environmental degrees of freedom. The resulting equation is expected to be similar in structureas Eq. (2), i.e., a first-order linear differential equation intime for ˆ ρ S , where the linearity ensures consistency withthe ensemble interpretation of the density matrix. Theresulting time-local and Markovian description of the re-duced density matrix dynamics is commonly referred toas (quantum) master equation. Its general form can beinferred by posing additional requirements to avoid un-physical behavior. In particular, this includes conser-vation of unit trace and positive semidefiniteness of thedensity matrix, as discussed above. Closer inspection re-veals that if there exists another system S (cid:48) , an evolutionequation for S which ensures positive semidefiniteness of ˆ ρ S can still lead to unphysical time evolution of the com-bined density matrix for S and S (cid:48) , even if S (cid:48) does notevolve and is completely decoupled from S . This prob-lem is cured by demanding complete positivity of theevolution, rather than only the preservation of positivesemidefiniteness of ˆ ρ S .From above requirements, the general form of the evo-lution equation can be inferred by invoking the Kraustheorem, characterizing completely positive trace pre-serving maps. The resulting master equation is calledLindblad equation. Dropping the subscript S fromhere on for ease of notation, it can be written as ∂ t ˆ ρ = − i (cid:126) (cid:104) ˆ H, ˆ ρ (cid:105) + (cid:88) k (cid:18) ˆ L k ˆ ρ ˆ L † k −
12 ˆ L † k ˆ L k ˆ ρ −
12 ˆ ρ ˆ L † k ˆ L k (cid:19) = L (ˆ ρ ) + (cid:88) k D k (ˆ ρ ) = L (ˆ ρ ) + D (ˆ ρ ) , (3)where ˆ H = ˆ H + ˆ H is the effective Hamiltonian of thereduced system. Here, the Hamiltonian ˆ H describesexternally induced perturbations, e.g., due to an inci-dent optical field. The description of light-matter in-teraction requires a time dependent Hamiltonian, which,although lifting the originally assumed time-homogeneityof the Lindblad equation, still gives a valid density ma-trix evolution. In addition, ˆ H may contain non-dissipative contributions stemming from the interactionwith the environment, such as energy shifts. The dis-sipation is described by the sum term, where the linearoperators ˆ L k are called Lindblad or (quantum) jump op-erators, which can in principle be chosen without furtherrestrictions in the Hilbert space of the reduced system.Equation (3) now includes both the coherent dynamicsdue to the Liouville superoperator L (ˆ ρ ) , correspondingto Eq. (2), and the incoherent dynamics induced by thedissipation superoperator D (ˆ ρ ) which contains the inter-action with the environment. Besides inferring the Lind-blad equation from the requirements given above, Eq. (3)can also be microscopically derived, assuming that thequantum system is weakly coupled to a large Markovianenvironment. As mentioned above, we allow for a time dependentHamiltonian in Eq. (3), which is required to include light-matter interaction as envisaged in this paper, and con-stitutes a slight generalization of the original equation. Occasionally, also time dependent Lindblad operators ˆ L k ( t ) are used, for example to model time dependentpumping rates. This also does not affect the physicalvalidity of Eq. (3), since conservation of trace and com-plete positivity are further guaranteed.
Moreover,Eq. (3) with time dependent operators ˆ H and ˆ L k is stilltime-local and also Markovian. A. Introduction of Basis States
In principle, the N basis states of the (reduced) quan-tum system can be freely selected as long as they spanthe entire Hilbert space of the N -level system. In mostcases, an orthonormal basis is the preferred option, sinceit results in more compact expressions and provides aclearer physical interpretation. The choice of energyeigenstates has the distinct advantage that the reducedsystem Hamiltonian ˆ H is diagonal. In certain cases,other choices may be preferable, such as a localized (ortight-binding) basis set for the description of tunneling,e.g., in double- or multiple-well systems. Assuming an orthonormal basis so that the unit op-erator becomes ˆ I = (cid:80) N(cid:96) =1 | (cid:96) (cid:105) (cid:104) (cid:96) | , Eq. (3) can be writtenas ∂ t ρ ij = − i (cid:126) (cid:88) (cid:96) ( H i(cid:96) ρ (cid:96)j − H (cid:96)j ρ i(cid:96) )+ (cid:88) k (cid:88) (cid:96)p (cid:20) L i(cid:96)k (cid:16) L jpk (cid:17) ∗ ρ (cid:96)p − (cid:0) L (cid:96)ik (cid:1) ∗ L (cid:96)pk ρ pj − (cid:16) L p(cid:96)k (cid:17) ∗ L pjk ρ i(cid:96) (cid:21) = − i (cid:126) (cid:88) mn L ijmn ρ mn + (cid:88) mn D ijmn ρ mn . (4)Here, H ij = (cid:104) i | ˆ H | j (cid:105) and L ijk = (cid:104) i | ˆ L k | j (cid:105) are the matrixelements of the operators ˆ H and ˆ L k . Also the superop-erators can be represented in form of a matrix, albeit ofsize N × N , with elements L ijmn = H im δ jn − H nj δ im , (5) D ijmn = (cid:88) k (cid:26) L imk (cid:16) L jnk (cid:17) ∗ − (cid:88) (cid:96) (cid:104)(cid:0) L (cid:96)ik (cid:1) ∗ L (cid:96)mk δ jn + (cid:0) L (cid:96)nk (cid:1) ∗ L (cid:96)jk δ im (cid:105) (cid:27) , (6)where δ denotes the Kronecker delta. We emphasize thatwhile Eq. (6) ensures that there is a matrix representa-tion D ijmn for any given set of Lindblad operators, theconverse is not necessarily true, and arbitrarily chosen D ijmn can produce unphysical results. B. Choice of Lindblad Operators
The choice of the ˆ L k for generating a certain time evo-lution is not unique. In particular, for a given set ˆ L k with k = 1 , . . . , K , the set ˆ L (cid:48) k = (cid:80) (cid:96) u k(cid:96) ˆ L (cid:96) (also with (cid:96) = 1 , . . . , K ) generates the same dynamics for an ar-bitrary unitary matrix with dimension K and elements u k(cid:96) . This can easily be verified by substituting the ˆ L k in Eq. (3) with above expression for ˆ L (cid:48) k , and con-sidering that (cid:80) k u k(cid:96) u ∗ km = δ (cid:96)m . Furthermore, the ˆ L k might also contain unitary contributions, which can al-ternatively be included into the Hamiltonian ˆ H . In par-ticular, replacing an operator ˆ L k by ˆ L (cid:48) k = ˆ L k + α k ˆ I where α k is an arbitrary complex constant with dimen-sion of inverse square root of time, and ˆ H by ˆ H (cid:48) = ˆ H +(i (cid:126) / (cid:16) α k ˆ L † k − α ∗ k ˆ L k (cid:17) generates the same dynamics. From the Kraus theorem it follows that it is alwayspossible to choose the Lindblad operators so that a givennon-unitary evolution can be represented by K ≤ N − operators [in addition to ˆ L ∝ ˆ I which gives a vanishingcontribution in Eq. (3)]. Formally, such a representationcan be constructed by starting from the Kossakowski–Sudarshan form of the Lindblad equation and apply-ing a unitary transformation to convert it to Eq. (3). However, it has been pointed out that the resulting stan-dard form does not give much insight into the underlyingphysical processes. From a practical point of view, itis more natural to choose the ˆ L k so that they representcertain physical effects. In the following, we will discussthe two most relevant mechanisms, i.e., incoherent transi-tions between states corresponding to hopping transport,and pure dephasing which affects the coherence betweentwo states but does not involve population transfer be-tween them.
1. Incoherent Transitions
For a transition from a given basis state | α (cid:105) to | β (cid:105) witha rate r α → β , the associated Lindblad operator is given by ˆ L α → β = r / α → β | β (cid:105) (cid:104) α | , (7)and Eq. (6) for the corresponding superoperator matrixelements yields D α → βijmn = r α → β (cid:20) δ iβ δ jβ δ mα δ nα −
12 ( δ iα δ jn δ mα + δ im δ jα δ nα ) (cid:21) . (8)Inserting Eq. (7) into Eq. (4), we obtain popula-tion changes [ ∂ t ρ ββ ] α → β = r α → β ρ αα , [ ∂ t ρ αα ] α → β = − r α → β ρ αα . The population relaxation is thus generallydescribed by rate equation terms [ ∂ t ρ αα ] relax = (cid:88) j (cid:54) = α r j → α ρ jj − r α ρ αα , (9)where r α = (cid:88) j (cid:54) = α r α → j (10)is the total outscattering rate from level α . Further-more, we see that apart from the population changes, ˆ L α → β also contains the associated lifetime contributionto dephasing, with [ ∂ t ρ αn ] α → β = − ( r α → β / ρ αn and [ ∂ t ρ nα ] α → β = − ( r α → β / ρ nα where n (cid:54) = α . This meansthat population transfer from a state | α (cid:105) to | β (cid:105) inducesdephasing not only for this transition, but also for othertransitions involving | α (cid:105) , and ignoring this fact mightlead to unphysical results. On the other hand, this im-plies that the total lifetime contribution to the dephas-ing rate for a transition α → β is with Eq. (10) given by ( r α + r β ) / , i.e., is obtained from the total outscatteringrates for levels α and β . We note that the operator inEq. (7) provides an elementary description of transitions,but does for example not take into account correlationsbetween different transition processes.
2. Pure Dephasing
In addition to above discussed population changes,there can be additional mechanisms which do not in-volve population transfer between the chosen basis states,but cause additional decoherence, resulting in a decay ofoff-diagonal density matrix elements only.
This so-called pure dephasing contribution between two levels α and β (cid:54) = α can be described as [ ∂ t ρ αβ ] pure = − γ (cid:48) αβ ρ αβ ,which also implies [ ∂ t ρ βα ] pure = − γ (cid:48) αβ ρ βα since ρ αβ = ρ ∗ βα . Here, γ (cid:48) αβ = γ (cid:48) βα ≥ denotes the pure dephasingrate. As can easily be seen, the corresponding dissipa-tion superoperator in Eq. (4) can be represented by thematrix elements D αβijmn = − γ (cid:48) αβ ( δ iα δ jβ δ mα δ nβ + δ iβ δ jα δ mβ δ nα ) . (11)The Lindblad operators for pure dephasing must bediagonal in the chosen basis. However, Eq. (11) doesnot generally ensure physical behavior, and thus a repre-sentation in terms of Lindblad operators does not alwaysexist. Notably, for N ≥ there are constraints on howto select the pure dephasing rates γ (cid:48) αβ ≥ to ensurecompatibility with Eq. (3), and an ill-considered choicecan for example easily result in a violation of positivesemidefiniteness for ˆ ρ . For example, γ (cid:48) + γ (cid:48) + γ (cid:48) ≤ γ (cid:48) γ (cid:48) + γ (cid:48) γ (cid:48) + γ (cid:48) γ (cid:48) ) / must hold in three-levelsystems, which is already violated if only one of the threepure dephasing rates is non-zero.In two-level systems, pure dephasing is described by asingle rate γ (cid:48) = γ (cid:48) = γ (cid:48) ≥ , and can for example berepresented by a Lindblad operator ˆ L = (2 γ (cid:48) ) / | (cid:105) (cid:104) | or ˆ L = (2 γ (cid:48) ) / | (cid:105) (cid:104) | , or also by the set ˆ L = ( γ (cid:48) ) / | (cid:105) (cid:104) | , ˆ L = ( γ (cid:48) ) / | (cid:105) (cid:104) | . More generally, if the same (typically empirical) pure dephasing rate γ (cid:48) is assumed for all tran-sitions of an N -level system, this case can always berepresented by Lindblad operators, for example by theset ˆ L k = ( γ (cid:48) ) / | k (cid:105) (cid:104) k | , k = 1 ..N . Taking into account the results of Section II B 1, thetotal phase relaxation due to pure dephasing plus life-time broadening associated with incoherent transitionsis described by the dissipation term [ ∂ t ρ αβ ] relax = − γ αβ ρ αβ = − (cid:2) ( r α + r β ) / γ (cid:48) αβ (cid:3) ρ αβ , (12)where γ αβ = ( r α + r β ) / γ (cid:48) αβ is the total dephasingrate and the r α,β are given by Eq. (10).
3. General Case
While physical dissipation channels can often berepresented by either incoherent transitions or puredephasing, see Sections II B 1 and II B 2, the Lind-blad operators should not a priori be restricted tothese two forms, but rather be found based on physi-cal considerations. Even more, the representation ofa dissipative channel as, e.g., incoherent transition orpure dephasing, only applies for the chosen basis.
For illustration, let’s assume an N -level system with or-thonormal basis states | n (cid:105) and dissipative channels de-scribed by a set of Lindblad operators ˆ L k . Alternatively,an orthonormal basis with states | n (cid:48) (cid:105) can be used, with | n (cid:105) = ˆ I | n (cid:105) = (cid:80) n (cid:48) (cid:104) n (cid:48) | n (cid:105) | n (cid:48) (cid:105) , which changes the char-acter of the Lindblad operators in the new basis system.As an illustrative example, we restrict ourselves to tworelevant levels | (cid:105) and | (cid:105) , which are assumed to be lo-calized in adjacent potential wells, and between whichtunneling through the separating barrier occurs. Thismechanism plays for example an important role in QCLs,which are frequently modeled with a density matrix ap-proach for a discrete quantum system, using localizedstates to describe the tunneling transport across thickbarriers. This tunneling process is criticallyaffected by dephasing between the two states involved,which can be modeled by Eq. (12).
We exem-plarily focus on the pure dephasing contribution, whichcan for a two-level system be described by the Lind-blad operator ˆ L = (2 γ (cid:48) ) / | (cid:105) (cid:104) | as discussed in Sec-tion II B 2. Changing to energy eigenstates | (cid:48) (cid:105) and | (cid:48) (cid:105) and for simplicity assuming near-degeneracy, we obtain | (cid:105) = 2 − / ( | (cid:48) (cid:105) + | (cid:48) (cid:105) ) and | (cid:105) = 2 − / ( | (cid:48) (cid:105) − | (cid:48) (cid:105) ) . Inthe energy basis, above Lindblad operator then becomes ˆ L = ( γ (cid:48) ) / ( | (cid:48) (cid:105) (cid:104) (cid:48) | + | (cid:48) (cid:105) (cid:104) (cid:48) | + | (cid:48) (cid:105) (cid:104) (cid:48) | + | (cid:48) (cid:105) (cid:104) (cid:48) | ) whichis not diagonal, i.e., does not represent pure dephasingin that basis.To summarize, the frequently used classification of dis-sipation channels in incoherent transitions and pure de-phasing is not always possible and additionally dependson the chosen basis system, but is frequently used since itallows for an intuitive physical interpretation. Thus, thisclassification might also be helpful for determining thecorresponding dissipative rates based on compact modelsor by comparison to experimental data. Conse-quently, for a given system a criterion for a convenientchoice of basis states might be that the dissipation chan-nels can reasonably well be described in terms of incoher-ent transitions and pure dephasing, which for examplemotivates the frequent use of localized states to describetunneling transport through thick barriers.
C. Conditions for Validity
As discussed in Section II B, the dissipation parame-ters must fulfill certain conditions to ensure physical be-havior of the density matrix, which is exactly true if arepresentation of the dissipation process in terms of Lind-blad operators exists. For example, the total dephasingrate of a given transition cannot be smaller than the life-time broadening contribution due to incoherent transi-tions, as can be seen from Eq. (12). Also, as discussedin Section II B 2, the pure dephasing rates cannot be in-dependently chosen for each transition, but must fulfillcertain conditions for N ≥ levels. Thus, if the exper-imentally obtained dissipation rates for a system do notsatisfy above conditions, this might indicate that the cho-sen model is not adequate, for example that not enoughlevels are considered. As noted above, the Lindblad equation can also bemicroscopically derived for a quantum system weaklycoupled to a large Markovian environment.
Theseassumptions require in particular that the coherentsystem dynamics and relaxation processes occur ona slower timescale than the memory decay of theenvironment.
These additional microscopic con-straints are not required to ensure completely positiveand trace preserving evolution of the density matrix,which is guaranteed by the Lindblad form of Eq. (3).However, disregarding the microscopic validity criteriamight result in a violation of other laws such as On-sager’s relation. On the other hand, it has been pointedout that some of the assumptions usually invoked in mi-croscopic derivations, such as the secular approximation,might be unnecessarily restrictive. Eventually, for a de-scription of realistic quantum systems where many de-grees of freedom affect the time evolution, there will al-ways be a trade-off between exactness and manageabilityof the model. From a practical point of view, Lindblad-type master equations, such as the MB system, oftenstill yield useful results on the verge of the microscopicvalidity range, for example in semiconductor structuresinteracting with high-intensity fields.
III. OPTICAL BLOCH EQUATIONS
The most basic quantum system is the two-level sys-tem with only N = 2 relevant states. This can be anatural two-level system with only two eigenstates such as a spin 1/2 particle, or a quasi-two-level system withtwo strongly coupled states, such as an optical transi-tion in resonance with an electromagnetic field or adriven double-well potential. In an early application ofthis model, Rabi investigated the interaction of a spin1/2 particle with a rotating magnetic field by solving thetime dependent Schrödinger equation. The term ”Blochequations”, in the narrow sense, refers to evolution equa-tions for a dissipative two-level system, first devised todescribe the evolution of the nuclear magnetic momentin a magnetic field. Here, the interaction with the en-vironment was taken into account by two phenomeno-logical relaxation time constants. This concept was ex-tended to other two-level systems, such as a pair of levelsin resonance with a classical optical field.
The re-sulting evolution equations are occasionally called opti-cal Bloch equations for distinction. The optical prop-agation can be considered by coupling the Bloch modelto Maxwell’s equations, resulting in the so-calledMaxwell-Bloch (MB) equations. In the following, we fo-cus on the interaction of a quantum system with an op-tical field, where the coupled MB equations have to beused for a combined description of the system dynam-ics and optical propagation. Here, we will not restrictourselves to two-level systems, but rather consider themore general case of N ≥ discrete levels. The resultingequations are for N ≥ states occasionally also referredto as a multilevel Bloch/MB model. Furthermore, forthe description of dissipative effects due to the system in-teraction with the environment, the Lindblad formalismintroduced in Section II will serve as a framework. Some-times the Lindblad equation, Eq. (3), is already referredto as Bloch equations.
In the following, the (optical)Bloch equations will be regarded as a special form ofEq. (3) containing an interaction Hamiltonian ˆ H ( t ) todescribe light-matter coupling. A. Dipole Approximation
We consider a Hamiltonian of the form ˆ H =( ˆp − q A ) / (2 m ) + qϕ + V , which models the system’sinteraction with a classical optical field, represented bya time and space dependent magnetic vector potential A = A ( ˆr , t ) and electric potential ϕ = ϕ ( ˆr , t ) . Here, ˆr and ˆp denote the position and (canonical) momentumoperators of the quantum system with the commutator [ˆ r i , ˆ p j ] = i (cid:126) δ ij , which are in position representation givenby ˆr = r and ˆp = − i (cid:126) ∇ , and V = V ( ˆr ) representsthe system’s potential energy. Furthermore, m and q denote the carrier mass and charge, which are for elec-trons given by m = m e and q = − e , with the elementarycharge e . Using the Coulomb gauge ∇ A = 0 , we have [ A , p ] = 0 . Furthermore assuming a radiation field with-out free charge contributions gives ϕ = 0 , and E = − ∂ t A for the corresponding electric field. Under these as-sumptions, we obtain ˆ H = ˆ H + ˆ H with the Hamiltonianof the unperturbed system ˆ H = ˆp / (2 m ) + V , and the0time dependent interaction Hamiltonian ˆ H = − ( q/m ) Aˆp + q A / (2 m ) . (13)The Bloch equations are then obtained from Eq. (4) bychoosing the energy eigenstates of the system Hamilto-nian ˆ H as basis, resulting in matrix elements H ,ij = E i δ ij where E i is the eigenenergy of state i , and H ,ij = − qm (cid:104) i | Aˆp | j (cid:105) + q m (cid:104) i | A | j (cid:105) . (14)Typically, the field varies on the scale of the optical wave-lengths involved, and the system dimensions are muchsmaller. The carriers do then not experience a spatialfield variation across the quantum system, and A inEqs. (13) and (14) can be represented by a space inde-pendent vector potential, evaluated at the macroscopicposition of the quantum system. In this case, it can beshown by a gauge transformation that the interactionHamiltonian in Eq. (13) is equivalent to ˆ H = − ˆdE ( t ) , (15)which corresponds to the interaction Hamiltonian in thewidely used (electric) dipole approximation. Here, ˆd = q ˆr denotes the system’s dipole operator, and theelectric field E ( t ) is taken at the system position. In-tuitively, the Hamiltonian in Eq. (15) corresponds to thepotential energy associated with the force q E ( t ) exertedby the electric field on the carriers.We note that under some special conditions, such ashigh harmonic generation or strong plasmonic confine-ment in nanophotonic structures, the field gradientmay become so large that the dipole approximation isnot applicable. In this context, we re-emphasize thatEq. (15) only assumes a spatially constant field withina given quantum system, but does not neglect the term ∝ A in Eq. (13) and is thus not restricted to weak fields,as is sometimes believed. For the interaction Hamilto-nian in Eq. (15), the mechanical and canonical momen-tum operators coincide, m ˆv = ˆp . The Hamiltonian ofthe unperturbed system ˆ H = ˆp / (2 m ) + V thus corre-sponds to the instantaneous energy operator, and a ma-trix element ρ ii in the eigenstate basis | i (cid:105) of H can beinterpreted as the measurable probability of finding thesystem in the corresponding energy eigenstate. For theinteraction Hamiltonian in Eq. (13), the mechanical mo-mentum operator is m ˆv = ˆp − q A . This complicatesthe physical interpretation of results, since, e.g., the in-stantaneous energy operator m ˆv / V is different from ˆ H , which prohibits an interpretation of ρ ii as a measur-able probability. These differences also explainwhy the matrix elements H ,ij in Eq. (14) for spatiallyconstant A and those obtained from Eq. (15) deviatefrom each other. While both versions of the interac-tion Hamiltonian lead to identical results for observablequantities as expected, it has been pointed out that theuse of approximations, such as the rotating wave approx-imation discussed in Section III E, can cause deviations between the two formulations.
In the following, wewill use the interaction operator of the form Eq. (15).
B. Optical Bloch Equations in Standard Form
From Eq. (4), we obtain with Eq. (15) in the dipoleapproximation the (multilevel) Bloch equations ∂ t ρ ij = − i ω ij ρ ij + i (cid:126) (cid:88) (cid:96) ( d i(cid:96) ρ (cid:96)j − d (cid:96)j ρ i(cid:96) ) E + (cid:88) mn D ijmn ρ mn , (16)with the transition frequencies ω ij = ( E i − E j ) / (cid:126) . Ifwe furthermore restrict the description of dissipative ef-fects to incoherent transitions and dephasing, Eqs. (9)and (12), Eq. (16) simplifies to ∂ t ρ ij = − i ω ij ρ ij + i (cid:126) (cid:88) (cid:96) ( d i(cid:96) ρ (cid:96)j − d (cid:96)j ρ i(cid:96) ) E − γ ij ρ ij , i (cid:54) = j, (17a) ∂ t ρ ii = i (cid:126) (cid:88) (cid:96) ( d i(cid:96) ρ (cid:96)i − d (cid:96)i ρ i(cid:96) ) E + (cid:88) j (cid:54) = i r j → i ρ jj − r i ρ ii . (17b)Although quantum optoelectronic devices can in princi-ple comprise a single isolated quantum system, for ex-ample a QD, in general they are based on ex-tended nanostructures such as quantum well structures,or an ensemble of many quantum systems such as multi-quantum-dot structures. This requires a position re-solved model, where the device is described by a represen-tative quantum system with density matrix ρ ij ( x , t ) ateach device position x . Furthermore, also the parameters ω ij , d ij , γ ij , r j → i and r i in Eqs. (16) and (17) generallydepend on x for inhomogeneous device structures, such as multi-section lasers. C. Optical Dipole Matrix Element
The Hamiltonian part of the Bloch equations, Eqs. (16)and (17), requires the dipole matrix element vectors d i(cid:96) ofthe optical transitions and the eigenenergies of the quan-tized states as an input. These can be computed frommodels derived from the stationary Schrödinger equation,such as the effective mass or k.p approach, as shortly dis-cussed in the following.In Fig. 3, the band structure of GaAs as an ex-emplary direct bandgap semiconductor material is dis-played. Shown is the conduction band (solid line) andthe valence band (dashed lines), consisting of heavy hole,light hole and split-off band. The holes tend to accu-mulate near the valence band maximum which is alwaysat the Γ point where the crystal wavevector is k = .1 FIG. 3. Band structure of gallium arsenide (GaAs), obtainedbased on a simple pseudo-potential tight-binding methodwithout spin-orbit coupling.
Shown is the valence band(dashed lines), the conduction band (solid line) and theparabolic dispersion relation assumed for the Γ valley in ef-fective mass approximation (dotted line). For direct bandgap semiconductors, the global conduc-tion band minimum where the electrons accumulate hap-pens to be also at the Γ point, and thus conservation ofcrystal momentum can be satisfied for radiative electron-hole recombination. This process is much less likely inindirect bandgap semiconductors, where the global con-duction band minimum is not at the Γ point and the pro-cess must additionally involve a phonon or crystal defectto achieve momentum conservation.Assuming a direct bandgap semiconductor, it is practi-cal to write the full wavefunction F i ( r ) of the initial andfinal state as a product of periodic Bloch function u v i ( r ) at the Γ point of band v i and an envelope wavefunction ϕ i ( r ) describing the slowly varying spatial modulation ofthe full wavefunction across the nanostructure. Whilea quantized state in a given band generally also containscontributions from neighboring bands, in a first approx-imation only the contribution of the dominant band isconsidered, F i ( r ) = u v i ( r ) ϕ i ( r ) . (18)In quantum well structures, the material compositionchanges only along the growth direction x . Here, quan-tum confinement only occurs in x direction, while thecarriers can move freely in the yz -plane. Thus, we canmake the ansatz ϕ i ( r ) = S − / ψ i ( x ) exp (i k y y + i k z z ) . (19)Here, S is the in-plane cross section area, k = [ k y , k z ] T denotes the in-plane wavevector in the yz -plane where T indicates the transpose, and ψ i ( x ) is the (gener-ally k dependent) one-dimensional envelope wavefunc-tion in confinement direction. In Fig. 4, the full wave-functions F i and corresponding envelope wavefunctions ϕ i are schematically illustrated for the two lowest con-duction band states and the valence band ground state Bandgap energyConduction bandValence band(heavy hole)
FIG. 4. Schematic illustration of full bound state wave-functions F i = u v i ϕ i (solid) and envelope wavefunctions ϕ i (dashed) in the conduction and valence band of a quantumwell. of a quantum well. Similar considerations apply to quan-tum wires, where quantum confinement occurs in twodimensions while the carriers can move freely along thethird coordinate. In QDs, the carriers are confined in allthree dimensions.
1. Computation of Envelope Wavefunction
Neglecting the coupling between conduction and va-lence bands, the simplest model for computing ϕ i ( r ) in aquantum structure is the Ben Daniel-Duke model, whichworks well for low-lying conduction band states in the Γ valley and also generally at the heavy hole valence bandmaximum. Here, we describe the dispersion relationbetween energy and wavevector around the Γ point by E ( k ) = V + (cid:126) | k | / (2 m ∗ ) which corresponds to a sec-ond order expansion, as illustrated by the dotted linein Fig. 3. The position dependent material compositionin nanostructures causes the effective mass m ∗ and bandedge energy V to depend on r , where V additionally con-tains the externally applied bias. Within this model, thestationary effective mass Schrödinger equation is givenby (cid:26) − (cid:126) (cid:20) ∇ m ∗ ( r ) ∇ (cid:21) + V ( r ) − E i (cid:27) ϕ i ( r ) , (20)where E i denotes the eigenenergy of state i . For thevalence band, commonly the hole picture is adopted toavoid a negative effective mass in Eq. (20). For the tran-sition between a conduction band electron state witheigenenergy E i and a valence band hole state with energy E j , the transition energy is then given by E i + E j + E g where E g denotes the bandgap energy, i.e., the energydifference between valence band maximum and conduc-tion band minimum. In quantum well systems, m ∗ and2 V only depend on the x coordinate, and Eq. (20) canbe reduced to the one-dimensional effective mass equa-tion by inserting Eq. (19). Similar considerations applyto quantum wires where m ∗ and V only depend on twocoordinates.The Ben Daniel-Duke model in Eq. (20) can be ex-tended, e.g., by accounting for band bending due to spacecharge effects in the potential, which is self-consistentlyincluded by solving Eq. (20) together with the Poissonequation. Furthermore, an energy dependent effec-tive mass can be introduced to include nonparabolicityeffects associated with the deviation of the dispersion re-lation from the parabolic form assumed above.
A further refined treatment of the conduction and va-lence bands, which accounts for band coupling, is usu-ally performed based on k.p theory, initially proposed byKane and Luttinger and Kohn.
Here the enve-lope wavefunctions are not scalar, but a multicomponentvector containing contributions from all the bands con-sidered. In many structures, strain arising from the lat-tice mismatch between the different semiconductor com-pounds plays an important role, and can be consideredbased on the Bir-Pikus model.
For modeling interbanddevices, eight-band k.p is a common option which consid-ers the top three valence bands and the lowest conductionband, along with spin orientation.
This approachis routinely applied to nanostructures, such as quantumdots, wires and wells.
If only valence band statesare considered, a restriction to six bands is possible.
This approach is sometimes also combined with Eq. (20)for the conduction band, assuming that it is decoupledfrom the valence bands. On the other hand, it has beenfound that for certain cases, eight-band k.p is not ac-curate enough. For example, a 14-band k.p approachwhich also includes the second conduction band in III-V semiconductors has been developed to obtain a moreaccurate conduction band dispersion relation at higherenergies, and 14-band k.p has also yielded improvedresults for SiGe/Si heterostructures.
2. Inter- and Intraband Dipole Matrix Elements
The dipole matrix element is best evaluated by com-puting the expectation value of the momentum operator ˆp = − i (cid:126) ∇ . Employing the product rule and exploitingthe fact that the periodic Bloch functions and envelopewavefunctions vary on two different length scales, we canwith Eq. (18) write (cid:104) F i | ˆp | F j (cid:105) ≈ (cid:104) ϕ i | ϕ j (cid:105) (cid:104) u v i | ˆp (cid:12)(cid:12) u v j (cid:11) + (cid:104) u v i (cid:12)(cid:12) u v j (cid:11) (cid:104) ϕ i | ˆp | ϕ j (cid:105) . (21)For transitions between conduction and valence bandstates, the first term dominates because the Bloch func-tions vary much more rapidly than the envelope wave-functions. Using (cid:104) F i | ˆp | F j (cid:105) = i m e E g (cid:104) F i | ˆr | F j (cid:105) / (cid:126) withthe electron mass m e and band gap energy E g , theinterband dipole matrix element can then in a first ap- proximation be written as (cid:104) F i | ˆr | F j (cid:105) = − i (cid:126) m − E − (cid:104) u v i | ˆp (cid:12)(cid:12) u v j (cid:11) (cid:104) ϕ i | ϕ j (cid:105) . (22)For intraband optical transitions, we have v i = v j , (cid:104) u v i | ˆp | u v i (cid:105) = 0 and (cid:104) u v i | u v i (cid:105) = 1 , and thus Eq. (21)yields (cid:104) F i | ˆp | F j (cid:105) ≈ (cid:104) ϕ i | ˆp | ϕ j (cid:105) and analogously (cid:104) F i | ˆd | F j (cid:105) ≈ (cid:104) ϕ i | ˆd | ϕ j (cid:105) . (23)In quantum well systems, confinement only occurs inthe growth direction x , and the envelope wavefunctionhas the form given by Eq. (19). For transitions betweena conduction band state | ψ i , k (cid:105) and a valence band state | ψ j , k (cid:48) (cid:105) , (cid:104) ϕ i | ϕ j (cid:105) = (cid:104) ψ i | ψ j (cid:105) δ k , k (cid:48) , i.e., the optical tran-sition is k conserving. The absolute value of the dipolematrix element can be approximately written as | e (cid:104) F i k | ˆr | F j k (cid:48) (cid:105)| = c ij E − P cv (cid:104) ψ i | ψ j (cid:105) δ k , k (cid:48) , (24)where e denotes the polarization direction of the elec-tric field, and P cv ≈ . .. × eV for most commonsemiconductors. For transitions between conductionband and heavy hole states, c ij = 2 − / for polariza-tion in in-plane direction and c ij = 0 for polarizationin growth direction. For transitions between conductionband and light hole states, c ij = 6 − / for polarizationin in-plane direction and c ij = 2 × − / for polarizationin growth direction. For intraband transitions occurringbetween quantized levels in the conduction band of quan-tum wells, as are for example employed for QCLs, theenvelope wavefunctions again assume the form Eq. (19).The dipole matrix element between an initial state | ψ i , k (cid:105) and a final state | ψ j , k (cid:48) (cid:105) is then with Eq. (23) given by d i k ,j k (cid:48) = d ij δ k , k (cid:48) , where d ij = (cid:104) ψ i | ˆd | ψ j (cid:105) = − e e x (cid:90) ψ ∗ i xψ j d x. (25)Here, e x denotes the unit vector in x direction, and onlythe dipole matrix element for polarization in growth di-rection x is nonzero. Notably, this is different from tran-sitions between conduction band and heavy hole statesin quantum wells where the x component of d ij is zero,as discussed above. In Fig. 5, the possible field polar-ization directions for interband [Fig. 5(a), (b)] and intra-band [Fig. 5(c)] transitions are indicated. For quantumwell lasers, Fig. 5(a), (b) and (c) correspond to the stan-dard edge-emitting, vertical-cavity surface-emitting andquantum cascade laser.In quantum dots, the uppermost valence band eigen-states usually exhibit heavy hole character. Thus,band coupling effects can often be neglected in Eq. (22)for interband transitions between the heavy-hole-likestates and low-lying conduction band states. Within theframework of these assumptions, only optical dipole tran-sitions between hole and electron states with equal quan-tum numbers are allowed, and the envelope wavefunctionoverlap in Eq. (22) typically approaches (cid:104) ϕ i | ϕ j (cid:105) ≈ for3 (a)(b)(c) FIG. 5. Polarizations with non-vanishing dipole matrix el-ements for (a) and (b) interband and (c) intraband transi-tions. The dashed arrow in (a) indicates that for this polar-ization, only transitions between conduction band and lighthole states contribute. the allowed transitions.
The symmetry of the wave-functions can however be affected by inhomogeneities inshape and composition of the quantum dots as well aspiezoelectric fields, resulting in additional weakly allowedtransitions.
Moreover, due to the strong confine-ment in quantum dots, Coulomb interactions tend toplay a pronounced role, causing energy shifts as wellas somewhat altered selection rules. Such effects canbe included in a more complete description based onthe electron-hole-pair picture, which replaces the single-carrier envelope wavefunctions ϕ i and ϕ j in Eq. (22) byexpressions for the excited electron-hole pair state andthe corresponding ground state. Intraband transi-tions, which are mainly relevant in the context of quan-tum dot infrared photodetectors, are again described byEq. (23).
D. Non-Redundant Density Matrix Representation
For a discrete-level system with N eigenstates | j (cid:105) , thedensity matrix contains N real diagonal elements and N − N complex off-diagonal elements which are relatedby ρ ij = ρ ∗ ji . Furthermore considering the trace condi-tion Tr { ˆ ρ } = 1 , the density matrix can be represented by N − non-redundant, real-valued elements, which areconveniently written as a vector S . This non-redundantrepresentation is for example achieved by the coherencevector (or pseudospin) representation, which has also been found useful for numerically efficient implementa-tions of the MB equations. For this purpose, thedensity matrix operator ˆ ρ is composed as ˆ ρ = 1 N ˆ I + 12 N − (cid:88) j =1 S j ˆ s j . (26)Here, ˆ s j are generators of the Lie algebra of SU( N )which are traceless Hermitian operators fulfilling thecondition Tr { ˆ s j ˆ s k } = 2 δ jk , and ˆ I is the identity op-erator. ˆ I and ˆ s j can be represented by correspond-ing N × N matrices. A possible choice for the gener-ators ˆ s = { ˆ u , ˆ u , . . . , ˆ v , . . . , ˆ w , . . . , ˆ w N − } consistsof N ( N − / generator pairs ˆ u jk = ˆ t jk + ˆ t kj , ˆ v jk = − i (cid:0) ˆ t jk − ˆ t kj (cid:1) , (27)and N − generators ˆ w l = − (cid:104) l − ( l + 1) − (cid:105) / (cid:32) l (cid:88) (cid:96) =1 ˆ t (cid:96)(cid:96) − l ˆ t l +1 ,l +1 (cid:33) , (28)where ˆ t jk = | j (cid:105) (cid:104) k | is the transition-projection operator,and the indices satisfy ≤ j < k ≤ N and ≤ l ≤ N − . For N = 2 and N = 3 these generators producethe Pauli and the Gell-Mann matrices, respectively.The elements of the coherence vector S are defined as S j = Tr { ˆ ρ ˆ s j } using the Hilbert-Schmidt inner product.Since both ˆ ρ and the generators ˆ s j are Hermitian, thevector elements are real. A similar transform can be ap-plied to the Lindblad equation. Inserting Eq. (26) intoEq. (3) and applying Tr {· ˆ s k } yields Tr { ∂ t ˆ ρ ˆ s k } = Tr N − (cid:88) j =1 ∂ t S j ˆ s j ˆ s k = 12 N − (cid:88) j =1 ∂ t S j Tr { ˆ s j ˆ s k } = ∂ t S k (29)for the left hand side. For the right hand side we canwrite Tr {L (ˆ ρ ) ˆ s k + D (ˆ ρ ) ˆ s k } = Tr {L (ˆ ρ ) ˆ s k } + Tr {D (ˆ ρ ) ˆ s k } , (30a) Tr {L (ˆ ρ ) ˆ s k } = Tr (cid:110) N − L ( ˆ I )ˆ s k (cid:111)(cid:124) (cid:123)(cid:122) (cid:125) =0 + N − (cid:88) j =1
12 Tr {L (ˆ s j ) ˆ s k } (cid:124) (cid:123)(cid:122) (cid:125) = L jk S j , (30b) Tr {D (ˆ ρ ) ˆ s k } = Tr (cid:110) N − D ( ˆ I )ˆ s k (cid:111)(cid:124) (cid:123)(cid:122) (cid:125) S eq j + N − (cid:88) j =1
12 Tr {D (ˆ s j ) ˆ s k } (cid:124) (cid:123)(cid:122) (cid:125) D jk S j , (30c)4since both superoperators L and D are linear. Notingthat L ( ˆ I ) = 0 and arranging Eqs. (29) and (30a)-(30c) inmatrix-vector form yields ∂ t S = ( L + D ) S + S eq , (31)where L and D are (cid:0) N − (cid:1) × (cid:0) N − (cid:1) real matricesand S eq denotes the equilibrium coherence vector.Alternatively, one can start from Eq. (4), where thesuperoperators L and D are represented as N × N ma-trices. This Liouville space representation was used forexample in, where column-major order was applied tomap the indices ( i, j ) (cid:55)→ k and ( m, n ) (cid:55)→ l . In this case,the density matrix is represented as N column vector R , and the Lindblad equation reads ∂ t R = (cid:18) − i (cid:126) L + D (cid:19) R , (32)where L = ˆ H T ⊗ ˆ I − ˆ I ⊗ ˆ H and D = (cid:88) k (cid:20) ˆ L T k ⊗ ˆ L k − (cid:16) ˆ L T k ˆ L ∗ k ⊗ ˆ I + ˆ I ⊗ ˆ L † k ˆ L k (cid:17)(cid:21) . (33)Here, ⊗ denotes the Kronecker product. Then, since theHilbert-Schmidt inner product reads Tr { ˆ a † ˆ b } = a † b = b † a in this representation, where the vectors a and b are the matrices ˆ a and ˆ b in column-major order, we canwrite the transform from Liouville space to the coherencevector representation as S j = Tr { ˆ ρ ˆ s j } = s † j R , S = T † R , (34)where the columns of the transformation matrix T arethe generators s j in column-major order. Conversely, thevector R can be recovered by R = 1 N I + 12 T S , (35)where I is the identity matrix in vectorized form. Usingthese transform relations, we can rewrite Eq. (32) as T ∂ t S = (cid:18) − i (cid:126) L + D (cid:19) T S + (cid:18) − i (cid:126) L + D (cid:19) N I , (36)and simplify the result by left-multiplication with T † to ∂ t S = (cid:18) − i2 (cid:126) T † L T + 12 T † D T (cid:19) S + 1 N T † D I , (37)where we used the orthogonality of the generators( T † T = I , where I is the N × N identity matrix)and the fact that the commutator of the identity is zero( L I = 0 ). This corresponds to Eq. (31).As we shall see in Section IV, the derivative of themacroscopic polarization ∂ t P q has to be calculated forthe Maxwell-Bloch equations. Naturally, it must be ex-pressed as function of the vector S . By replacing the trace operation and inserting the transformation rule, we canwrite for this term ∂ t P q = n Tr (cid:110) ∂ t ˆ ρ ˆd (cid:111) = n u † ∂ t R = 12 n u † T [( L + D ) S + S eq ] , (38)where u is the vectorized dipole moment operator and n denotes the carrier number density. Note that theelements of u could be vectors themselves, depending onwhether one ore more dimensions are considered.Using the dipole approximation, Eq. (15), we plug inthe Hamiltonian ˆ H = ˆ H − ˆdE ( t ) , which can be repre-sented with two matrices L and L in Liouville space(and two matrices L and L in coherence vector repre-sentation, respectively). Since u † L = − u † (cid:104) ˆd T E ( t ) ⊗ ˆ I − ˆ I ⊗ ˆdE ( t ) (cid:105) = − (cid:104) ˆd T E ( t ) ⊗ ˆ I u − ˆ I ⊗ ˆdE ( t ) u (cid:105) † = − (cid:104) vec (cid:16) ˆ I ˆdˆdE ( t ) (cid:17) − vec (cid:16) ˆdE ( t ) ˆd ˆ I (cid:17)(cid:105) † = 0 , (39)where vec denotes the vectorization of an operator, vec( ˆd ) = u , and the Hermitian property of the opera-tors involved as well as the properties of the Kroneckerproduct have been exploited, the polarization does notdepend on the electric field and Eq. (38) can be refinedas ∂ t P q = 12 n u † T [( L + D ) S + S eq ] . (40) E. Rotating Wave Approximation
The Bloch equations (17) are solvable only under spe-cial conditions, like | ∆ M | = 1 transitions in hydrogen-like atoms excited with circularly polarized light. Inparticular, closed analytical solutions do not exist forthe basic and very important case of excitation with amonochromatic, linearly polarized field.
Furthermore,the numerical solution of the Maxwell-Bloch equationsrequires high spatiotemporal resolution since the fieldsas well as the off-diagonal density matrix elements inEq. (17) oscillate with the optical period. For these rea-sons, the rotating wave approximation (RWA) is com-monly invoked, which significantly reduces the numericalburden and enables an analytical treatment of the Blochequations, at least for incident monochromatic radiationand some other relevant cases.The RWA is only applicable for not too broadbandoptical fields, which can then be separated into aslowly varying amplitude, given in complex notation by E ( x , t ) = | E ( x , t ) | exp [i φ ( x , t )] , and a rapidly oscillat-ing carrier with frequency ω c > . We note that there isno unique definition of ω c , but rather any choice which5ensures that all relevant spectral components are closeto ω c will suffice (for optical fields with symmetric powerspectra, it obviously makes sense to pick the center fre-quency). In complex notation, the electric field can thenbe written as E = 12 E exp ( − i ω c t ) + c.c., (41)where c.c. denotes the complex conjugate. Furthermoreassuming that all transitions between pairs of states i and j with non-negligible coupling to the optical fieldare in near-resonance, | ω ij | ≈ ω c , the corresponding off-diagonal density matrix elements are transformed into arotating reference frame, ρ ij = η ij exp [ − sgn ( ω ij ) i ω c t ] , (42)where sgn denotes the sign function. Inserting Eqs. (41)and (42) in Eq. (17), multiplying both sides of Eq. (17a)with exp [sgn ( ω ij ) i ω c t ] and applying the RWA, i.e., dis-carding all rapidly oscillating terms ∝ exp ( ± i ω c t ) and exp ( ± ω c t ) , we obtain ∂ t η ij = i∆ ij η ij + i2 (cid:126) ( ρ jj − ρ ii ) d ij (cid:26) EE ∗ (cid:27) − γ ij η ij , ω ij (cid:26) > < (cid:27) , (43a) ∂ t ρ ii = 1 (cid:126) (cid:88) jω ij > (cid:61) { d ji η ij E ∗ } + 1 (cid:126) (cid:88) jω ij < (cid:61) { d ji η ij E } + (cid:88) j (cid:54) = i r j → i ρ jj − r i ρ ii , (43b)with ∆ ij = sgn ( ω ij ) ( ω c − | ω ij | ) . As discussed above, theRWA is only applicable if the near-resonance condition isfulfilled, i.e., all significant spectral components E ( ω ) ofthe field are close to resonance with all relevant opticaltransitions at frequencies ω ij , | ω − | ω ij || (cid:28) | ω ij | . As asecond condition, the interaction energy must be so smallthat the eigenfrequencies of the quantum system are notconsiderably perturbed, i.e., | d ij E | / (cid:126) (cid:28) | ω ij | . IV. MAXWELL-BLOCH EQUATIONS
The optical field propagation in the device is classicallydescribed in terms of Maxwell’s equations. Assumingthat the magnetization is negligible at optical frequen-cies, we can write Faraday’s and Ampère’s law for theelectric field E and magnetic field H as ∇ × E = − µ ∂ t H , (44a) ∇ × H = (cid:15) (cid:15) r ∂ t E + σ E + J q = (cid:15) (cid:15) r ∂ t E + σ E + J f + ∂ t P q . (44b) E , H , and J q are functions of both t and x . J q = J f + ∂ t P q denotes the total current density contribution ZX FIG. 6. Schematic illustration of modeling based on Maxwell-Bloch equations. For an exemplary optical resonator, the dis-tribution of representative quantum systems is shown alongwith the cavity field. of the quantum systems. Here, J f and ∂ t P q correspondto the current density due to free carrier motion and thepolarization current density, respectively, where P q is themacroscopic polarization. In Fig. 6, the coupled modelingof the field propagation and the quantum system dynam-ics is schematically illustrated. For extended nanostruc-tures such as quantum well structures or ensembles ofQDs, the medium must be described by a representativequantum system at each position x . The Bloch equa-tions, Eq. (17), are coupled to Eq. (44) via E . On theother hand, Eq. (44) is coupled to Eq. (17) via J q . Forpractical reasons, we consider the background polariza-tion due to the host medium separately by the (generally x dependent) dielectric constant (cid:15) r ( x ) , where we have fornow neglected any frequency dependence, and assumedlinearity and isotropy of the host. Likewise, we includethe absorption of the host medium by a scalar conduc-tivity σ ( x ) , which gives rise to an ohmic current contri-bution σ E in Eq. (44b). Furthermore, (cid:15) and µ denotethe vacuum permittivity and permeability, respectively. A. Macroscopic Polarization And Current Density
Here, as above, we use the position variable r to re-solve microscopic behavior, while the variable x describesthe position in the modeled device or geometry andrefers to macroscopic dependencies, obtained from mi-croscopic models by adequate ensemble averaging. InMaxwell’s equations (44), the total macroscopic currentdensity contribution of the quantum systems is given by J q ( x , t ) = J f ( x , t ) + ∂ t P q ( x , t ) , where the free chargecurrent density J f and polarization current density ∂ t P q J f is for exam-ple induced by electrical pumping or generated by thephotovoltaic effect, while ∂ t P q is associated with thebound charge oscillations induced by the optical field.Microscopically, in nanostructured devices the carriers inbound or quasi-bound states may contribute to ∂ t P q bycoherent or incoherent interaction with the optical field,as well as to J f via coherent transport such as tunnel-ing and incoherent transport such as scattering-inducedhopping. Thus, it makes sense to treat polarization andfree current density together.The macroscopic polarization P q can be obtained fromthe dipole moment of the quantum system, given bythe expectation value of the dipole moment operator (cid:68) ˆd (cid:69) ( t ) = Tr (cid:110) ˆ ρ ( t ) ˆd (cid:111) . P q at a position x is then ob-tained by summing over the quantum systems in a vol-ume V p around x , P q = V − (cid:88) i (cid:68) ˆd i (cid:69) , (45)where V p is chosen big enough to obtain a smooth depen-dence of P q , but small enough so that spatial variationson classical length scales can still be resolved. For a largeensemble of identical systems with carrier number den-sity n , the polarization is then given by P q = n (cid:68) ˆd (cid:69) = n q Tr { ˆr ˆ ρ } , (46)where ˆ ρ is the density operator of a representative quan-tum system at position x . On the other hand, the electriccurrent in the quantum system can be computed fromthe expectation value of carrier velocity (cid:104) ˆv (cid:105) in the sys-tem, where the velocity operator is defined in the Heisen-berg picture by the time derivative of the position op-erator ˆr H , ˆv H = d t ˆr H . For the coherent contributioncorresponding to the Hamiltonian part in Eq. (3), wethen obtain with the Ehrenfest equation (cid:104) ˆv (cid:105) = d t (cid:104) ˆr (cid:105) , where we have dropped the index H since expectation val-ues for physical observables are independent of the cho-sen picture. In the Schrödinger picture, we thus obtain (cid:104) ˆv (cid:105) = d t Tr { ˆr ˆ ρ } = Tr { ˆr d t ˆ ρ } , which is also valid for theincoherent contribution induced by the Lindblad opera-tor term in Eq. (3). Thus, I = q |(cid:104) ˆv (cid:105)| /L corresponds tothe current through an individual (single-carrier) quan-tum system, where L here indicates the system length inthe direction of current flow, and L/ |(cid:104) ˆv (cid:105)| is the transittime of the carrier through the system. Again averag-ing over a large ensemble of identical systems, we obtainthe macroscopic current density J f = n q Tr { ˆr d t ˆ ρ } . Wenote that this result is the same as for the polarizationcurrent density, obtained by taking the time derivativeof Eq. (46), which reflects the fact that the carriers of thequantum system are responsible for both the free chargecurrent and polarization current. Even more, from a mi-croscopic standpoint, this distinction is inappropriate forour case. Thus we can write J q = n q Tr { ˆr ∂ t ˆ ρ } , (47) where ˆ ρ again describes a representative quantum systemat position x , and ∂ t ˆ ρ , given by Eq. (3), contains boththe coherent and incoherent dynamics. Assuming an N -level system with orthonormal basis states | i (cid:105) , Eq. (47)can with the dipole matrix element d ij = q (cid:104) i | ˆr | j (cid:105) bewritten as J q = n (cid:88) i,j d ji ∂ t ρ ij , (48)with ∂ t ρ ij given by Eq. (4) or (17).A widely used criterion to distinguish between themacroscopic free charge and polarization current contri-butions in Eq. (48) is the frequency range, where com-monly ∂ t P q is expected to contain frequencies in therange of the driving optical field spectrum, while J f cov-ers the low-frequency and direct current contributions.In this context, we point out that due to nonlinear op-tical mixing, the polarization generally contains up- anddown-converted components. This especially ap-plies to nanostructured optoelectronic devices where gi-ant optical nonlinearities can be artificially engineered,and are actively exploited in both the optical and tera-hertz regime.
On the other hand, the electric cur-rent can contain components up to tens of GHz due toexternal modulation or back-coupling of the optical dy-namics to the electrical circuitry. Notably, in QCLs em-bedded into a micro-strip line, strong coupling of the co-propagating microwave current modulation and opticalwaveform has recently been found, indicating thata clear differentiation between free and polarization cur-rent contributions is not always possible. However, aspointed out above, such a distinction is also not neces-sary since the current density and polarization appear as J q = J f + ∂ t P q in Maxwell’s equations. Ultimately, thefrequency range of the measured electrical current willbe limited by both the measurement setup itself and theelectrical properties of the device, such as its intrinsiccapacitance.
1. Coherent Contribution
Using the Ehrenfest equation d t (cid:104) ˆr (cid:105) = i (cid:126) − (cid:68)(cid:104) ˆ H, ˆr (cid:105)(cid:69) , we can write the coherent part of the current density as J coh = i n q (cid:126) − (cid:68)(cid:104) ˆ H, ˆr (cid:105)(cid:69) . (49)In the following, we assume an effective mass Hamilto-nian of the form ˆ H = (1 / ˆp [ m ∗ ( ˆr )] − ˆp + V ( ˆr , t ) asused in Eq. (20), yielding J coh = (1 / n q (cid:68) [ m ∗ ( ˆr )] − ˆp + ˆp [ m ∗ ( ˆr )] − (cid:69) . (50)Using an orthonormal basis as for Eq. (48), and insertingthe unit operator ˆ I = (cid:82) d r | r (cid:105) (cid:104) r | in Eq. (50), we can ex-7press the result in terms of wavefunctions ϕ i ( r ) = (cid:104) r | i (cid:105) , J coh = (cid:126) n q (cid:88) i,j (cid:90) (cid:61) (cid:110) ρ ji ϕ ∗ i ( r ) [ m ∗ ( r )] − ∇ ϕ j ( r ) (cid:111) d r. (51)Equation (51) can also be interpreted as the current den-sity contribution of a representative individual (single-carrier) quantum system at the corresponding position,averaged over the associated volume V p = n − , wherethe microscopically resolved current density is given bythe familiar expression J coh = (cid:126) q (cid:88) i,j (cid:61) (cid:110) ρ ji ϕ ∗ i ( r ) [ m ∗ ( r )] − ∇ ϕ j ( r ) (cid:111) . (52)
2. Incoherent Contribution
The incoherent contribution to the current density isgiven by J inc = n q [ (cid:104) ˆv (cid:105) ] inc = n q Tr { ˆr [d t ˆ ρ ] inc } , which yields with Eqs. (3) and (4) J inc = n q Tr (cid:40) ˆr (cid:88) k (cid:18) ˆ L k ˆ ρ ˆ L † k −
12 ˆ L † k ˆ L k ˆ ρ −
12 ˆ ρ ˆ L † k ˆ L k (cid:19)(cid:41) = n (cid:88) ij d ji (cid:88) mn D ijmn ρ mn . (53)For an incoherent transition from a state α to β (cid:54) = α , weobtain with the corresponding Lindblad operator givenin Eq. (7) J α → β inc = r α → β n ( d ββ − d αα ) ρ αα − (cid:88) i (cid:54) = α (cid:60) { d iα ρ αi } . (54)Furthermore, inserting Eq. (12) in Eq. (53) yields thepure dephasing contribution between two levels α and β J αβ inc = − γ (cid:48) αβ n (cid:60) { d βα ρ aβ } , (55)with the pure dephasing rate γ (cid:48) αβ . The current contri-butions from incoherent transitions due to Eqs. (54) and(55) can also be rearranged so that J αβ hop = ( r α → β ρ αα − r β → α ρ ββ ) n ( d ββ − d αα ) (56)is the net current due to the hopping transport betweenstates α and β which corresponds to the classical rateequation description, and J αβ dep = − γ αβ n (cid:60) { d βα ρ aβ } (57)is the dephasing contribution due to the decay of thecorresponding off-diagonal matrix elements ρ aβ and ρ βa .Here, γ αβ = ( r α + r β ) / γ (cid:48) αβ is the total dephasingrate, including lifetime broadening and pure dephasing, and r α,β is given by Eq. (10). The total incoherent cur-rent density, resulting from incoherent transitions andpure dephasing, is then obtained by summing over alltransitions. With Eqs. (56) and (57), we obtain J inc = N − (cid:88) α =1 N (cid:88) β = α +1 (cid:16) J αβ hop + J αβ dep (cid:17) = n (cid:88) α (cid:88) β (cid:54) = α (cid:2) r α → β ρ αα ( d ββ − d αα ) − γ αβ (cid:60) { d βα ρ aβ } (cid:3) . (58) B. Slowly Varying Amplitude Approximation
Although the Bloch equations in RWA, Eq. (43), aresometimes solved in combination with the full Maxwell’sequations, Eq. (44), typically the RWA is combinedwith an envelope propagation equation, derived fromMaxwell’s equations under the assumption of a slowlyvarying field amplitude. In this way, above mentionedadvantages of the RWA, namely a significantly reducednumerical burden and a larger number of analytical so-lutions, also applies to the coupled Maxwell-Bloch sys-tem. Taking the curl of Eq. (44a) and eliminating H us-ing Eq. (44b) yields ∇ × ∇ × E = − n (1 − n ) c ∂ t E − µ σ∂ t E − µ ∂ t J f − µ ∂ t P q . (59)Here, c = ( µ (cid:15) ) − / is the vacuum speed of light.Furthermore, the background permittivity of the hostmaterial (cid:15) r ( x ) is here modeled as (cid:15) r = n (1 − n ) ,where ∆ n ( x, y ) (with the minimum value ) describesa transverse refractive index profile, as widely employedin waveguiding structures. For no free space charges, Gauss’s law dictates that ∇ D = 0 where D = P q + (cid:15) n (1 − n ) E is thedisplacement field in Eq. (59). Assuming an isotropicmedium, we can thus set ∇ ( ∇ E ) ≈ in the case of weaknonlinearity and weak inhomogeneity, or gener-ally if the field intensity transverse to the propagationdirection is slowly varying over an optical wavelength. This assumption is only fulfilled for weak waveguid-ing, i.e., if the relative changes of the refractive index | ∆ n b | /n b and its gradient | ∆ ( ∇ n b ) | / |∇ n b | over the dis-tance of a wavelength in the medium is small againstunity, where n b = n (1 − n ) / in Eq. (59). Fur-thermore, also the polarization contribution P q of thequantum structure must be compatible with the assump-tion of weak inhomogeneity. As discussed in SectionIII C, quantum structures, as modeled by the Bloch equa-tions, can be highly anisotropic; e.g., the dipole momentelement vector d ij of inter-conduction band transitions inquantum wells only has a nonzero component in growthdirection. If the optical field is however also polarized in8this direction, which is for example often the case in laserssince only the corresponding field component gets ampli-fied, then ∇ E ≈ can still hold for weak nonlinearity andinhomogeneity. Using ∇× ( ∇ × E ) = ∇ ( ∇ E ) −∇ E andsubsequently neglecting the term ∇ ( ∇ E ) , we obtain thegeneralized inhomogeneous wave equation ∇ E = n (1 − n ) c ∂ t E + µ σ∂ t E + µ ∂ t J f + µ ∂ t P q . (60)For deriving the slowly varying amplitude approxima-tion (SVAA), E and P q are written as a product of itsenvelope and carrier, as done above for the derivation ofthe RWA. However, in contrast to Eq. (41), we also takeinto account the spatial dependence of the carrier, wherewe assume that the direction of the optical energy flowat every position is close to a reference direction definedby the carrier wavevector k c , which corresponds to theparaxial approximation. This assumption is for exampletypically fulfilled in laser resonators or optical fibers. In-troducing the complex-valued field and polarization am-plitudes, E ( x , t ) and P ( x , t ) , and assuming propagationalong the z direction, we have E ( x , t ) = 12 E ( x , t ) exp (i k c z − i ω c t ) + c.c., (61a) P q ( x , t ) = 12 P ( x , t ) exp (i k c z − i ω c t ) + c.c., (61b)with | k c | = n ω c /c . We note that although Eq. (61a)contains the term exp (i k c z ) not included in Eq. (41),the Bloch equations in RWA, Eq. (43), remain unchangedsince exp (i k c z ) cancels out. To apply the SVAA, we in-sert Eq. (61) in Eq. (60). Just as for the RWA, we assumethat all significant spectral components of the field areclose to ω c , i.e., at frequencies ω c + ∆ ω with | ∆ ω | (cid:28) ω c . This implies that ∂ t E can be neglected against − ω c ∂ t E , as can be seen in Fourier domain wherethe two terms become − ∆ ω E (∆ ω ) and − ω c ∆ ω E (∆ ω ) .Similarly, also σ∂ t E and ∆ n ∂ t E can be dropped against − i ω c σ E and − i ω c ∆ n E . The polarization amplitude P ,introduced in Eq. (61b), couples the optical propaga-tion equation to the Bloch equations, Eq. (43), as fur-ther discussed in Section IV B 1. The RWA impliesthat also P is narrowband, which means that for ex-ample harmonic or difference frequency generation can-not be included. Thus, similarly as for the field, ∂ t P and − ω c ∂ t P can be neglected against − ω P . In addi-tion, the paraxial approximation implies that ∂ z E can beneglected against i k c ∂ z E . Finally multiplying all termswith exp (i ω c t − i k c z ) and discarding all rapidly oscillat-ing terms, which also eliminates J f since it is assumedto contain only low frequency components (see SectionIV A), we arrive at ∂ t E ± cn ∂ z E = − i ω c ∆ n E + 12 n (cid:18) i c ω c ∇ E + i ω c (cid:15) P − σ(cid:15) E (cid:19) . (62) Here, ∇ = ∂ x + ∂ y denotes the transverse Laplace op-erator. The ”+” and ”-” signs in Eq. (62) are for forwardand backward propagation corresponding to k c > and k c < , respectively. For counterpropagating fields whichfor example arise in Fabry-Pérot resonators, the standingwave pattern causes a position dependent inversion grat-ing, also referred to as spatial hole burning. This effectis not yet included in Eq. (62), and its implementation isdiscussed in Section V C 2.
1. Polarization in Rotating Wave Approximation
In the RWA, the off-diagonal density matrix elements ρ ij that are associated with near-resonant optical transi-tions are represented in terms of transformed elements η ij in a rotating reference frame, as obtained with Eq. (42).Writing the total current as J q = J f + ∂ t P q as in Eq. (47),and assigning the low-frequency contributions to J f andthe optical contributions to P q , we see from Eq. (48) thatthe transformation into the rotating frame only affectsthe evaluation of the polarization P q . With Eq. (46) andEq. (42), we obtain P q = n (cid:88) ω ij > d ji η ij exp ( − i ω c t ) + c.c. + n (cid:88) i d ii ρ ii . (63)For inclusion of optical propagation, the RWA is often notcoupled to the full Maxwell equations, but rather solvedtogether with Eq. (62) in SVAA which contains the polar-ization in terms of the amplitude P . As discussed above,we have to replace exp ( − i ω c t ) by exp (i k c z − i ω c t ) inEq. (63) since the SVAA also takes into account the spa-tial dependence of the carrier. Comparing the resultingequation with Eq. (61b), and neglecting the quasi-staticdipole moment contribution (cid:80) i d ii ρ ii which does not os-cillate at the optical excitation frequency and thus dropsout in the SVAA, we obtain P = 2 n (cid:88) ω ij > d ji η ij . (64) C. Initial Conditions
The Bloch equations without or with RWA, Eq. (17)or Eq. (43), have to be supplemented by correspond-ing initial conditions at time t = t . Apart from spe-cial cases where the quantum system may be coherentlyprepared in a certain initial state such as a coherentsuperposition, the system will be initially in equilib-rium. The corresponding density matrix elements arethen obtained by setting ∂ t = 0 in Eq. (17) or (43) andassuming a vanishing optical field for t ≤ t , whichgives rise to a mixed state with off-diagonal elements ρ ij ( t = t ) = 0 and η ij ( t = t ) = 0 , respectively. Thediagonal elements ρ ii ( t = t ) = ρ eq ii are given by the9equilibrium occupation probabilities ρ eq ii , which can beobtained by setting d t = 0 in Eq. (9). This yields for asystem with N levels the linear equation system (cid:88) j (cid:54) = i r j → i ρ eq jj − r i ρ eq ii , i = 1 .. ( N − , (65) N (cid:88) i =1 ρ eq ii , (66)where r i is given by Eq. (10). The ρ eq ii do not necessarilycorrespond to a thermal distribution, but are rather de-termined by the transition rates which may for exampleinclude the pumping process in lasers. For inhomoge-neous device structures, the rates r i and r j → i generallydepend on position x , giving rise to x dependent ρ eq ii .Suitable initial conditions have also to be definedfor the Maxwell equations, Eq. (44), or the propagationequations in SVAA derived thereof, Eq. (62). Here wecannot choose identically vanishing fields, since the op-tical field would then remain zero throughout the MBsimulation. Laser seeding by spontaneous emission is of-ten mimicked by initializing the electric field with whiteGaussian amplitude noise, which can also be added atevery time step of the simulation to model spontaneousemission noise. In the SVAA, Eq. (62), the electric fieldis represented by its complex envelope function, and thuscomplex white Gaussian noise is used in this case.
D. Two-Level Approximation
In most cases, the simplest model with only two rele-vant states (e.g., an upper laser level and lower laserlevel ) is considered. The corresponding density matrixcontains the elements ρ , ρ and ρ = ρ ∗ . Assuming aclosed system, we obtain ρ + ρ = 1 . The dissipationin the Bloch equations, Eq. (17), is then parametrizedby the three rates γ , r → and r → , where γ = γ as described in Section II B 2, and Eq. (10) gives r = r → , r = r → . Introducing the population inversion w = ρ − ρ , we can substitute ρ = (1 − w ) / and ρ = (1 + w ) / in Eq. (17). Furthermore neglectingstatic dipole moments, d − d ≈ , we obtain ∂ t ρ = − i ω ρ − i w Ω − γ ρ , (67a) ∂ t w = 2i ( ρ ∗ Ω − ρ Ω ∗ ) − γ w + Γ . (67b)Here, Ω = (cid:126) − d E denotes the instantaneous Rabi fre-quency, γ = r → + r → and γ = γ are the pop-ulation inversion relaxation and dephasing rates, and Γ = r → − r → represents the (net) pumping ratefrom the lower to the upper level. Equation (67b) is of-ten also written as ∂ t w = 2i ( ρ ∗ Ω − ρ Ω ∗ ) − γ ( w − w eq ) , (68)where w eq = Γ /γ denotes the equilibrium populationinversion for Ω = 0 . t/T-6 -4 -2 0 2 4 6-1-0.500.51 uvw w (a)(b) FIG. 7. (a) Bloch vector components and (b) Bloch vectortrajectory for a dissipationless two-level system excited witha sech pulse.
A real-valued, redundance-free representation can beobtained by applying Eqs. (27) and (28), yielding thethree real-valued quantities u := (cid:104) ˆ u (cid:105) = ρ + ρ =2 (cid:60) { ρ } , v := − (cid:104) ˆ v (cid:105) = − i ( ρ − ρ ) = − (cid:61) { ρ } , and w := (cid:104) ˆ w (cid:105) = ρ − ρ . These are usually represented interms of the so-called Bloch vector S = [ u, v, w ] T , where aminus sign has been added to the definition of v in orderto obtain the usual convention for the Bloch vector. Separating Ω in its real and imaginary part Ω r + iΩ i ,Eq. (67) then becomes ∂ t u = − ω v + 2Ω i w − γ u, (69a) ∂ t v = ω u + 2Ω r w − γ v, (69b) ∂ t w = − i u − r v − γ w + Γ , (69c)which can also be written as ∂ t S = − r i ω × S − γ uγ vγ ( w − w eq ) . (70)The polarization term in Eq. (44) is then with Eq. (48)obtained as ∂ t P q = 2 n (cid:60) { d ∂ t ρ } = n ( (cid:60) { d } ∂ t u + (cid:61) { d } ∂ t v ) . (71)The time evolution of the Bloch vector S ( t ) can bevisualized in the Bloch sphere representation, where the0Bloch vector trajectory is displayed in a Cartesian coordi-nate system with axes u , v and w . For γ = γ = 0 , | S ( t ) | is conserved over time, as can be seen by mul-tiplying Eqs. (69a), (69b) and (69c) with u , v and w ,respectively, and adding the resulting equations, whichyields ∂ t (cid:0) u + v + w (cid:1) = 0 . For pure states, | S ( t ) | = 1 ,i.e., the tip of the Bloch vector moves along the surfaceof a unit sphere, the so-called Bloch sphere. For mixedstates, the tip is located within the Bloch sphere, cor-responding to | S ( t ) | < . In Fig. 7, the time evolutionof the Bloch vector components and the correspondingBloch vector trajectory are shown for a two-level systemwith γ = γ = 0 and initial conditions u = v = 0 , w = − . The optical field is assumed to be a sechpulse Ω = 2 T − sech ( t/T ) , which corresponds to a self-induced transparency soliton as further discussed in Sec-tion VI A 2, with T = 10 /ω chosen for this example.A further representation of the Bloch equations is ob-tained by assuming a real d and thus Ω = Ω r . SolvingEq. (69a) for v and using the result to eliminate v inEqs. (69b) and (69c) yields with Eq. (71) (cid:2) ∂ t + 2 γ ∂ t + (cid:0) ω + γ (cid:1)(cid:3) P q = − ω d (cid:126) n w d | d | E , (72a) ∂ t w = 2 ∂ t P q + γ P q (cid:126) n ω E − γ w + Γ . (72b)This representation can be seen as an extension ofthe classical Lorentz model for resonant polarization indielectrics, assuming the same mathematical form asEq. (72a) if we set w constant. Accordingly, Eq. (72)is mainly used in computational electrodynamics, es-pecially in combination with the finite-difference time-domain method, as a substitute for more basic classicalpolarization models.
1. Rotating Wave/Slowly Varying Amplitude Approximation
In the RWA, we obtain from Eq. (43) with
Ω = (cid:126) − d E ∂ t η = i∆ η −
12 i w Ω − γ η , (73a) ∂ t w = i ( η ∗ Ω − η Ω ∗ ) − γ ( w − w eq ) , (73b)where ∆ = ω c − ω denotes the detuning of the opti-cal field from the resonance frequency ω . In analogyto above, we can introduce the Bloch vector s for theoff-diagonal density matrix elements in RWA, with com-ponents s = η + η = 2 (cid:60) { η } , s = − i ( η − η ) = − (cid:61) { η } , s = w , and obtain in analogy to Eq. (69)with Ω = Ω r + iΩ i ∂ t s = ∆ s + Ω i w − γ s , (74a) ∂ t s = − ∆ s + Ω r w − γ s , (74b) ∂ t w = − Ω i s − Ω r s − γ ( w − w eq ) . (74c) The polarization term in the SVAA propagation equa-tion, Eq. (62), is then with Eq. (64) obtained as P = 2 n d η = n d ( s − i s ) . (75) V. REDUCTION TO ONE-DIMENSIONAL MODEL
Although the MB equations are sometimes solved intwo or even three spatial dimensions, the model isfrequently reduced to a single spatial coordinate in or-der to minimize the numerical load. This is usuallyachieved by assuming plane wave propagation in theMaxwell equations, Eq. (44), or the corresponding prop-agation equations in SVAA, Eq. (62).
For extendedbeams propagating in a homogeneous medium such asa gas or bulk solid-state medium, the plane wave ap-proximation may be a reasonable assumption. For op-toelectronic devices which are the focus of this paper,the light is usually strongly guided, often with sub-wavelength confinement in at least one dimension. Here,the plane wave approximation is clearly too simplistic.However, optoelectronic devices such as semiconductor-based lasers often employ waveguiding structures whichare invariant in propagation direction z , in particularschemes with a suitable transverse refractive index pro-file or metal cladding. Such geometries provide lat-eral field confinement and give rise to guided mode so-lutions, i.e., field solutions which are at a given fre-quency ω characterized by a propagation constant and a z independent transverse field distribution. While someone-dimensional plane wave treatments have included alltransverse field components to describe elliptically or cir-cularly polarized light, we assume linearlypolarized waveguide modes in the following, and thusconsider a single transverse component of the electric andmagnetic fields. In Fig. 8, an exemplary waveguide struc-ture is schematically illustrated. A. Full Maxwell Equations
We employ the full Maxwell equations, Eq. (44), cou-pled to the Bloch equations, Eq. (16) or (17), to describethe carrier-light interaction and optical propagation ina waveguide geometry which is invariant with respectto the propagation direction z . Our goal is to extracta one-dimensional MB model with a single electric andmagnetic field component, as typically used in simula-tions due to the associated computational burden. Wefocus on guided mode solutions, which are at a given fre-quency ω characterized by a (generally complex) prop-agation constant β , and z invariant transverse field de-pendencies E t x,y,z ( x, y ) and H t x,y,z ( x, y ) for the electricand magnetic field components. Thus, we can make the1 x yz xL h FIG. 8. Schematic illustration of one-dimensional propaga-tion model for a waveguide structure. (a) The white rectangleon the front facet denotes the cross section of the quantumstructure, and the optical intensity distribution is indicated.(b) The refractive index profile is indicated, with darker colorscorresponding to higher refractive indices. ansatz E p ( x , t ) = (cid:60) (cid:8) E t p ( x, y ) exp (cid:0) i βz − i ωt (cid:1)(cid:9) = (cid:60) (cid:8) E ωp ( x ) exp ( − i ωt ) (cid:9) , (76a) H p ( x , t ) = (cid:60) (cid:8) H t p ( x, y ) exp (cid:0) i βz − i ωt (cid:1)(cid:9) = (cid:60) (cid:8) H ωp ( x ) exp ( − i ωt ) (cid:9) , (76b)with p = x, y, z . By inserting Eq. (76) into Eq. (44), thecomputation of the transverse mode profile in the xy -plane can be decoupled from the z coordinate and re-duces to a two-dimensional problem. For example, byeliminating the electric field, we obtain (cid:0) ∂ p + ∂ q (cid:1) H t p + ∂ q (cid:15) r (cid:15) r (cid:0) ∂ p H t q − ∂ q H t p (cid:1) = (cid:18) β − ω c (cid:15) r (cid:19) H t p (77)with p = x, y and q = y, x . Waveguiding is for instanceobtained by surrounding the optically active region withanother dielectric material featuring a lower refractiveindex as illustrated in Fig. 8(b), or with a metal cladding.Both cases can be described by a transversely dependentcomplex background permittivity (cid:15) r = (cid:15) r + i σ/ ( ω(cid:15) ) , (78)where (cid:15) r and σ generally depend on x , y and ω , and σ ac-counts for the conductivity or dielectric losses. Together with the boundary condition H t p,q → for x + y →∞ , Eq. (77) constitutes a complex eigenvalue problem.Equation (77) can for example be solved with the filmmode matching method, which is especially suitable forwaveguides with a rectangular cross section. The po-larization contribution P q in Eq. (44) due to the quan-tum systems is not yet considered in Eq. (77) since itis assumed to be small enough to be included in first-order perturbation theory, with negligible influence onthe transverse field distribution. Using ∇ H = 0 , we cancalculate the longitudinal component H t z ( x, y ) from H t x and H t y as H t z = i β − (cid:0) ∂ x H t x + ∂ y H t y (cid:1) . (79)Furthermore, the electric field components are withEq. (44b) obtained as ω(cid:15) (cid:15) r E t x = i ∂ y H t z + βH t y , (80a) ω(cid:15) (cid:15) r E t y = − βH t x − i ∂ x H t z , (80b) ω(cid:15) (cid:15) r E t z = i ∂ x H t y − i ∂ y H t x . (80c)For general solutions of Eq. (77), the polarizationvaries over the waveguide cross section. As indicated inFig. 8(b), in many optoelectronic devices, such as typi-cal standard edge-emitting and quantum cascade lasers,rectangular waveguides are used where the width in lat-eral y direction significantly exceeds its thickness in x di-rection. This allows an approximate treatment as a slabwaveguide structure, which is assumed to be infinitelyextended in y direction and thus can, to first order, bedescribed by (cid:15) r ( x ) . The field components are then as-sumed to be constant in y direction, which correspondsto setting ∂ y = 0 . The guided field solutions can bedivided into two classes: Transverse electric (TE) modesare characterized by E ωz = 0 , where for ∂ y = 0 all compo-nents except E ωy , H ωx and H ωz vanish as can be seen fromEqs. (79) and (80); similarly, transverse magnetic (TM)modes, characterized by H ωz = 0 , have for ∂ y = 0 onlynon-vanishing H ωy , E ωx and E ωz components. The y dependence of the field distribution may then be reintro-duced using the effective refractive index approximationmethod, which preserves the TE or TM character ofthe solution. From the discussion of dipole matrix el-ements for quantum well structures in Section III C, itfollows that standard edge-emitting lasers, which utilizeinterband transitions, preferably operate in TE mode [seealso Fig. 5(a)]. On the other hand, QCLs, which rely onintraband transitions, only operate in TM mode [see alsoFig. 5(c)].As pointed out above, simulations typically employ aplane-wave-type propagation model which only dependson the propagation coordinate z and time t , and considersa single transverse electric and transverse magnetic fieldcomponent. Our goal is to derive such equations, witha form equivalent to the Maxwell equations, for guidedrather than plane-wave propagation, as applies to manyphotonic devices and systems.2
1. Transverse Electric Mode
For TE modes in slab waveguides, Eq. (77) yields with H t y = 0 and ∂ y = 0 ∂ x H t x = (cid:18) β − ω c (cid:15) r (cid:19) H t x , (81)and the boundary conditions are given by H t x ( x → ±∞ ) → . From Eq. (44), we furthermoreobtain − ∂ z E y = − µ ∂ t H x , (82a) ∂ z H x − ∂ x H z = (cid:15) (cid:15) r ∂ t E y + σE y . (82b)The polarization contribution of the quantum systemsis not contained in Eq. (82b) since it will subsequentlybe included in a perturbative manner. Equation (82b)does not yet have the desired form since it contains an x derivative and the longitudinal field component in theterm ∂ x H z . With Eqs. (76b), (79) and (81), we obtain ∂ x H ωz = i β − (cid:18) β − ω c (cid:15) r (cid:19) H ωx , (83)where (cid:15) r = (cid:15) r + i σ/ ( ω(cid:15) ) . In the following, it is practicalto switch to the frequency domain, where Eq. (82a) iswith Eq. (76) given by ∂ z E ωy = i βE ωy = − i ωµ H ωx . (84)From Eq. (84), we see that the electric and magneticfields have the same transverse distribution. InsertingEq. (83) into Eq. (82b) in frequency domain, and employ-ing Eq. (84), we arrive at ∂ z H ωx = − i ω − µ − β E ωy . (85)In the following, the polarization contribution of thequantum system will be included as a perturbation. Re-deriving Eq. (81) from Maxwell’s equations Eq. (44),but now with the polarization contribution due to thequantum systems included, we see that the perturbationgenerated by the polarization on the E ωy field compo-nent is formally equivalent to an additional backgroundpermittivity ∆ (cid:15) r = (cid:15) − P ωy /E ωy , where P ωy contains thepolarization contribution of the quantum systems in fre-quency domain. In the following, we assume that thedevice operates in a single transverse mode with themagnetic field distribution H t x , possibly the fundamentalmode. Using the similarity of Eq. (81) to the Schrödingerequation in quantum mechanics, we can apply perturba-tion theory in an analogous matter. To first order, H t x remains unchanged, and for the eigenvalue β we obtainthe correction ∆ β = ω (cid:15) E ωy c (cid:82)(cid:82) ∞−∞ (cid:12)(cid:12) H t x (cid:12)(cid:12) P ωy d x d y (cid:82)(cid:82) ∞−∞ (cid:12)(cid:12) H t x (cid:12)(cid:12) d x d y . (86) For completeness, we also include integration over the y coordinate in Eq. (86) since the y dependence of H t x maybe reintroduced based on above mentioned effective re-fractive index method. It should be mentioned that if (cid:15) r in Eq. (81) has a non-vanishing imaginary part, theeigenvalue problem is non-Hermitian and strictly speak-ing, a biorthogonal basis set must be used. In thiscase, Eq. (86) serves as an approximation to the ex-act perturbation term. Furthermore, it is practical tosplit the unperturbed propagation constant β into a realand an imaginary part, β = β + i β (cid:48) . Here, β (cid:48) is re-lated to the power loss coefficient a by β (cid:48) = sgn ( β ) a/ ,with the sign function sgn . Assuming | β | (cid:29) | β (cid:48) | , wecan write β ≈ β + i | β | a + ∆ β . Introducing theeffective waveguide refractive index n eff ( ω ) defined by β = sgn ( β ) ωn eff /c , we then obtain ∂ z H ωx = − i ω(cid:15) n E ωy + (cid:15) cn eff aE ωy − i ω (cid:82)(cid:82) ∞−∞ (cid:12)(cid:12) H t x (cid:12)(cid:12) P ωy d x d y (cid:82)(cid:82) ∞−∞ (cid:12)(cid:12) H t x (cid:12)(cid:12) d x d y . (87) a. Field Confinement Factor Equations (84) and(87) effectively reduce the complexity of the propaga-tion problem from three spatial dimensions to a singlecoordinate z . However, for computing the integral inthe polarization term of Eq. (87), P ωy must be obtainedby solving the Bloch equations in the whole device vol-ume, using the full spatial field dependence given byEq. (76). This greatly impedes the numerical efficiencyof the one-dimensional propagation model. As indicatedin Fig. 8(a), frequently the transverse field distributiondoes not vary significantly across the quantum nanos-tructure, e.g., because the nanostructure covers only partof the waveguide cross section, preferably at the positionof maximum intensity. Consequently, also P ωy is approx-imately constant over the quantum system cross sectionand can be taken out of the integral in Eq. (86), whichcan then be written as Γ P ωy . Here, Γ denotes the fieldconfinement factor, which gives the overlap of the quan-tum nanostructure with the mode profile and is thus alsoreferred to as overlap factor. With Eq. (84), Γ can bewritten as Γ = (cid:82)(cid:82) A q (cid:12)(cid:12) H t x (cid:12)(cid:12) d x d y (cid:82)(cid:82) ∞−∞ (cid:12)(cid:12) H t x (cid:12)(cid:12) d x d y = (cid:82)(cid:82) A q (cid:12)(cid:12) E t y (cid:12)(cid:12) d x d y (cid:82)(cid:82) ∞−∞ (cid:12)(cid:12) E t y (cid:12)(cid:12) d x d y . (88)Here, the enumerator contains an integration over thecross section area A q of the active region formed bythe quantum systems. The intensity distribution in thewaveguide is given by the time-averaged magnitude ofthe z component of the Poynting vector, which is withEq. (84) obtained as I = |(cid:104) S z (cid:105)| = (cid:12)(cid:12) (cid:60) (cid:8) E ωy ( H ωx ) ∗ (cid:9)(cid:12)(cid:12) / (cid:15) cn eff (cid:12)(cid:12) E ωy (cid:12)(cid:12) / . (89)3 | 〈 S z 〉 | x A q A q / G FIG. 9. Illustration of the transverse electric mode profile ofthe waveguiding structure in Fig. 8. The solid curve showsthe mode profile (cid:12)(cid:12) E t y (cid:12)(cid:12) , and the dashed curve represents theequivalent rectangular mode profile. With Eqs. (84) and (88), we then arrive at the usualdefinition
Γ = (cid:82)(cid:82) A q |(cid:104) S z (cid:105)| d x d y (cid:82)(cid:82) ∞−∞ |(cid:104) S z (cid:105)| d x d y . (90)The meaning of Γ is visualized in Fig. 9: We can repre-sent the field confinement factor as Γ = A q /A eff where A eff is the area covered by an equivalent mode which con-serves (cid:82)(cid:82) ∞−∞ |(cid:104) S z (cid:105)| d x d y , but has a rectangular intensitydistribution with |(cid:104) S z (cid:105)| fixed to the value in the quantumnanostructure. Thus, the optical power can with Eq. (89)be written as P = [ I ] q A eff = (cid:12)(cid:12)(cid:12)(cid:2) E ωy (cid:3) q (cid:12)(cid:12)(cid:12) A eff (cid:15) cn eff / , (91)where [ I ] q and (cid:2) E ωy (cid:3) q refer to the values of I and E ωy inthe quantum nanostructure. b. One-Dimensional Maxwell Equations In the fol-lowing, we regard the E ωy and H ωx fields at a frequency ω as spectral components of time dependent fields E y and H x , and transform Eq. (87) into time domain. Forconvenience, we do not consider the ω dependence of E t y and H t x in Eq. (88), but rather evaluate Γ at thecenter frequency ω c of the optical field. Furthermore,to obtain a form compatible with Eq. (82b), we divide n ( ω ) into a constant part, e.g., the value n ( ω c ) at ω = ω c , and a frequency dependent part ∆ (cid:15) eff = n ( ω ) − n ( ω c ) which describes chromatic waveguidedispersion and gives rise to an extra polarization contri-bution. Considering that multiplications with ω in fre-quency domain correspond to operators i ∂ t in time do-main, we obtain from Eqs. (84), (87) and (48) ∂ z E y = µ ∂ t H x , (92a) ∂ z H x = (cid:15) n ( ω c ) ∂ t E y + σ (i ∂ t ) E y (92b) + Γ n (cid:88) i,j d y,ji ∂ t ρ ij + (cid:15) ∂ t [∆ (cid:15) eff (i ∂ t ) E y ] . Here, the generally frequency dependent conductivity σ ( ω ) = (cid:15) cn eff ( ω ) a ( ω ) (93) is often approximated by σ = σ ( ω c ) . Obviously, ∆ (cid:15) eff and σ must be even functions f ( − ω ) = f ( ω ) to pre-serve the real-valued character of Eq. (92). Furthermore,causality requires that the real and imaginary parts ofthe complex permittivity defined in Eq. (78) fulfill theKramers-Kronig relation, which is, strictly speaking,already violated when modeling a medium as a lossless,frequency independent dielectric with (cid:15) r = (cid:15) r (cid:54) = 1 . Notably Eq. (92) does not explicitly depend on x and y anymore. Thus, it is practical to identify H x ( z, t ) and E y ( z, t ) with the field strengths at the transverse posi-tion of the nanostructure, because then E y can directlybe used in Eq. (16) or (17) to evaluate ∂ t ρ ji ( z, t ) . Forcompleteness, we mention that the longitudinal magneticfield component H z can be obtained from ∇ H = 0 , i.e., ∂ z H z = − ∂ x H x .
2. Transverse Magnetic Mode
For TM modes in slab waveguides, Eq. (77) yields with H t x = 0 and ∂ y = 0 (cid:15) r ∂ x (cid:0) (cid:15) − ∂ x H t y (cid:1) = (cid:18) β − ω c (cid:15) r (cid:19) H t y , (94)and the boundary conditions are given by H t y ( x → ±∞ ) → . The mutual dependence ofthe field components is given by Eq. (80), which yieldswith Eq. (76) i ω(cid:15) (cid:15) r E ωx = ∂ z H ωy = i βH ωy , (95) i ω(cid:15) (cid:15) r E ωz = − ∂ x H ωy . (96)From Eq. (44a), we furthermore obtain with Eq. (76) ∂ z E ωx = i ωµ H ωy + ∂ x E ωz . (97)Using Eq. (96) to eliminate E ωz in Eq. (97) yields withEq. (94) ∂ z E ωx = i ω(cid:15) (cid:15) r β H ωy . (98)In Fig. 10, the fundamental TM mode of a QCL waveg-uide structure is shown. The magnetic field distribution H t y is continuous, while E t x exhibits jumps at interfaces oflayers with different (cid:15) r , as can also be seen from Eq. (95).In the following, we treat the background refractiveindex n of the host material for the quantum systemsas a real constant, assuming that the main frequencydependence and gain/loss in the nanostructure is pro-vided by the quantum systems rather than the host ma-terial. Furthermore, although for example in a quan-tum well structure the barrier and well materials willhave different refractive indices as indicated in Figs. 8(b)and 9, the nanostructured region can still be approxi-mately described by a single effective (cid:15) r since the indi-vidual layers are too thin to be resolved by the opti-cal field. For TM modes, the effective permittivity is4 x [ µ m] -2 0 2 4 6 8 10 12 14 F i e l d s t r eng t h [ a r b . u .] R e f r a c t i v e i nde x FIG. 10. Fundamental mode of a QCL waveguide structure at
28 THz . The solid and dashed curves show (cid:12)(cid:12) E t x (cid:12)(cid:12) and (cid:12)(cid:12) H t y (cid:12)(cid:12) ,respectively. Furthermore, the refractive index profile, i.e., (cid:60) (cid:110) (cid:15) / ( x ) (cid:111) , is displayed. then obtained as the harmonic mean of the individualpermittivity values, i.e., we obtain in the host material n − = (cid:0) ∆ (cid:15) − , + ∆ (cid:15) − , (cid:1) / (∆ + ∆ ) where ∆ i and (cid:15) r ,i denote the total thicknesses and permittivities of the re-gions made from material i = 1 , . Generally, the treatment of TM modes is known to bemore complex than for TE modes.
Specifically, in con-trast to the TE case the derivation of one-dimensionalMaxwell-type equations for the transverse field compo-nents is not as straightforward as in Section V A 1. Thiscan for example be seen from Eq. (97), which only as-sumes a Maxwell-type form in analogy to Eq. (92a) if ∂ x E ωz can be neglected. Similarly as in Section V A 1, weassume that the fields are approximately constant overthe transverse cross section of the nanostructured region,and identify in the following E ωx and H ωy with the fieldstrengths at the transverse position of the nanostructure.With (cid:15) r = n , we obtain from Eq. (98) ∂ z E ωx = i ωµ n ( ω c ) n H ωy , (99)where we have assumed that | β (cid:48) | (cid:28) | β | and approximated β ≈ β = n ω /c as in Section V A 1, with the ef-fective waveguide refractive index n eff . Furthermore, wehave neglected the frequency dependence of n eff , evaluat-ing it at the center frequency ω = ω c of the optical field,so as to formally obtain a Maxwell-type equation witha frequency independent effective relative permeability µ eff = n ( ω c ) /n . In order to complete our model,Eq. (99) must be complemented by a second Maxwell-type equation with a form similar to Eq. (85) in frequencydomain, i.e., Eq. (92b) in time domain. Importantly, thisequation has to include the losses and frequency depen-dence omitted in Eq. (99), so that the correct field prop-agation dynamics is obtained. Specifically, from Eq. (76)we obtain the field propagation equations in frequency domain ∂ z E ωx = − β E ωx , ∂ z H ωy = − β H ωy . This requiresthat ∂ z H ωy = i n ωµ n ( ω c ) β E ωx , (100)as can be verified by differentiating Eq. (99) with respectto z and eliminating H ωy with Eq. (100), or alternativelyeliminating E ωx in an analogous way.As in Section V A 1, the polarization due to the quan-tum systems is again perturbatively included in terms ofa change ∆ β to β . To this end, we re-derive Eq. (94)from Maxwell’s equations, Eq. (44), but now keep thepolarization contribution, which yields on the left side ofEq. (94) the perturbation term βωP ωx . Here, P ωx is the x component of P q in frequency domain, while a pos-sible additional z component has been neglected. WithEq. (95) and β ≈ β = n ω /c , the perturbation termcan then be written as βωP ωx ≈ ∆ LH t y with ∆ L ≈ ω c n n P ωx (cid:15) E ωx , (101)and Eq. (94) becomes (cid:16) ˆ L + ∆ L (cid:17) H t y = β H t y with ˆ L = (cid:15) r ∂ x (cid:15) − ∂ x + ω c − (cid:15) r . Similarly as in Section V A 1,we assume that the device operates in a given trans-verse mode with propagation constant β , and use thatfirst order perturbation theory does not affect the cor-responding eigenfunction H t y . Since ˆ L is non-Hermitian,a biorthogonal basis set must be used, and the changeof β is given by ∆ β = (cid:10) ˜ φ (cid:12)(cid:12) ∆ L (cid:12)(cid:12) φ (cid:11) / (cid:10) ˜ φ (cid:12)(cid:12) φ (cid:11) . Here φ = H t y , while ˜ φ denotes the corresponding eigenfunction ofthe adjoint problem ˆ L † ˜ φ ( x ) = (cid:0) β (cid:1) ∗ ˜ φ ( x ) , with ˆ L † = ∂ x (cid:0) (cid:15) − (cid:1) ∗ ∂ x ( (cid:15) ∗ r . . . ) + ω c − (cid:15) ∗ r and ˜ φ ( x → ±∞ ) → . Ascan be seen by inserting Eq. (95) into Eq. (94), ˜ φ sim-ply corresponds to the conjugate complex electric fielddistribution of the mode (cid:0) E t x (cid:1) ∗ , and in analogy to Sec-tion V A 1 a we then obtain the field confinement factor Γ = (cid:82)(cid:82) A q E t x H t y d x d y/ (cid:82)(cid:82) ∞−∞ E t x H t y d x d y . Generally, thisexpression is complex, but for real (cid:15) r it coincides withthe previous result Eq. (90). Thus, the expression forthe confinement factor given in Eq. (90) corresponds tothe perturbative expression for both TE and TM modesin the case of real (cid:15) r , and is commonly also used forcomplex (cid:15) r where it can be seen as a real-valued ap-proximation to the perturbative result. Similarly as inSection V A 1, we insert β ≈ β + i | β | a + ∆ β with β ( ω ) = sgn ( β ) ωn eff ( ω ) /c into Eq. (100), which yields ∂ z H ωy =i ω (cid:15) n n ( ω c ) n ( ω ) E ωx − (cid:15) n n ( ω c ) cn eff ( ω ) aE ωx + i ω Γ P ωx . (102)Here, we have neglected a possible frequency dependenceof n eff and Γ in the last term. Obviously, our two derivedMaxwell-type equations, Eqs. (99) and (102), have thesame form as Eqs. (84) and (87) for TE modes, as can5be seen by substituting H ωx → − H ωy , E ωy → E ωx , (cid:15) → (cid:15) n /n ( ω c ) , µ → µ n ( ω c ) /n . Thus, the Maxwell-type equations in time domain can be obtained in thesame way as Eq. (92), yielding ∂ z E x = − µ n ( ω c ) n ∂ t H y , (103a) ∂ z H y = − (cid:15) n ∂ t E x − σ (i ∂ t ) E x (103b) − Γ n (cid:88) i,j d x,ji ∂ t ρ ij − (cid:15) ∂ t [∆ (cid:15) eff (i ∂ t ) E x ] , where the generally frequency dependent conductivity σ ( ω ) = (cid:15) n cn eff ( ω ) a ( ω ) /n ( ω c ) (104)is often approximated by σ = σ ( ω c ) . Furthermore, in thelast term describing chromatic waveguide dispersion, wenow have ∆ (cid:15) eff = n n ( ω ) /n ( ω c ) − n . As discussedbelow Eq. (93), certain conditions apply to ∆ (cid:15) eff and σ .In particular, they must be even functions f ( − ω ) = f ( ω ) to preserve the real-valued character of Eq. (103). Identi-fying E x ( z, t ) and H y ( z, t ) with the field strengths at thetransverse position of the nanostructure, E x can directlybe used in Eq. (16) or (17) to evaluate ∂ t ρ ij ( z, t ) . Asstated above, Eq. (103) has been constructed to assumethe form of one-dimensional Maxwell equations and toyield the correct propagation behavior for E x and H y .On the other hand, the relation between E x and H y ,given by Eq. (95) for P ωx = 0 , is in Eq. (103) for the gen-eral case of waveguide loss and dispersion only approxi-mately fulfilled. B. Slowly Varying Amplitude Approximation
As in Section V A, we assume a waveguiding struc-ture which is invariant in propagation direction z , andnow employ the slowly varying amplitude approxima-tion. The guided mode solutions at a given frequency ω are characterized by the propagation constant and a z independent transverse field distribution F ( x, y ) . Here F is induced by the refractive index profile ∆ n ( x, y ) ,while the polarization of the quantum systems and othernonlinear effects are assumed to act as perturbationswhich do not significantly affect the transverse fielddistribution. We start from Eq. (60) and introduce the slowly vary-ing field envelopes by inserting Eq. (61), where we how-ever replace k c by the propagation constant of the guidedmode β = β ( ω c ) . In the following, it is advantageous toswitch to the spectral domain, where the slowly varyingenvelopes depend on the frequency variable ∆ ω = ω − ω c ,corresponding to the frequency offset from ω c . Neglect-ing higher order derivatives of t and z in the spirit ofthe SVAA and the paraxial approximation as describedin Section IV B, and considering that time derivatives ofthe envelopes are replaced by multiplications with − i∆ ω in frequency domain, we obtain − β ∂ z E + β E = ∇ E + ω c n E + ω µ P . (105)For the term ω c − n E , we have retained the fullfrequency dependence of the complex refractive in-dex n ( x, y, ω ) to include chromatic dispersion, asdescribed further below. More specifically, n = n (1 − n )+i σ/ ( ω(cid:15) ) contains the refractive index pro-file via ∆ n ( x, y ) and losses via the conductivity σ , where n , ∆ n and σ may be treated as frequency dependent.Assuming a guided mode solution, we can use the sepa-ration ansatz E ( x , t ) = e E ( z, t ) F ( x, y ) , (106)with the polarization direction of the electric field e andmodal distribution F . Inserting Eq. (106) into Eq. (105),multiplying by e and introducing the separation constant β , the right side of the resulting equation becomes (cid:18) β − ω c n (cid:19) F = ∇ F, (107)which does not depend on the propagation coordinate z .As in Section V A 1, P is subsequently included basedon first order perturbation theory. Equation (107),together with the boundary condition that F → for x + y → ∞ , constitutes an eigenvalue equation for F with complex eigenvalues β , featuring multiple eigen-solutions which correspond to the different transversewaveguide modes. In the following, we assume that thedevice operates in a single transverse mode, possibly thefundamental mode. As in Section V A 1, we split thecomplex propagation constant β into a real and an imag-inary part β = β +sgn ( β ) i a/ with power loss coefficient a , and assume that | a | (cid:28) | β | . Including P in first or-der perturbation theory in analogy to Section V A 1 doesnot alter F , but yields a modified propagation constant β + ∆ β with ∆ β = ω c β (cid:82)(cid:82) ∞−∞ ∆ (cid:15) r | F | d x d y (cid:82)(cid:82) ∞−∞ | F | d x d y ≈ ω c β Γ∆ (cid:15) r , (108)where ∆ (cid:15) r = eP / ( (cid:15) E ) and Γ = (cid:82)(cid:82) A q | F | d x d y/ (cid:82)(cid:82) ∞−∞ | F | d x d y in agreement withEq. (88). Here we have evaluated ∆ β at the carrierfrequency ω c , in accordance with the SVAA. For arealistic description of guided mode propagation, thefrequency dependence of β itself should however beretained, giving rise to chromatic dispersion. Thiseffect is commonly described in terms of a Taylor series, β (∆ ω ) = (cid:80) n ( β n /n !) ∆ nω with β n = [d nω β ] ω = ω c . Whilefrequency dependent waveguide loss can be includedin a similar manner by an ω dependent coefficient a ,we ignore this effect since usually the spectral gainor loss profile is dominated by the contribution of the6quantum systems, contained in P . Furthermore assum-ing that | ∆ β | (cid:28) | β | and β ≈ β , we can approximate β − β ≈ β ( β − β ) + i | β | a + 2 β ∆ β . With thisresult and Eq. (108), the separation ansatz yields for theleft-hand side of Eq. (105) in time domain v g ∂ t E + ∂ z E =i (cid:88) n ≥ β n n ! (i ∂ t ) n E − sgn ( β ) a E + i ω (cid:15) β c Γ eP , (109)where v g = β − denotes the group velocity at ω c .The guided field solution is characterized by a linearlypolarized field distribution, with the electric field point-ing in direction e as reflected by the ansatz for the elec-tric field, Eq. (106). The corresponding modes aretransverse electromagnetic, i.e., with transverse, perpen-dicular electric and magnetic fields. Notably, this ap-proach always yields two degenerate modes, orthogonallypolarized in transverse x and y directions. In reality, thisapplies for example to an ideal, cylindrically symmet-ric single-mode fiber, while irregularities such as randomvariations in the core shape already break the degener-acy. Within the assumptions of weak waveguiding, theoptical power is given by the corresponding expressionfor the TE mode, Eq. (91).Above approach is commonly used to model coher-ent propagation effects in optical fibers like self-inducedtransparency, where the dopants, such as erbium ions,take the role of the quantum systems modeled by theBloch equations, and the host material is for exampleglass. Here, in addition to the refractive index pro-file, fiber loss and chromatic dispersion, also other effectsrelated to the host material are commonly considered.This in particular includes optical nonlinearity due to anintensity dependent refractive index of the host material,which induces an intensity dependent phase shift of theoptical field and is thus referred to as self-phase modu-lation. This effect can be included in Eq. (105) by sub-stituting n with (cid:16) n + n | E | (cid:17) ≈ n + 2 n n | EF | . Treating the nonlinear component as a perturbation, wecan again use Eq. (108) with ∆ (cid:15) r = 2 n n | EF | and in-clude this effect in a similar manner as discussed above.With Eqs. (109), (64) and (88), we finally obtain thepropagation equation v g ∂ t E + ∂ z E = i (cid:88) n ≥ β n n ! (i ∂ t ) n E − sgn ( β ) a E + i γ | E | E + i n ω (cid:15) β c Γ (cid:88) ω ij > d ji η ij , (110)with the self-phase modulation coefficient γ = n n ω β c (cid:82)(cid:82) ∞−∞ | F | d x d y (cid:82)(cid:82) ∞−∞ | F | d x d y . (111) We note that in Eq. (111), the nonlinearity is assumedto extend over the whole fiber cross section since boththe core and cladding typically consist of the same hostmaterial. The MB equations are then obtained by cou-pling Eq. (110) to the Bloch equations in RWA, Eq. (43).Here it is practical to normalize F in Eq. (106) so that F = 1 at the transverse position of the dopants actingas quantum systems; then the field in Eq. (43) is directlygiven by E = e E ( z, t ) .Typically, the MB equations are stepped in time toobtain the temporal evolution of the optical field in agiven geometry. For the case of unidirectional prop-agation along a fiber where the input at z = 0 is agiven time-limited optical waveform such as a pulse, itis more practical to propagate the field in z direction. Itis then convenient to introduce the retarded time variable τ = t − z/v g , which is defined with respect to a time framewhich co-propagates with the waveform. Denoting theposition variable in the new coordinate system as ζ = z ,we then obtain the partial derivatives ∂ z = ∂ ζ − (1 /v g ) ∂ τ , ∂ t = ∂ τ . Thus, Eq. (110) becomes ∂ z E =i (cid:88) n ≥ β n n ! (i ∂ τ ) n E − sgn ( β ) a E + i γ | E | E + i n ω (cid:15) β c Γ (cid:88) ω ij > d ji η ij , (112)where we have resubstituted ζ with z . For very shortpulses with durations of only a few optical cycles, ad-ditional corrections might have to be included on theright side of Eq. (112). In particular, this includes theself-steepening term − γω − ∂ τ (cid:16) | E | E (cid:17) which isa higher order term dropped in the SVAA, and theRaman-induced frequency shift term − i γT R E∂ τ (cid:16) | E | (cid:17) with Raman-response time T R . In Eqs. (110)and (112), we have assumed that d ji and E are aligned inthe same direction or, as is more realistic for an opticalfiber, d ji are effective dipole moments which average overthe different orientations of the dopant ions with respectto the field. Equation (112) is solved together with theBloch equations in RWA, Eq. (43), which are expressed inthe retarded time frame simply by substituting t with τ in the density matrix elements and derivative operators.The resulting equation system is sometimes also referredto as Hirota–Maxwell–Bloch system. C. Fabry-Pérot Type Resonator
For lasers, optical feedback has to be provided, whichis in semiconductor lasers typically achieved by using aFabry-Pérot type waveguide resonator. Here, the cleavedend facets provide natural reflection due to the refractiveindex jump between the semiconductor material and air.7
1. Boundary Conditions at the End Facets
In the Fabry-Pérot type resonator, the ansatz for theoptical field, Eqs. (61a), is extended to include a forwardand a backward propagating component, with ampli-tudes E + and E − , respectively. Furthermore assuminga guided mode solution as in Eq. (106), we obtain E ( z, t ) = 12 (cid:2) E + ( z, t ) exp (i β z − i ω c t )+ E − ( z, t ) exp ( − i β z − i ω c t ) + c.c. (cid:3) , (113)where β is the real part of the propagation constantat ω = ω c . With the (generally complex) field reflectioncoefficients r and r of the facets, assumed to be locatedat z = 0 and z = L where L is the resonator length, weobtain E + ( z = 0 , t ) = r E − ( z = 0 , t ) ,E − ( z = L, t ) = r E + ( z = L, t ) , (114)where we have neglected a possible frequency dependenceof r , .For the full Maxwell equations, a decomposition of thefield into a forward and a backward propagating com-ponent is not practical. Here, reflecting boundary con-ditions can in principle be implemented by position de-pendent parameters, e.g., by setting n eff = 1 , ∆ (cid:15) eff = 0 , σ = 0 and d y,ji = 0 in Eq. (92b) for z < and z > L if we assume air outside of the resonator region and ne-glect modal effects. In this context, care has to be takento suppress unwanted spurious reflections at the sim-ulation domain boundaries, which can be achieved byimplementing absorbing boundary conditions. However,this is not quite trivial, and various methods with dif-ferent degrees of complexity have been developed. Asimplified treatment, which works best for highly reflect-ing facets, is to use perfectly reflecting boundary condi-tions by setting the transverse electric field componentat the facet positions to zero. The mirror loss, i.e., thedecay of the optical field in the cavity due to outcouplingthrough the mirrors, can then be considered by a dis-tributed power loss coefficient a m , which is is obtainedfrom | r r | = exp ( − a m L ) as a m = − ln ( | r r | ) /L. (115)Using Eq. (93) or (104), σ in Eq. (92b) or (103b) can thenbe determined from the total power loss coefficient a = a m + a w , where a w = 2sgn (cid:0) (cid:60) (cid:8) β (cid:9)(cid:1) (cid:61) (cid:8) β (cid:9) denotes thewaveguide loss. a. Reflection Coefficient Using special reflectivestructures, such as reflection/antireflection coatings ordistributed Bragg reflectors, r , can be custom-tailored.In the following, we focus on the highly relevant casewhere the bare end facets are used as reflective elements.For sufficiently large transverse waveguide dimensions,Fresnel’s formula for normal incidence can be used to es- timate the field reflection coefficient at the facet as r = n eff − n eff + 1 = β − k β + k , (116)with n eff = β/k and k = ω/c . While Eq. (116) is usu-ally valid for weak waveguiding assumed in the deriva-tion of Eqs. (60) and (110), modal effects can result inincreased reflection at the facets. Various methods areavailable to compute the reflectance R = | r | from thetransverse mode profile. For TE polarization, it is practical to decompose thewaveguide mode, characterized by its complex propaga-tion constant β and magnetic field distribution H t x whichcan be computed from Eq. (81), into plane waves, usingthe Fourier transform Φ x ( k x ) = (cid:90) ∞−∞ H t x ( x ) exp ( − i k x x ) d x. (117)Then, a generalized version of Eq. (116) for tilted inci-dence is applied to each plane wave in order to calculatethe reflection coefficient. The reflectance R , i.e., theratio of the optical power reflected at the facet to theincident power, is obtained by integrating over all com-ponents, yielding R = 12 π (cid:20)(cid:90) ∞−∞ (cid:12)(cid:12) H t x (cid:12)(cid:12) d x (cid:21) − (cid:90) ∞−∞ (cid:12)(cid:12)(cid:12)(cid:12) β − κβ + κ (cid:12)(cid:12)(cid:12)(cid:12) | Φ x | d k x , (118)where κ ( k x ) = (cid:112) k − k x and the square root is cho-sen so that (cid:61) { κ } > . From Eq. (84), we see thatEq. (118) can also be evaluated by replacing H t x with E t y in Eqs. (118) and (117). Making the reasonable as-sumption that (cid:12)(cid:12) (cid:60) (cid:8) β (cid:9)(cid:12)(cid:12) (cid:29) (cid:12)(cid:12) (cid:61) (cid:8) β (cid:9)(cid:12)(cid:12) , we can approximatelytreat the reflection coefficient r as real-valued. Further-more assuming that the share of reflected power goinginto other waveguide modes is negligible, we obtain r = R / .The case of TM polarization is somewhat more com-plex and can be treated based on the boundary valuemethod. Starting from the magnetic field distribution H t y given by Eq. (94), we first evaluate the power trans-mittance T through the facet, T = 2 π (cid:12)(cid:12) β (cid:12)(cid:12) (cid:20)(cid:90) ∞−∞ | (cid:15) r | − (cid:60) (cid:8) β ∗ (cid:15) r (cid:9) (cid:12)(cid:12) H t y (cid:12)(cid:12) d x (cid:21) − × (cid:90) k − k κ | Φ y | (cid:12)(cid:12) Φ (cid:48) y (cid:12)(cid:12) (cid:12)(cid:12) κ Φ y + β Φ (cid:48) y (cid:12)(cid:12) d k x . (119)Here, Φ y ( k x ) denotes the Fourier transform, Eq. (117),of H t y ( x ) , and Φ (cid:48) y ( k x ) is the Fourier transform appliedto the function (cid:15) − ( x ) H t y ( x ) , with the complex relativepermittivity profile of the slab waveguide structure (cid:15) r ( x ) .Again neglecting modal effects and assuming a real r , thereflection coefficient is then given by r = (1 − T ) / .8 (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8) z (cid:9)(cid:4)(cid:7) (cid:10)(cid:6)(cid:7)(cid:8) FIG. 11. Illustration of spatial hole burning in a Fabry-Pérotcavity.
Since the field distributions and β in Eqs. (118) and(119) depend on ω , this also applies to the obtained re-flection coefficients. Usually, this frequency dependenceis neglected in the formulation of the boundary condi-tions, and r is taken at the center frequency ω c of theoptical field.
2. Spatial Hole Burning
In a Fabry-Pérot resonator, the reflection at the endfacets gives rise to counterpropagating waves, which pro-duce a standing wave pattern with a periodicity corre-sponding to the wavelength. At the field node positions,there is no interaction of the optical field with the quan-tum systems. This also implies that the population in-version and resulting optical gain, as provided by thequantum systems in the active region of a semiconduc-tor laser, do not get saturated at those positions. Thus,other modes at slightly different frequencies which havetheir maxima close to these unsaturated regions can alsostart lasing. In Fig. 11, this effect is illustrated, which isreferred to as (longitudinal) spatial hole burning (SHB).The resulting multimode lasing can be desired or unde-sired, depending on the envisaged application. For exam-ple, the broadening of the lasing spectrum is beneficial inapplications such as the generation of frequency combs inQCLs, which are comb-like optical spectra used for pre-cision metrology and sensing. On the other hand, spa-tial hole burning tends to introduce optical instabilitiesin form of irregular variations in the mode amplitudesand phases.
In a similar way as just discussed forthe propagation direction, SHB can also occur along thetransverse directions, and has been shown to affect thespatiotemporal dynamics especially in broad-area semi-conductor lasers
The inversion grating is smoothed out by carrier dif-fusion processes, and SHB can even be neglected in afirst approximation if diffusion is strong enough.
Dif-fusion can be generically described by adding a term ∂ t ρ ii | diff = ∇ D i ∇ ρ ii to the Bloch equations, Eq. (16) or(17). Here D i ( x ) is the diffusion coefficient associ-ated with level i . In the following, we focus on longitudi-nal SHB. Furthermore assuming constant coefficients D i ,the diffusion term added to the Bloch equations, Eq. (16) or (17), becomes [ ∂ t ρ ii ] diff = D i ∂ z ρ ii , (120)with zero-flux boundary conditions ∂ z ρ ii = 0 at the res-onator ends. In a two-level description of bulk semicon-ductors, the levels i correspond to the conduction andvalence bands, and the diffusion process is largely medi-ated by carrier-phonon and carrier-carrier scattering be-tween the k states in the bands. In quantum wells,the levels correspond to the subbands formed by one-dimensional carrier confinement, and the diffusion pro-cess is mediated by scattering between the k states inthe subbands. In multi-quantum-dot structures, theSHB dynamics is often modeled by taking into accountthe carrier diffusion in the wetting layer, as well as carriercapture and escape processes to and from the quantumdots.
Since these processes effectively reduce thediffusion length, SHB can have a strong effect, similarlyas in intersubband devices like QCLs which typically fea-ture a very fast gain recovery dynamics.
On the otherhand, the inversion grating is usually eliminated in in-terband bulk and quantum well lasers due to effectivediffusion. a. Slowly Varying Amplitude Approximation
As forEq. (42), we assume that all transitions between pairs ofstates i and j with non-negligible coupling to the opticalfield are in near-resonance, | ω ij | ≈ ω c . For the corre-sponding off-diagonal density matrix elements, we nowmake the ansatz ρ ij = η + ij ( z, t ) exp [sgn ( ω ij ) i ( β z − ω c t )]+ η − ij ( z, t ) exp [sgn ( ω ij ) i ( − β z − ω c t )] . (121)The periodicity of the inversion grating corresponds tothat of the optical intensity, i.e., half the wavelength.Thus, for the populations we make the ansatz ρ ii = ρ ii ( z, t ) + ρ + ii ( z, t ) exp (2i β z )+ ρ − ii ( z, t ) exp ( − β z ) , (122)where ρ + ii = (cid:0) ρ − ii (cid:1) ∗ correspond to the inversion grating’samplitudes. An analogous ansatz with ρ ij and ρ ± ij is alsochosen for off-diagonal density matrix elements which areassociated with two closely aligned levels and are thusnot treated in RWA, such as resonant tunneling transi-tions in QCLs. Inserting Eqs. (113), (121) and (122)into Eq. (17) with the diffusion term Eq. (120) added toEq. (17b), we obtain in a similar way as described in Sec-9tion III E ∂ t η ± ij = i∆ ij η ± ij + i2 (cid:126) (cid:0) ρ jj − ρ ii (cid:1) d ij (cid:26) E ± E ∗± (cid:27) + i2 (cid:126) d ij (cid:26) (cid:0) ρ ± jj − ρ ± ii (cid:1) E ∓ (cid:0) ρ ∓ jj − ρ ∓ ii (cid:1) E ∗∓ (cid:27) − γ ij η ± ij , ω ij (cid:26) > < (cid:27) ,∂ t ρ ii = 1 (cid:126) (cid:88) jω ij > (cid:61) (cid:8) d ji (cid:0) η + ij E ∗ + + η − ij E ∗− (cid:1)(cid:9) + 1 (cid:126) (cid:88) jω ij < (cid:61) (cid:8) d ji (cid:0) η + ij E + + η − ij E − (cid:1)(cid:9) + (cid:88) j (cid:54) = i r j → i ρ jj − r i ρ ii ,∂ t ρ + ii = i2 (cid:126) (cid:88) jω ij > (cid:0) d ij η − ji E + − d ji η + ij E ∗− (cid:1) + i2 (cid:126) (cid:88) jω ij < (cid:0) d ij η + ji E ∗− − d ji η − ij E + (cid:1) + (cid:88) j (cid:54) = i r j → i ρ + jj − r i ρ + ii − β D i ρ + ii , (123)with ∆ ij = sgn ( ω ij ) ( ω c − | ω ij | ) . Furthermore, we have ρ − ii = (cid:0) ρ + ii (cid:1) ∗ , η ± ji = (cid:0) η ± ij (cid:1) ∗ .Assuming that the coefficients β n in Eq. (110) refer toforward propagation, i.e., β > and β = v − > , thebackward propagating field is described by coefficients − β n . Furthermore deriving the polarization term analo-gously to Eq. (64), we can summarize the equations forthe forward and backward propagating fields as v g ∂ t E ± ± ∂ z E ± =i (cid:88) n ≥ β n n ! (i ∂ t ) n E ± − a E ± + i n ω (cid:15) β c Γ (cid:88) ω ij > d ji η ± ij , (124)where we have neglected the self-phase modulation term.Equation (124) has to be complemented by the boundaryconditions Eq. (114), which together with Eq. (123) con-stitute the MB model in RWA and SVAA for a waveguideresonator.This treatment of SHB can be extended by consider-ing higher spatial frequencies of the inversion grating inEq. (122). Furthermore, it has been suggested to con-sider the formation of a grating, and its relaxation dueto diffusion, also for the off-diagonal density matrix ele-ments in Eq. (121).
VI. ANALYTICAL SOLUTIONSA. Rotating Wave Approximation
1. Monochromatic Excitation
The Bloch equations in RWA, Eq. (43), are in princi-ple analytically solvable for monochromatic excitation,corresponding to a time-constant field envelope.
This is usually achieved by using the Laplace transform,which takes a time dependent function f ( t ) to a func-tion L { f } ( s ) of a complex frequency variable s . Themain advantage is that differentiation becomes a mul-tiplication with s , i.e., L { ∂ t f } = s L { f } − f ( t = 0+) .Restricting ourselves to a two-level system with initialconditions w = w ( t = 0) , η = η ( t = 0) , and con-sidering that the transform is linear and L { } = s − ,Eq. (73) becomes in Laplace domain L { η } = s − − (cid:0) − i L { w } Ω / η (cid:1) , (125a) L { η } = s − (cid:104) i L { w } Ω ∗ / (cid:0) η (cid:1) ∗ (cid:105) , (125b) L { w } = 2 ss + s + s − ( γ w eq + sw ) s (cid:104) ss + s − + ( s + γ ) | Ω | + γ s + s − (cid:105) , (125c)with s ± = ( s + γ ± i∆) and s = (cid:61) (cid:8) η Ω ∗ (cid:9) ( s + γ ) +∆ (cid:60) (cid:8) η Ω ∗ (cid:9) . Here, the off-diagonal matrix elementshave already been eliminated in Eq. (125c) by insert-ing Eqs. (125a) and (125b). Back-transformation ofEq. (125c) is achieved by performing a partial fraction de-composition, which results in a sum of simpler fractionswith known inverse Laplace transforms. This requiresfinding the poles of the rational function in Eq. (125c),given by s = 0 and the roots of the cubic function of s inthe square brackets of the denominator, for which closedanalytical expressions are readily available. The solutionfor f ( t ) = w ( t ) and f ( t ) = η ( t ) is of the form f ( t ) = A + B exp ( − at )+[ C cos ( ωt ) + D sin ( ωt )] exp ( − bt ) , where the decay constants a , b and the oscillation fre-quency ω are the same for w and η , but the coefficients A , B , C and D are different and also depend on the initialconditions. a. Rabi Oscillations In the following, we considerthe case of dissipationless light-matter interaction, i.e., γ = γ = 0 . Then, Eq. (125c) becomes L { w } = 2 s + (cid:0) s + ∆ (cid:1) w s ( s + iΩ g ) ( s − iΩ g )= As + ( w − A ) s + 2 (cid:61) (cid:8) η Ω ∗ (cid:9) s + Ω , (126)where Ω g = (cid:16) ∆ + | Ω | (cid:17) / denotes the generalizedRabi frequency for detuned excitation, and A = Time | Ω | t I n v e r s i on w -1-0.500.51 FIG. 12. Time dependent inversion for monochromatic exci-tation with a field of amplitude Ω and a detuning of ∆ = 0 (solid line), ∆ = | Ω | / (dashed line), ∆ = | Ω | (dash-dottedline), and ∆ = 2 | Ω | (dotted line). Ω − (cid:0) (cid:60) (cid:8) η Ω ∗ (cid:9) + ∆ w (cid:1) . Inverse Laplace transfor-mation of Eq. (126) yields w ( t ) = A +( w − A ) cos (Ω g t )+2Ω − (cid:61) (cid:8) η Ω ∗ (cid:9) sin (Ω g t ) . (127)As can be seen from Eq. (127), a monochromatic, near-resonant light field interacting with an ideal, dissipation-less two-level system causes Rabi flopping, i.e., an os-cillation of the population between states and withfrequency Ω g , as predicted by I. I. Rabi for the analo-gous case of a two-level system in a rotating magneticfield. In Fig. 12, w ( t ) is shown for the initial condi-tions w = − , η = 0 and different detunings. As canbe seen, complete population inversion with w = 1 isachieved only for resonant excitation.Due to the presence of dissipation, above presentedanalytical treatment of Rabi flopping can rarely be di-rectly used for the description of optoelectronic deviceoperation. However, under favorable conditions, signa-tures of Rabi oscillations have been observed in nanos-tructured optoelectronic systems and devices. Theseincludes quantum well structures, QCLs, singlequantum dots, and nanowire lasers at cryogenictemperatures, as well as quantum dot and quan-tum dash amplifiers at room temperature. Besidesthe usually strong influence of dissipation, effects be-yond the two-level dynamics and inherent restrictionsof models based on macroscopic Maxwell-Bloch equa-tions, the applicability of Eq. (126) is also limited bythe validity range of the RWA. In particular, for verystrong optical excitation where the Rabi frequency ap-proaches the optical resonance frequency, effects beyondthe RWA have been observed in bulk and nanostructuredsemiconductors. b. Steady State Solution In the following, we con-sider the steady state behavior for a dissipative two-levelsystem under monochromatic optical excitation, i.e., fora field at frequency ω c with constant amplitude. In the presence of dissipation, the coherent transients associatedwith above discussed Rabi oscillation, Eq. (127), decay,and the system approaches the steady state for t → ∞ .The steady state solution can be obtained by setting ∂ t = 0 in Eq. (73). Alternatively, we can apply the finalvalue theorem to Eq. (125), stating that if lim t →∞ f ( t ) exists, it is identical to lim s → s L { f } . Introducing therelaxation times T , = γ − , , we then obtain η = 12 Ω T (∆ T − i) w eq T + T T | Ω | , (128a) w = (cid:0) T (cid:1) w eq T + T T | Ω | . (128b)From Eq. (128), an expression for the relative per-mittivity (cid:15) r and susceptibility χ = (cid:15) r − due to thequantum systems can be derived. Setting the classi-cal expression for the complex polarization amplitude P = (cid:15) χ E equal to Eq. (75), we obtain with Eq. (128a), Ω = (cid:126) − d E and Eq. (89) the frequency and intensitydependent susceptibility χ = n | d | T (cid:15) (cid:126) (∆ T − i) w eq T + I/I s , (129)with the saturation intensity at zero detuning I s = (cid:126) (cid:15) n eff c T T | d | . (130)For arbitrary detuning, the saturation intensity is thengiven by I s (cid:0) T (cid:1) , i.e., non-resonant fields inter-act less strongly with the two-level system, and thus sat-uration occurs for higher intensities. As in Section V B,it is here assumed that d and E are aligned in the samedirection, or that d is an effective value averaged overthe different orientations of, e.g., dopant ions in an opti-cal fiber.In Fig. 13, the real and imaginary parts of χ are shownas a function of detuning from the optical resonance fre-quency for various optical intensities. (cid:60) (cid:8) χ (cid:9) , which con-tains chromatic dispersion, changes sign at the resonancefrequency. (cid:61) (cid:8) χ (cid:9) describes gain for w eq > , i.e., positivepopulation inversion, and loss for w eq < , where the fre-quency dependence is given by a Lorentzian profile withthe full width at half-maximum (FWHM) bandwidth γ (1 + I/I s ) / . For increased intensities, the profilethus gets broadened which is known as power broaden-ing, and also the peak value at resonance frequency isreduced by a factor of I/I s , which corresponds togain saturation for w eq > and saturable absorption for w eq < .In the following, we investigate optical power amplifi-cation or absorption by two-level systems. For TE modesor generally in the limit of weak waveguiding, the poweris given by Eq. (91). Multiplying Eq. (110) from left with E ∗ and adding the complex conjugate, we obtain withEqs. (128a) and (91) and β = ω c n eff /c∂ z P = − aP + gP, (131)1 ℜ { χ } / χ -0.500.5 T ∆ -6 -3 0 3 6 ℑ { χ } / χ s =0I/I s =1I/I s =3(a)(b) FIG. 13. Complex susceptibility of the two-level system, nor-malized to χ = − n | d | T (cid:15) − (cid:126) − w eq , as a function ofthe detuning frequency ∆ . Shown is the (a) real and (b)imaginary part for various optical intensities. where a is the waveguide loss coefficient, and the two-level power gain coefficient is given by g = α w eq T + P/P s , (132)with α = Γ ω c n | d | T (cid:15) c (cid:126) n eff (133)and the saturation power at zero detuning P s = A eff I s .With the help of the Lambert W function, defined by x =W ( x ) exp [W ( x )] , we can write the solution of Eqs. (131)and (132) for zero waveguide loss, a = 0 , as P ( z ) = P G ( z ) with the power gain factor G ( z ) = P (cid:48) s P W (cid:20) P P (cid:48) s exp (cid:18) P P (cid:48) s (cid:19) exp (cid:18) αw eq T z (cid:19)(cid:21) , (134)where P = P ( z = 0) and P (cid:48) s = (cid:0) T (cid:1) P s . WithEq. (134), the steady state field solution of Eq. (110) canthen for a = γ = 0 be written as E ( z ) = E ( z = 0) [ G ( z )] (1+i∆ T ) / , (135)and the density matrix elements are with Ω = (cid:126) − d E given by Eq. (128). In the exponent of Eq. (135), α H = − ∆ T corresponds to the Henry or linewidth enhance-ment factor, which relates phase changes to changesin the optical gain.In Fig. 14, the optical power gain and phase shift ofthe electric field are shown as a function of propaga-tion distance for different values of initial power P for w eq > , i.e., amplification. In the small signal limit, G = P / P z/ [(1 + ∆ T ) α − w − eq ] [ φ − φ ( ) ] / ( ∆ T ) P /P ′ s = 0 . P /P ′ s = 0 . P /P ′ s = 1 (a)(b) FIG. 14. (a) Optical power gain and (b) phase shift of theelectric field as a function of propagation distance for differentvalues of initial power. P (cid:28) P (cid:48) s , the typical exponential increase in power isobserved. This can also be seen from Eq. (134), whichyields with W ( x ) ≈ x for x (cid:28) the usual exponentialamplification (for G > ) or loss (for G < ) characteris-tics. In the saturation regime, the power increases onlylinearly. Physically, this is a consequence of the fact thatthe growth in optical power is ultimately limited by thesupplied pump power.The expressions Eqs. (129) and (132) for the sus-ceptibility χ and optical gain g are widely used tomodel the optical properties of homogeneously broadenedatomic and nanostructured optical media. Inparticular for interband transitions in bulk semiconduc-tor and quantum well media, the electron wavevectormust be explicitly considered, along with additional cor-rections due to Coulomb interactions. Above deriva-tion of an expression for the susceptibility from the Blochequations can be extended to more than two levels, whichis for example relevant for the investigation of slow lightpropagation. This is usually achieved based on electro-magnetically induced transparency (EIT), where a con-trol laser beam induces a narrow transparency windowwith an extremely low group velocity in the absorptionspectrum of a suitable medium.
EIT requires athree-level configuration, and expressions for the suscep-tibility have been derived in a similar way as above.
A reduction of the group velocity to subsonic speeds,as well as complete halting of light, has been demon-strated in an ultra-cold atomic vapor.
Possible appli-cations include optical buffers, imaging, and quan-tum memory. In view of a future commercializationof these technologies, a compact solid-state based im-plementation is desirable, and slow light propagationas well as light trapping has meanwhile been demon-strated in doped crystals.
The realization of slow2light in suitably engineered semiconductor structures isespecially attractive. Here, the exploitation of tunnelinginduced transparency is highly promising, which differsfrom EIT in that it does not require an optical controlfield, but utilizes strong tunneling coupling between apair of states. For this case, the susceptibility has beenanalytically derived for quantum dot and intersubbandquantum well systems.
Furthermore, nonlinear opti-cal mixing effects which involve optical field contributionsat two or more frequencies, and often rely on more thantwo energy levels, are exploited in many semiconductor-based applications, requiring a description by higher or-der susceptibilities.
The corresponding expressionscan for example be obtained from the Bloch equations byemploying time dependent perturbation theory.
2. Self-Induced Transparency
In addition to above presented steady state solutionto the Maxwell-Bloch equations, also dynamic solutionsare available for some special cases. An important ex-ample is self-induced transparency (SIT), where a spe-cial optical pulse solution exists which can propagatethrough the two-level medium without being attenuatedor disturbed. This effect was theoretically predicted,and first experimentally demonstrated in ruby, by Mc-Call and Hahn.
SIT is based on coherent interac-tion with the medium, which requires that the pulse du-ration must be much shorter than the relaxation pro-cesses described by γ and γ , and thus we can set γ = γ = 0 . Furthermore we assume that the fieldenvelope Ω = (cid:126) − d E is real-valued, and initially re-strict ourselves to resonant excitation, i.e., ∆ = 0 .Then, the solution of Eq. (74) for the initial condition s ( z, −∞ ) = s ( z, −∞ ) = 0 , w ( z, −∞ ) = − canbe written as s = 0 , s = − sin ϑ and w = − cos ϑ with ϑ ( z, t ) = (cid:82) t −∞ Ω ( z, t (cid:48) ) d t (cid:48) , as can easily be veri-fied by re-insertion of the solution into Eq. (74). Thisanalysis can be extended to incorporate inhomogeneousbroadening in media consisting of quantum systems withslightly different resonance frequencies. Assumingthat non-resonant systems with ∆ (cid:54) = 0 essentially re-spond in the same way to Ω as the resonant ones, apartfrom a change in amplitude, we can make the factoriza-tion ansatz s , ∆ ( z, t ) = F (∆) s ( z, t ) , which again yieldsclosed analytical solutions to Eq. (74), s , ∆ ( z, t ) = − F (∆) sin ϑ ( z, t ) , (136) w ∆ ( z, t ) = − F (∆) [cos ϑ ( z, t ) − − . (137)Taking the second derivative of Eq. (74b) for ∆ (cid:54) = 0 andinserting Eqs. (74a), (136) and (137) yields ∂ t ϑ = ∆ F (∆)1 − F (∆) sin ϑ := T − sin ϑ. (138)Since the electric field envelope, and hence also ϑ , is ∆ in-dependent, this must also apply to ∆ F (∆) / [1 − F (∆)] Ω / Ω Normalized retarded time τ /T -4 -2 0 2 4 I n v e r s i on w -101 ℑ { η } -0.500.5 (b)(a) FIG. 15. Self-induced transparency soliton. Shown are (a)the normalized electromagnetic field and (b) inversion w andimaginary part of the off-diagonal matrix element η , where (cid:60) { η } = 0 . which we have thus set equal to a constant T − inEq. (138). This yields a Lorentzian dependence F =1 / (cid:0) T ∆ (cid:1) . Equation (138) corresponds to the pen-dulum problem, where the solutions are given by ellip-tic functions. Here we require Ω = ∂ t Ω = 0 and thus ∂ t ϑ = ∂ t ϑ = 0 at t = ±∞ , yielding the unique solution ϑ = 4 arctan { exp [( t − t ) /T ] } . Introducing t = z/v with the pulse propagation velocity v , we thus obtain Ω ( z, t ) = 2 T − sech ( τ /T ) (139)with the retarded time variable τ = t − z/v .In the optical propagation equation, inhomogeneousbroadening can approximately be included by substitut-ing s − i s with (cid:82) ∞−∞ g (∆) ( s , ∆ − i s , ∆ ) d∆ in the po-larization, Eq. (75), where g (∆) with (cid:82) ∞−∞ g (∆) d∆ = 1 gives the distribution of quantum systems as a functionof the detuning ∆ from ω c . Here, a possible dependenceof the dipole matrix element on ∆ has been neglected.Using s , ∆ = F (∆) s and above result for F (∆) , weobtain s , − ∆ = s , ∆ , and with Eq. (74a) we see that then s , − ∆ = − s , ∆ for Ω i = γ = 0 . Often g (∆) is an evenfunction as further discussed in Section VIII B, and thenthe contribution of s , ∆ cancels out. Under this assump-tion, Eq. (110) becomes without dispersion ( β n = 0 for n ≥ ), loss ( a = 0 ) and self-phase modulation ( γ = 0 ) ∂ t Ω + v g ∂ z Ω = 2 − αγ v g s (cid:90) ∞−∞ g (∆) F (∆) d∆ , (140)where α is given by Eq. (133). Inserting Eqs. (139) and(136) into Eq. (140) yields the pulse propagation velocity v = v g
22 + αγ v g T (cid:82) ∞−∞ g (∆) (1 + T ∆ ) − d∆ . (141)In Fig. 15, the pulse shape, inversion w and the imag-inary part of the off-diagonal matrix element (cid:61) { η } = − s / is shown. Notably, for w eq < when the two-levelmedium normally absorbs light [see Eq. (132)], the opti-cal energy absorbed during the first half of the SIT pulse3and stored in the inversion, is re-emitted during the sec-ond half, which delays the pulse so that v is smaller thanthe group velocity v g without the coherent interaction ascan be seen from Eq. (141), but does not change its shapeor amplitude. This is accompanied by a Rabi flop of thepopulation inversion from w = − to w = 1 and backagain.Generally, based on the area theorem it was foundthat for coherent propagation, the pulse area Θ = (cid:82) Ωd t evolves towards the closest even multiple of π ( Θ =0 , π, π, . . . ) for absorbing media ( w eq < ), and to theclosest odd multiple of π for gain media ( w eq > ). Importantly, the area theorem only makes a statementabout Θ , but does not indicate if the pulse envelopechanges. As discussed above, the SIT pulse Eq. (139),which has a pulse area Θ = 2 π , is the only finite en-ergy solution of Eq. (138) where the pulse envelope ispreserved. However, analytical solutions of the MB equa-tions with changing pulse shapes can also be obtained forother cases of coherent propagation. With regards to novel practical applications, SIT is forexample a highly interesting candidate for the generationof ultrashort optical pulses in various types of lasers withsufficiently long coherence times, such as quantum dotand quantum cascade lasers.
This SIT (orcoherent) mode-locking approach requires a laser designwith one or multiple gain and absorber regions, wherean SIT soliton with a pulse area of d (cid:82) E d t/ (cid:126) = 2 π isapproximately realized in the absorber sections. In or-der to obtain a stable pulse area of π in the gain regions,they are engineered to have half the dipole moment d ofthe absorber sections. Instead of sequential gain and ab-sorber regions, another option is to stack the gain and lossregions in transverse direction, i.e., perpendicular to thepropagation axis. This approach is for example compati-ble with the manufacturing process of QCLs, and an ana-lytical solution has been derived for the steady state pulsesolution. Despite its great promise, SIT mode-lockinghas not been experimentally demonstrated to date.
B. Full-Wave Bloch Equations
Without employing the RWA, the Bloch equations(17) are solvable only for very special conditions. Inparticular, for | ∆ M | = 1 transitions in hydrogen-likeatoms where the dipole matrix element is given by d =2 − / | d | ( e x − i e y ) , excitation with circularly polarizedlight where E x = E ( t ) cos ( ωt ) and E y = E ( t ) sin ( ωt ) leads to dE ( t ) = 2 − / | d | E ( t ) exp ( − i ωt ) . Further-more using Eq. (42) to substitute the off-diagonal den-sity matrix elements in Eq. (17), the resulting equationformally corresponds to the RWA Bloch equation, withthe analytical solutions discussed in Section VI A. Closedanalytical solutions are not available for the relativelysimple, but very important case of monochromatic ex-citation with a linearly polarized wave. Some approxi-mate corrections have been derived, such as the Bloch- Siegert shift describing the change in the system’s res-onance frequency for strong driving, and the Mollowtriplet which refers to the emergence of satellite peaks inthe spectrum of resonantly excited systems.
Interest-ingly, the full Bloch equations can be solved analyticallyif the linearly polarized electric field has the form of an N -soliton. This is also true for the so-called reducedMB equations, which combine the full-wave Bloch equa-tions with a first-order unidirectional optical propagationequation. VII. NUMERICAL SCHEMES
As discussed in Section VI, the full-wave MB equa-tions have known analytical solutions only for very spe-cial cases, and also in the RWA/SVAA approximation, nogeneral analytical solution exists. Therefore computersimulations are in general necessary. From a practicalpoint of view, the numerical scheme should be stable,accurate, and efficient, and a naive discretization will of-ten fail. The goal of this section is to introduce well-established approaches which are straightforward to im-plement, and give a critical discussion of their properties.Furthermore, an overview of recent developments in thefield will be given. Since the RWA/SVAA problem andthe full MB equations are not of the same mathematicalform, their numerical implementation has to be treatedseparately.Several software projects have been published that areable to solve the Maxwell-Bloch equations. For example,the Freetwm tool is an open-source MATLAB codethat simulates the dynamics of semiconductor lasers us-ing the 1D MB equations in rotating wave approximation.The Electromagnetic Template Library (EMTL) is a freeC++ library with Message Passing Interface (MPI) sup-port, which has for example been used to model quan-tum emitters with the full-wave MB equations in two di-mensions.
Another solver library for the full-wave MBequations is the open-source MEEP project, using asimilar representation of the Bloch equations as given inEq. (72). The mbsolve project solves the full-wave MBequations using different parallel acceleration techniquesand features an open-source codebase. Finally, a com-mercial MB solver has been announced.
A. Rotating Wave/Slowly Varying Amplitude Approximation
1. Finite Difference Discretization of the One-DimensionalPropagation Equation
In the following, the numerical solution of the one-dimensional optical propagation equation in the SVAA isdiscussed. Neglecting chromatic dispersion, i.e., setting β n = 0 for n ≥ , we write Eq. (124) in the form ∂ t E ± = ∓ v g ∂ z E ± + f ± ( z, t ) − (cid:96)E ± . (142)4An obvious choice is to use a finite difference discretiza-tion approach where a full spatiotemporal discretizationof E ± onto an equidistant grid with z m = m ∆ z , t n = n ∆ t is imposed. In the following, E ± m,n and f ± m,n denote thenumerical solution of E ± and f ± on the grid. The start-ing point is a Taylor series expansion of E ± ( z m , t n +1 ) around the point ( z m , t n ) , yielding up to second order E ± m,n +1 = E ± m,n + ∆ t ∂ t E ± m,n + ∆ t ∂ t E ± m,n . (143)Then ∂ t E ± m,n and ∂ t E ± m,n are replaced by space deriva-tives: Differentiating Eq. (142) with respect to z , multi-plying the result by ∓ v g and adding it to the time deriva-tive of Eq. (142) yields (cid:0) ∂ t − v ∂ z (cid:1) E ± = (cid:0) ∂ t ∓ v g ∂ z (cid:1) ( f ± − (cid:96)E ± ) . Plugging the result into Eq. (143) and furthermoreusing Eq. (142), we obtain E ± m,n +1 = E ± m,n + ∆ t (cid:110) ∓ v g (cid:2) ∂ z E ± (cid:3) m,n + f ± m,n − (cid:96)E ± m,n (cid:111) + ∆ t (cid:110) v (cid:2) ∂ z E ± (cid:3) m,n + (cid:2) ∂ t f ± (cid:3) m,n ∓ v g (cid:2) ∂ z f ± (cid:3) m,n (cid:111) + (cid:96) ∆ t (cid:110) ± v g (cid:2) ∂ z E ± (cid:3) m,n − f ± m,n + (cid:96)E ± m,n (cid:111) . (144)For finite difference discretization, there are differentpossibilities such as the well known and widely used nd order Lax-Wendroff method, or the Risken-Nummedalfinite differences (RNFD) scheme which was specifi-cally developed in the context of MB simulations. In both cases, (cid:2) ∂ z E ± (cid:3) m,n is approximated by the stan-dard finite difference approximation (cid:0) E ± m +1 ,n − E ± m,n + E ± m − ,n (cid:1) / ∆ z . Here we will treat in detail the RNFDscheme, since it has some advantageous properties as dis-cussed further below. The main difference as comparedto the Lax-Wendroff method is that rather than employ-ing centered differences, depending on the propagationdirection backward/forward finite differences are used,with [ ∂ z E ± ] m,n ≈ ± (cid:0) E ± m,n − E ± m ∓ ,n (cid:1) / ∆ z , [ ∂ z f ± ] m,n ≈± (cid:0) f ± m,n − f ± m ∓ ,n (cid:1) / ∆ z . Furthermore, a time step of ∆ t = ∆ z /v g is imposed. From Eq. (144), we then ob-tain the RNFD scheme E ± m,n +1 = E ± m ∓ ,n + ∆ t (cid:18) f ± m,n + 12 f ± m ∓ ,n − (cid:96)E ± m ∓ ,n (cid:19) + ∆ t (cid:110)(cid:2) ∂ t f ± (cid:3) m,n − (cid:96)f ± m,n + (cid:96) E ± m,n (cid:111) . (145)The term [ ∂ t f ± ] m,n is not substituted with a correspond-ing finite difference approximation since ∂ t η ij can directlybe obtained from the Bloch equations, Eq. (43). In aFabry-Pérot type resonator, Eq. (145) is complementedby the boundary conditions Eq. (114). a. Numerical Properties of the RNFD Scheme Fora numerical scheme to be useful, an important require-ment is that round-off and truncation errors do not get amplified during the computation, since this will even-tually lead to numerical instability. The stability of fi-nite difference discretization schemes can be investigatedbased on a von Neumann stability analysis.
It turnsout that the RNFD scheme is stable for (cid:96) ≥ , which isalso true for the Lax-Wendroff method for a sufficientlysmall Courant number v g ∆ t / ∆ z . On the other hand, forpositive linear gain, i.e., (cid:96) < , we obtain unconditionallyunstable behavior for both schemes. Furthermore, likethe Lax-Wendroff method, the RNFD scheme is secondorder accurate in space and time. This guarantees thatthe numerical scheme converges to the original partialdifferential equation as the grid spacing approaches zero,with a convergence order of two. However, this does notyet guarantee that the numerical solution for finite gridspacing has a physically meaningful behavior, e.g., satis-fies certain physical conservation laws. Thus, additionalconditions might be desirable for a finite difference dis-cretization of Eq. (142), which has the form of an inhomo-geneous scalar convection equation and thus allows us todraw from related work.
Specifically, it has been estab-lished that second and higher order linear finite differenceschemes tend to introduce artificial numerical dispersion,yielding phase errors and numerical oscillations near ex-trema or discontinuities of the solution.
The numeri-cal solution is less prone to phase errors for monotonicitypreserving schemes, which guarantee that for every non-decreasing (non-increasing) initial condition E ± m,n =0 , thenumerical solution at all later instants n > is also non-decreasing (non-increasing). A sufficient condition forthe RNFD scheme to be monotonicity preserving for thehomogeneous propagation equation, i.e., Eq. (142) with avanishing source term f ± ( z, t ) ≡ , can be easily derived:Formulating Eq. (145) for E ± m +1 ,n +1 and subtracting theresulting expression from Eq. (145), we arrive at E ± m +1 ,n +1 − E ± m,n +1 = (1 − (cid:96) ∆ t ) (cid:0) E ± m +1 ∓ ,n − E ± m ∓ ,n (cid:1) + ∆ t (cid:96) (cid:0) E ± m +1 ,n − E ± m,n (cid:1) , which yields ∆ t ≤ /(cid:96) as sufficient condition for mono-tonicity preservation in the stability regime (cid:96) ≥ . This isa unique feature for a second order finite difference prop-agation scheme which is directly related to the choice oftime step ∆ t = ∆ z /v g . Also, this constitutes an im-portant advantage of the RNFD scheme over the Lax-Wendroff method, which does not have this property innumerically stable regions, as can be shown in a similarway as above or directly from Godunov’s order barriertheorem. In Fig. 16, the Lax-Wendroff and the RNFD schemeare compared for lossless propagation of an initially rect-angular pulse without interaction with a quantum sys-tem. For the Lax-Wendroff scheme, spurious oscillatoryfeatures arise in the vicinity of the field discontinuities,which are absent in the RNFD scheme due to its mono-tonicity preserving nature.5
Time [arb. u.]2.8 3 3.2 3.4 3.6 3.8 F i e l d en v e l ope [ a r b . u .] Lax-WendroffRNFD
FIG. 16. Comparison between Lax-Wendroff and RNFDschemes for the forward propagation of a rectangular pulsewith f + = (cid:96) = 0 .
2. Density Matrix Equations
The numerical scheme for the optical propagationequation has to be coupled to a time-propagation schemefor the Bloch equations, Eq. (43). These constitute an or-dinary differential equation (ODE) system describing thetemporal evolution of the density matrix, which has tobe solved for each spatial grid point. In principle, moststandard methods should do the job although they willdiffer in numerical stability, accuracy and efficiency, andwell-established schemes such as Runge-Kutta andAdams-Bashforth have successfully been used. Mostresearch on the suitability of different numerical schemesin literature has focused on the full-wave Bloch equa-tions, as detailed in Section VII B. Since they are iden-tical in structure to the RWA Bloch equations [com-pare for example Eqs. (67a) and (68) to Eq. (73)], theobtained insights should in principle also be valid forthe RWA Bloch equations. The Runge-Kutta methodis further described in Section VII B 4 in the context offull-wave Bloch equations. Here we exemplarily discussthe explicit Adams-Bashforth scheme as an especiallystraightforward to implement and numerically highly effi-cient method. The RNFD scheme for the MB equations,Eq. (145), is strongly coupled, i.e., requires an evaluationof the density matrix and the electric field at the sametime value. The k − step Adams-Bashforth method forthe solution of an ODE system is given by ˆ ρ n +1 = ˆ ρ n + ∆ t k − (cid:88) m =0 c m F n − m (cid:0) ˆ ρ n − m (cid:1) . (146)Here, n corresponds to the time t n = n ∆ t , ∆ t is the timestep size, and F (ˆ ρ ) = L (ˆ ρ ) + D (ˆ ρ ) represents the righthand side of the Lindblad equation (3), and specificallyin our case of the RWA Bloch equations, Eq. (43). Fur-thermore, the c m are suitably chosen coefficients so that maximal accuracy is reached in the approximation.A k − step Adams-Bashforth method has a global numer-ical error on the order of O (∆ kt ) . In this context, itmust be considered that the overall numerical accuracycannot be arbitrarily improved by choosing a high valueof k , since it is also limited by the numerical discretiza-tion of the optical propagation equation, e.g., based onthe RNFD method. As discussed in Section IV C, theBloch equations are initialized by the starting values ofthe density matrix elements at a given time, while theAdams-Bashforth method would require k initial valuesas can be seen from Eq. (146). This problem can for ex-ample be solved by doing the first k − time steps witha different numerical scheme such as the Runge-Kuttamethod, or by initializing the simulations with two-stepAdams-Bashforth on a finer grid. In simulations of laseroperation which are typically started from noise, theexact choice of initial conditions is not critical and thusthe initialization steps required by Adams-Bashforth donot pose a problem. The main advantage is the reducednumerical load as compared to the Runge-Kutta method(see Section VII B 4), which however requires initializa-tion only at a single time point.
3. Generalizations and Alternative Methods
In Section VII A 1, one-dimensional propagation hasbeen assumed, neglecting the transverse coordinates inthe SVAA propagation equation, Eq. (62). In reality,the field dependence, and thus also the temporal evo-lution of the quantum systems, is varying along the x and y coordinates, which must be explicitly consideredfor an inclusion of diffraction and other effects. As long as no transverse boundary conditions or mate-rial dependencies have to be considered, i.e., ∆ n and σ in Eq. (62) are constant, the most straightforward ap-proach is to Fourier-transform Eq. (62) with respect to x and y before the time propagation step is carried out. The resulting equation then depends on z , t and the spa-tial Fourier frequencies k x and k y , converting the deriva-tive operator ∇ into a multiplication with − (cid:0) k x + k y (cid:1) .Thus a one-dimensional propagation method can be used,such as the one discussed in Section VII A 1. Since thisprocedure requires a Fourier transform before and an in-verse transform after each propagation step, the numer-ically efficient fast Fourier transform method is usuallyemployed.As discussed in the context of Eq. (112), for the model-ing of unidirectional fiber or beam propagation often theinitial field at z = 0 is given, and the solution at a certaindistance z = L is required. Then it is more practical topropagate the field in z direction rather than in time, andto introduce the retarded time variable τ = t − z/v g whichsimplifies the propagation operator (cid:0) v − ∂ t E + ∂ z E (cid:1) to ∂ z E . In the absence of other time derivatives, e.g., due tochromatic dispersion, this effectively reduces the propa-gation equation to an ODE. The solution is then marched6in z direction in dependence of τ (and k x,y if applica-ble), and the density matrix is updated after every prop-agation step. The propagation along z can be per-formed with a conventional ODE scheme where for ex-ample the Adams-Moulton method (with the trapezoidalrule as a widely used special case) or Adams-Bashforthmethod, Eq. (146), have been employed, in both casescombined with fourth-order Runge-Kutta for the Blochequations. In the more general case where timederivatives have to be considered in Eq. (112), for ex-ample to incorporate chromatic dispersion, these can behandled in Fourier domain, similarly as for the x and y derivatives discussed in the previous paragraph. One op-tion is to process all terms in Fourier domain, whichhowever complicates the treatment of expressions whichare nonlinear in the field, such as the self-phase modula-tion term in Eq. (112). Another strategy might be to cou-ple the Bloch equations to the split-step Fourier method,which treats only the terms containing time derivativesin Fourier domain, and the others in time domain. B. Full-Wave Simulation
While the RWA significantly reduces the computa-tional workload, care must be taken in cases where its ba-sic assumptions are not fulfilled. For example, the RWAassumes that the electric field intensity is small and thefield spectrum narrow. However, in a scenario where ul-trashort pulse generation is simulated (e.g., mode-lockingoperation in quantum cascade lasers), the electric fieldfeatures high peak intensity and a broad spectrum. Insuch cases, the full electromagnetic wave might have tobe considered in the simulation, and a suitable numeri-cal scheme has to be used. In the following, we describethe methods for the Maxwell and full-wave Bloch equa-tions, Eqs. (44) and (17), which are most widely usedin related literature, and address the coupling betweenthe updates of the electric field and the density matrix.Finally, we assess the advantages and drawbacks of thedifferent methods.
1. Numerical Schemes for Maxwell’s Equations
Out of the many numerical methods that solveMaxwell’s Equations, mainly two – namely the finite-difference time-domain (FDTD) and the pseudo-spectraltime-domain (PSTD) method – are used in the contextof Maxwell-Bloch equations.The FDTD method is one of the standard methods forMaxwell’s equations, and is widely used in combinationwith the optical Bloch equations. Here, the derivatives with respect to time and spaceare approximated using central differences. Hence, themethod has second order accuracy. In order to facilitatethe calculation of the central differences, the Yee grid isused where the discretization points are staggered by half of the respective step size.
Figure 17 depicts an exam-ple of a Yee grid in one spatial dimension. The mainadvantage of the FDTD scheme is its simplicity. Theimplementation of the method as well as boundary con-ditions or sources is straightforward. Additionally, itcan be executed efficiently in parallel, although the naiveimplementation will not yield the maximum performanceand a more advanced approach must be used.
Themajor drawback is the introduced numerical dispersionwhich can only be avoided by using very fine discretiza-tion sizes. Otherwise, artifacts in the simulation resultscould be the consequence. In the context of MB simu-lations, different values (and value ranges) for the max-imal spatial discretization size ∆ z have been found ade-quate for the FDTD scheme. Namely, λ /20 to λ /100, λ /50, λ /100, and λ /200 have been used, where λ represents the smallest occuring wavelength. The max-imum time step ∆ t is, similarly as in Section VII A 1,determined by the Courant number, which leads for theFDTD scheme to the condition v ∆ t < ∆ z [or v ∆ t ≤ (cid:0) ∆ − x + ∆ − y + ∆ − z (cid:1) − / for three spatial dimensions]. Here, the velocity is obtained from the parameters inEq. (44) as v = ( µ (cid:15) (cid:15) r ) − / . In related literature, choos-ing v ∆ t = ∆ z / was found to be adequate. To reduce the numerical burden, different approachesusing the pseudo-spectral time-domain method havebeen presented.
This method calculates the spa-tial derivatives using the fast Fourier transform in space.As long as Nyquist-Shannon theorem is not violated, themethod is exact in space (and the introduced numeri-cal dispersion minimal). However, the time derivative isstill approximated with finite differences that cause nu-merical error and dispersion. Nevertheless, fewer spatialgrid points are required to achieve reasonable accuracy(for example, the spatial discretization size ∆ z = λ/ has been used ). Thereby, the computational work-load is reduced. These advantages come at the price ofa more complicated implementation. In particular, ab-sorbing boundary conditions must be implemented in or-der to avoid the wrap-around effect. Furthermore, sharpmaterial parameter changes and the implementation ofsources are not trivial anymore.
2. Coupling Electric Field Updates and Density MatrixUpdates
Since the electric field in Maxwell’s equations and thedensity matrix in the Bloch equations depend on eachother, this coupling must be treated appropriately for anynumerical method that solves Maxwell’s equations. Bidé-garay distinguishes between strongly and weakly coupledmethods.
The difference is the discretization of thedensity matrix in time and in relation to the electric field.Strongly coupled methods discretize the density matrixand the electric field at the same time value, weakly cou-pled methods apply a discretization which is staggered (ahalf time step difference between density matrix and elec-7 (a) PC (b) RK (c) MEt z
FIG. 17. Discretization and data dependencies of the finite-difference time-domain (FDTD) method combined with the(a) predictor-corrector (PC), (b) Runge-Kutta (RK), and (c)matrix exponential (ME) method for the density matrix up-dates. The crosses and circles denote the electric and mag-netic fields, respectively, while the squares represent the den-sity matrix. (a) The PC approach updates the electric fieldand the density matrix in parallel. (b) The RK method firstupdates the electric field (solid line) and then the density ma-trix (dashed line), where in the latter both the old and thenew electric field values are used. (c) For the ME scheme,the updates of density matrix and the magnetic field are per-formed in parallel. The electric field discretization is shiftedby a half time step. tric field). In the following, we discuss various approachesto update the density matrix with different forms of cou-pling.
3. Crank-Nicolson Scheme/Predictor-Corrector Method
The pioneering work by Ziolkowski et al. treatsthe Bloch equations with the Crank-Nicolson scheme ˆ ρ n +1 = ˆ ρ n + ∆ t (cid:2) F n +1 (cid:0) ˆ ρ n +1 (cid:1) + F n (ˆ ρ n ) (cid:3) , (147)where n corresponds to the time t n = n ∆ t , ∆ t is the timestep size, and F (ˆ ρ ) = L (ˆ ρ ) + D (ˆ ρ ) represents the righthand side of the Lindblad equation (3). Since this im-plicit scheme requires solving a linear system of equationsat every time step, usually modifications are employed toreduce the numerical load, such as keeping the field at afixed value while advancing the density matrix by a timestep. A widely used variant is based on the predictor-corrector technique, where the update step first initializes ˆ ρ PC = ˆ ρ n , then executes the procedure ˆ ρ PC ← ˆ ρ n + ∆ t F (cid:18)
12 ˆ ρ PC + 12 ˆ ρ n (cid:19) (148)four times, and finally assigns the result to the value ˆ ρ n +1 = ˆ ρ PC . In Fig. 17(a), the coupling of the method to the FDTDscheme is illustrated. It should be noted that this is astrongly coupled method and the electric field is updatedwith the same procedure (of course, F is replaced withthe right hand side of Ampere’s law) and in parallel tothe density matrix update.
4. Runge-Kutta Method
Several research groups use the fourth-orderRunge-Kutta (RK) method to solve the Blochequations.
As illustrated in Fig. 17(b),the method is strongly coupled since electric field anddensity matrix are discretized at the same time steps.The exact procedure is not always described in relatedwork, but can be outlined as follows:
First, the electricfield is updated using the standard FDTD update step.Then, the update of the density matrix using the rule ˆ ρ n +1 = ˆ ρ n + ∆ t ( k + 2 k + 2 k + k ) / (149)follows, where k = F n (ˆ ρ n ) , k = F n +1 / (ˆ ρ n + ∆ t k / , k = F n +1 / (ˆ ρ n +∆ t k / , and k = F n +1 (ˆ ρ n +∆ t k ) . Since the F n contains the electric field E n at time t n = n ∆ t , not only the old and updated field values are re-quired, but also the value at the half time step. The lattercan be approximated by averaging between the old andthe updated field value, i.e., E n +1 / ≈ (cid:0) E n + E n +1 (cid:1) / .
5. Matrix Exponential Methods
The methods of this group aim to solve the Blochequations exactly for one time step. As illustrated inFig. 17(c), the updates of electric field and density ma-trix are weakly coupled, i.e., their updates are performedalternatingly. The density matrix update reads ˆ ρ n +1 / = exp ( F n ∆ t ) ˆ ρ n − / , (150)where F n may depend on the electric field E n and exp( F n ∆ t ) represents the exact solution of the Lindbladequation. After that, the standard FDTD update rulecalculates E n +1 using E n and ˆ ρ n +1 / .If an analytical expression for the solution superop-erator exp( F t ) exists, this method is clearly the mostaccurate one. However, finding such an analytical ex-pression is far from trivial. In fact, the exact form of theexponential depends on the representation. In Liouvilleor coherence vector representation described in SectionIII D, the solution superoperator has the form exp( F t ) ,where F is a matrix. While this is straightforward tosolve, the size of the matrix is in the order N × N for a N × N density matrix. Since the exponential of a N × N matrix would generally need O ( N ) operations, calculat-ing the exponential in Liouville representation requires O ( N ) operations and becomes unfeasible for large N .In regular representation, a solution for the Lindbladequation must be found first. The Strang splitting tech-nique can help here to separate the effects of the Li-ouvillian L and the dissipation superoperator D . Thesolution for the Liouvillian requires the calculation of exp( − i (cid:126) − ˆ Ht ) , where the Hamiltonian ˆ H is a N × N Hermitian matrix. The calculation requires O ( N ) oper-ations, which is still quite intensive.8The Strang splitting introduces an additional error of O (∆ t ) in general. Furthermore, F is generally time de-pendent due to its dependence on the time-varying elec-tric field, in which case the resulting matrix exponentialscontain an integral in the exponent. Commonly, the in-tegral is approximated using the midpoint rule. Thisleads to the conclusion that in reality the accuracy ofmatrix exponential methods is comparable to other ap-proaches. Nevertheless, this group of methods preservescertain matrix properties and despite their limited per-formance they have attracted the focus of many researchgroups.Several techniques have been applied in order to im-prove the performance of matrix exponential methods.The already mentioned Strang splitting has not onlybeen used to allow analytical solutions, but also to sep-arate the time dependent and time independent part of F . This has the advantage that a part ofthe solution can be precalculated and applied at everytime step, while for the remaining part efficient eval-uation techniques exist in some cases. For example,we discovered that the coherence vector representationleads to a real skew-symmetric matrix in the exponen-tial. This expression can be evaluated efficiently usingthe generalized Rodrigues’ formula.
Other techniquesto calculate the matrix exponential have been appliedin related work: An approximation based on the Cay-ley transform,
Magnus expansion via Sylvester’sformula, diagonalization of the matrix, the scal-ing and squaring method as well as a Krylov subspacemethod, and Chebyshev polynomials.
6. Comparison of Numerical Methods for the BlochEquations
As already outlined above, the matrix exponentialmethods are the most computationally expensive ones.In fact, this was confirmed in a detailed investigation, where both the Runge-Kutta and the predictor-correctorimplementation outperformed the matrix exponentialmethod. In this comparison, the predictor-correctormethod demonstrated the best performance.In terms of accuracy, Runge-Kutta methods have thehighest order. However, the accuracy alone is not thecrucial criterion. In particular, it was demonstrated thatthe Crank-Nicolson scheme does not preserve the posi-tivity of the density matrix in the general case (at leastwhen more than two energy levels are considered) andtherefore might yield unrealistic results, e.g., negativepopulations.
Furthermore, it was found that both thepredictor-corrector and Runge-Kutta method yield neg-ative populations in certain cases (e.g., long simulationend time combined with unfortunate choices for the timestep size), while the matrix exponential method preservesthe properties of the density matrix independent of thesimulation settings.
7. Alternative Methods
Besides the full-wave Bloch equations of the formEq. (17), also related formalisms are used to model quan-tum systems interacting with a semiclassical optical field,requiring adapted numerical schemes which are oftencombined with the FDTD method for Maxwell’s equa-tions. For example, MB simulations which replace theBloch equations by an equivalent evolution equationfor the polarization vector, Eq. (72), require modifiedschemes adapted to the second-order differential form ofEq. (72).
If the dissipation term in the Bloch equations canbe neglected and the quantum system is in a purestate, the time evolution can be described in a simpli-fied manner with the time dependent Schrödinger equa-tion, Eq. (1), for which suitable numerical schemes havebeen developed.
In analogy to the MB equations,the Schrödinger and Maxwell’s equations can be com-bined to the Maxwell-Schrödinger approach, which is forexample used to model nanoelectronic systems, orto describe the interaction of atoms with intense laserfields.
As for the MB equations, such a coupledsimulation complicates the numerical treatment, and var-ious numerical schemes have been developed, e.g., com-bining FDTD or transmission line matrix simulations ofMaxwell’s equations with a spatial grid representation oreigenstate expansion of the wavefunction.
Inthis context, a recent interest has been on algorithmspreserving the symplectic structure of the Maxwell-Schrödinger equations, thus ensuring energy conservationof the coupled system.
VIII. INCLUSION OF FURTHER EFFECTSA. Local-Field Correction
In principle, the current/polarization contribution ofan individual quantum system at a given position can bedirectly represented in Maxwell’s equations by a pointsource, without using ensemble averaging as in Sec-tion IV A. However, the complexity of such an approachincreases significantly with the number of quantum sys-tems to be included.
Moreover, care must be takenthat the field which drives the quantum system does notcontain the divergent self-field contribution of the systemitself, which further adds to the numerical load.
Analternative approach, which is especially suitable for alarge ensemble of quantum systems as considered in thispaper, is based on macroscopic MB equations. Ratherthan setting up Bloch equations for each of the quan-tum systems, the ensemble is here modeled by repre-sentative density matrices ρ ij ( x , t ) distributed over thedevice volume, e.g., placed on the spatial grid points,where x is the macroscopic position coordinate. Like-wise, Maxwell’s equations then contain the macroscopiccurrent densities (see Section IV A) and fields, defined as9 FIG. 18. Illustration of dipoles arranged in a cubic lattice. ensemble averages over the individual microscopic con-tributions. From the macroscopic electric field E whichis averaged over local variations associated with the indi-vidual dipoles, the local microscopic field which interactswith a given physical quantum system can be determinedbased on a compact Clausius-Mosotti type model. Thislocal-field correction can for example lead to frequencyshifts, and becomes relevant for densely spaced quantumsystems as discussed in Section VIII A 1. A relatedeffect emerges in tightly localized artificial quantum sys-tems consisting of multiple semiconductor atoms, see Sec-tion VIII A 2.
It has been pointed out that local-fieldeffects can be exploited as an additional design degree offreedom in nanostructures.
1. Near-Dipole-Dipole Effects in Dense Media
The macroscopic field E comprises contributions ofexternal sources as well as an internal contribution E p due to the induced dipoles in the material, which is re-lated to P q . In the following, we consider a mediumsuch as a gas or crystal lattice which consists of a densecollection of atoms, molecules or other quantum sys-tems, as illustrated in Fig. 18. The local field E L atthe position of the considered quantum system is deter-mined by replacing the volume-averaged field E p withthe microscopic contribution E near due to the nearbydipole moments, E L = E − E p + E near . Based on amicroscopically large, but macroscopically small probevolume, which is conveniently chosen to be spherical,the macroscopic polarization contribution is obtainedas E p = − P q / (3 (cid:15) ) . On the other hand, it can be shown that for dipoles arranged in a cubic lattice as illus-trated in Fig. 18, E near vanishes for a particle on a latticesite, where the particle’s self-field is not included. This yields E L = E + P q / (3 (cid:15) ) , which is also approx-imately true for other reasonably isotropic media andcompletely random arrangements, such as amorphousmedia or gases. Retardation effects are here neg-ligible since the probe volume diameter is assumed tobe much smaller than the optical wavelength.
Fortwo-level systems where the polarization is given byEq. (71), the local-field corrections can thus be includedin the Maxwell-Bloch equations by formally substituting
Ω = d E / (cid:126) with Ω L = d E L / (cid:126) = Ω + 2 ω L (cid:60) { ρ } (151)in Eq. (67), where the static Lorentz shift ω L = d n / (3 (cid:126) (cid:15) ) has the dimension of frequency. Here we have assumeda real-valued d for simplicity. Applying the RWA, weobtain Eq. (73) where we have to substitute ˆΩ by ˆΩ L = ˆΩ + 2 ω L η , (152)which changes Eq. (73a) but not Eq. (73b) since there thelocal-field correction term cancels out. If the quantum systems are embedded in a hostmedium, such as dopant ions in a crystal, they inter-act not only with particles of the same species, but alsowith those of the host material, and the local-field correc-tion must be suitably extended. Accordingly, above ap-proach has been generalized to a dense collection of two-level atoms embedded in a linear, potentially dispersiveand absorptive host medium, and to multicomponentmedia in general.
Furthermore, above concept can bestraightforwardly extended to more than two levels.
We note that above correction to the MB equationshas mainly been considered for ensembles of atoms ormolecules, where typically much higher number densi-ties are obtained than for artificial systems such as quan-tum dots. This model has enabled the analytical andnumerical investigation of numerous effects, suchas solitonic and ultrashort pulse propagation or opticalswitching. For artificial, tightly localized quantum sys-tems, local-field effects are typically governed by the de-polarization field, as discussed below.
2. Depolarization Field in Tightly Localized QuantumSystems
Artificial quantum systems, such as quantum dots, arenanostructures built from a larger number of semicon-ductor atoms. Thus, the quantum systems only ”feel”an averaged polarization contribution of the individualatoms, which can be described by the background dielec-tric constant of the host material (cid:15) r . Here, the local-field0 Applied fieldDepolarization field
FIG. 19. Illustration of a localized polarizable object in anelectric field. correction accounts for the deviation of the mesoscopicaverage field inside the localized polarizable nanostruc-ture from the macroscopic average field E in the en-tire composite material, which enters the macroscopicMaxwell equations. In the following, we assume that (cid:15) r is real-valued and frequency independent, and is identi-cal for the nanostructure and the surrounding material.As illustrated in Fig. 19, the field E generates a polar-ization inside the polarizable object, and the uncompen-sated surface dipoles give rise to bound surface chargeswhich induce an electric field in the object, the so-calleddepolarization field E d , thus altering the local field insideto E L = E + E d . We assume that the object’s dimen-sions are much smaller than the wavelength of the ex-citing field. The (spatially averaged) depolarization fieldis then obtained as E d = − NP q / ( (cid:15) (cid:15) r ) , where N is thedepolarization tensor which only depends on the geome-try of the polarizable object. For a localized two-levelsystem of volume V s , the (volume-averaged) polarizationis given by P q = d ( ρ + ρ ) /V s . Proceeding in thesame way as in Section VIII A 1, we include the local-fieldcorrection into the MB equations by using Eq. (151) orEq. (152), where we now have ω L = − d L d / ( (cid:126) (cid:15) (cid:15) r V s ) . The depolarization factor L d = d Nd /d with ≤ L d ≤ accounts for the anisotropy of the object, suchas an ellipsoidal quantum dot, and becomes L d = 1 / for spherical geometries. The resulting equations havebeen used for both analytical and numerical studies of local field effects in quantum dots.
B. Inhomogeneous Broadening
Homogeneous broadening is naturally considered inthe Bloch equations. This can best be seen from thesteady state solution in RWA, Eq. (128), yielding aLorentzian lineshape with the width given by the dephas-ing rate of the corresponding transition, see Eqs. (129)and (132) as well as Fig. 13(b). In addition, inhomoge-neous broadening arises if the optically active mediumconsists of quantum systems with slightly different reso-nance frequencies, as for example frequently arisesin ensembles of quantum dots due to size fluctuations.Another example is Doppler broadening in a gas, caused by the Doppler shift due to the thermal motion ofthe atoms or molecules. For a given transition betweentwo states | i (cid:105) and | j (cid:105) , the distribution of the resonancefrequency ω ij is commonly described by a distributionfunction g ij ( ω ij − ω ij ) with (cid:82) ∞−∞ g ij ( ω ) d ω = 1 , where ω ij − ω ij is the deviation from the average resonance fre-quency ω ij . For thermal Doppler broadening, g ij is givenby a Gaussian distribution g ij ( ω ij − ω ij ) = 1 √ πσ ij exp (cid:34) − ( ω ij − ω ij ) σ ij (cid:35) . (153)In this case, the standard deviation becomes σ ij =( ω ij /c ) ( k B T /m ) / , where k B is the Boltzmann con-stant, T indicates the temperature, and m is the massof the atom or molecule. The Gaussian distribution,Eq. (153), is also frequently used as a generic model ifthe distribution of resonance frequencies is not exactlyknown, e.g., to describe above mentioned inhomogeneousbroadening in ensembles of quantum dots due to sizefluctuations.
If the individual quantum systems con-tain more than one relevant optical transition with dis-tributed resonance frequencies, in principle joint distri-bution functions have to be used.Numerically, the full-wave or RWA Bloch equations,Eq. (17) or (43), have to be solved separately for eachpossible value of the resonance frequency (or each possi-ble combination of resonance frequency values if the indi-vidual quantum systems contain more than one relevantoptical transition).
This requires discretizing the dis-tribution function into a finite number of N inh bins withresonance frequencies ω sij , s = 1 ..N inh . Each of thesebins is represented by a corresponding quantum systemsubensemble with density matrix ρ sij , where the fractionof carriers p s is proportional to the weight of that binand (cid:80) s p s = 1 . The polarization current density is thenobtained from a generalized version of Eq. (48), J q = n (cid:88) i,j d ji (cid:88) s p s ∂ t ρ sij . (154)In SVAA, the polarization amplitude is given by a gen-1eralized form of Eq. (64), P = 2 n (cid:88) ω ij > d ji (cid:88) s p s η sij . (155)For certain broadening mechanisms, such as fluctuationsin quantum dot size, also d ij can in principle vary whichcould be considered by introducing a quantity d sij inEqs. (154) and (155) in analogy to ω sij , but this effectis usually neglected.In certain cases, inhomogeneous broadening can alsobe considered in analytical solutions of the MB equationsbased on the RWA. These usually invoke the factoriza-tion ansatz, which assumes that non-resonant systemswith a finite frequency detuning ∆ = ω c − ω (cid:54) = 0 essen-tially respond in the same way to the optical field as theresonant ones, apart from a detuning dependent ampli-tude F (∆) . In Section VI A 2, this approach has beendemonstrated in the context of self-induced transparency.
C. Noise
Noise in optoelectronic devices arises for example fromspontaneous emission and from processes in the semi-conductor host, such as lattice vibrations. Noise andfluctuations can generally be included into the semiclas-sical MB equations by adding stochastic terms.
Numerically, the stochastic terms are typically imple-mented by using a pseudorandom number generator toobtain uncorrelated, Gaussian distributed random num-bers for every grid point.
The MB equations, com-plemented by these additional stochastic terms, are thennumerically solved as discussed in Section VII, for exam-ple with an FDTD-based approach.
The stochas-tic terms are systematically obtained from quantumLangevin equations, which are then represented byequivalent stochastic c-number equations, i.e.,evolution equations for operator expectation values withadditional stochastic terms.Spontaneous emission obviously plays an importantrole in optoelectronic devices. While the resulting recom-bination can simply be included by nonlinear rate termsfor the carrier occupations in Eq. (17b) or (43b), the noise contribution is not included in the MB modeldue to its semiclassical nature. This effect can howeverbe considered in terms of a Gaussian white noise sourcein the optical propagation equation.
In a dif-ferent model, also dipole fluctuations are included byadding Langevin noise terms not only to the propagationequation, but also to Eq. (17a) for the off-diagonal den-sity matrix elements.
By virtue of the fluctuation-dissipation theorem, a decay of populations, coherencesor the optical field is generally accompanied by fluc-tuations, and an MB equation model which includessuch decay-induced fluctuations has been presented.
Furthermore, an extension of the stochastic c-numberapproach to incorporate nonclassical effects has beendiscussed.
IX. APPLICATION TO OPTOELECTRONIC DEVICESA. Bulk and Quantum Well Interband OptoelectronicDevices
For interband optoelectronic devices based on semi-conductor bulk or quantum-well media, the conductionand valence band states are given by | c , k (cid:105) and | v , k (cid:105) ,respectively. Here, k is the three-dimensional crystalwavevector in bulk media or the two-dimensional in-planewavevector for quantum well structures, see Section III C.In the following, we define ρ c , k and (1 − ρ v , k ) as the elec-tron occupation probability of a conduction band state | c , k (cid:105) and a valence band state | v , k (cid:105) , respectively, i.e., ρ v , k corresponds to the hole occupation probability ofa valence band state. Restricting ourselves to directbandgap semiconductors in a two-band approximationand using that the optical transitions are k conserving,as can for quantum wells be seen from Eq. (24), theBloch equations, Eq. (17), become ∂ t ρ cv , k = − i ω cv , k ρ cv , k − iΩ k ( ρ c , k + ρ v , k − ∂ t ρ cv , k ] col , (156a) ∂ t ρ α, k = i (cid:0) Ω k ρ ∗ cv , k − Ω ∗ k ρ cv , k (cid:1) + [ ∂ t ρ α, k ] col , (156b)with α = c , v . Here, the sum over k implicitly also in-cludes summation over the two possible spin orientations.This equation applies to quantum wells, where the carri-ers are treated as a two-dimensional gas, as well as bulkmedia. The dissipation processes are here included bygeneral collision terms [ ∂ t ρ cv , k ] col and [ ∂ t ρ α, k ] col . Thesecan be modeled on a microscopic level, in particularaccounting for carrier-carrier and carrier-phonon scat-tering, or under certain approximations by relaxationrate terms similar to those in Eq. (17). Many-bodyCoulomb interactions can be taken into account basedon the Hartree-Fock approximation, which results in theso-called semiconductor Bloch equations, which have theform of Eq. (156) but feature renormalized transition andRabi frequencies ω cv , k = 1 (cid:126) E cv , k − (cid:126) (cid:88) q (cid:54) = k V | k − q | ( ρ c , q + ρ v , q ) , (157a) Ω k = d cv , k E (cid:126) + 1 (cid:126) (cid:88) q (cid:54) = k V | k − q | ρ cv , q , (157b)leading to a coupling of the states with differ-ent k through the time dependent renormalizationterms. In Eq. (157a), E cv , k is the energy of freeelectron-hole pairs, which can in a simple model be de-scribed as E cv , k = 12 (cid:126) k (cid:18) m ∗ e + 1 m ∗ h (cid:19) + E g . (158)Here, E g denotes the band gap energy, or for quantumwells the energy difference between the electron and holestate energies. Furthermore, m ∗ e and m ∗ h are the electron2and hole effective masses, assuming a parabolic disper-sion relation near the conduction and valence band edges,respectively. V q is the Coulomb potential in Fourier rep-resentation, which is in a bulk with probe volume V p andbackground permittivity (cid:15) r given by V q = e / (cid:0) (cid:15) (cid:15) r V p q (cid:1) ,and in a quantum well structure with in-plane cross sec-tion S by V q = e / (2 (cid:15) (cid:15) r qS ) . Screening effects, whichresult from the response of the other carriers and weakenthe potential, are then incorporated as corrections tothe Hartree-Fock equations in the form of a modified V q . Summing in Eq. (48) over the initial and finalstates | i (cid:105) = | c , k (cid:105) and | j (cid:105) = | v , q (cid:105) where we consider thatthe total number density of electron-hole pairs is givenby n cv ( t ) = V − (cid:88) k ρ c , k ( t ) = V − (cid:88) k ρ v , k ( t ) , (159)and using the k conservation of optical transitions, d ij = d cv , k δ k , q , we obtain for the polarization term ∂ t P q = V − (cid:88) k ( d cv , k ∂ t ρ vc , k + c.c. ) . (160)With Eq. (160), Eqs. (156) and (157) can be coupled tothe Maxwell equations, Eq. (44), resulting in the semi-conductor Maxwell-Bloch equations. These equationshave been extensively used for the simulation of semi-conductor lasers and related devices. Further-more, they have been adapted to the modeling of quan-tum wire structures, as well as graphene andcarbon nanotubes.
As stated in Section I, a further discussion of thismodel is beyond the scope of this paper. Rather, wewill focus here on approaches which reduce above two-band model to macroscopic two- or N -level Bloch equa-tions. As a first step, the collision terms in Eq. (156)are modeled by relaxation rate terms similar to those inEq. (17). Extending the rate equation model Eq. (9)to states | α, k (cid:105) results in a Boltzmann-type collision termfor the populations, which can in consideration of Pauliblocking be written as [ ∂ t ρ α, k ] col = − ρ α, k (cid:0) W ααβ, k + W ααα, k (cid:1) + (1 − ρ α, k ) (cid:16) W βαβ, k + W βαα, k (cid:17) . (161)Here, β (cid:54) = α , i.e., for the conduction band collision term( α = c ) we have β = v , while β = c for α = v . The Boltz-mann rates are related to the electron transition ratesof Section II B 1 by W c αα (cid:48) , k = (cid:80) k (cid:48) r α, k → α (cid:48) , k (cid:48) (1 − f α (cid:48) , k (cid:48) ) and W v αα (cid:48) , k = (cid:80) k (cid:48) r α (cid:48) , k (cid:48) → α, k f α (cid:48) , k (cid:48) , where f α, k de-notes the electron occupation probability, i.e., f c , k = ρ c , k and f v , k = 1 − ρ v , k . For example, spontaneousemission and carrier-phonon scattering can be modeledby Eq. (161) with adequately chosen transition rates r α, k → α (cid:48) , k (cid:48) , while the inclusion of carrier-carrier inter-actions beyond Hartree-Fock effectively requires rateswhich themselves depend on the carrier distribution. Furthermore modeling the dephasing in analogy toEq. (12), we obtain [ ∂ t ρ cv , k ] col = − γ cv , k ρ cv , k , (162a) [ ∂ t ρ α, k ] col = Γ βα, k − r αβ, k ρ α, k + Γ αα, k − r αα, k ρ α, k . (162b)In Eq. (162a), γ cv , k indicates the dephasing rate. Re-arranging the contributions in Eq. (161), Eq. (162b) isobtained, where the first two terms on the right handside with β (cid:54) = α describe interband processes, while theother two terms model the intraband transitions. Here, r αα (cid:48) , k = W ααα (cid:48) , k denotes the interband (for α (cid:48) = β ) or in-traband (for α (cid:48) = α ) recombination rate due to nonradia-tive transitions and spontaneous emission. Furthermore, Γ α (cid:48) α, k = (1 − ρ α, k ) W βαα (cid:48) , k describes the filling of state | α, k (cid:105) , and can for α (cid:48) = β be interpreted as a pump rate,where carriers are induced for example by an injectioncurrent or optical pumping. Summation of Eq. (161) orEq. (162b) over k yields (cid:80) k [ ∂ t ρ c , k ] col = (cid:80) k [ ∂ t ρ v , k ] col ,as expected from Eq. (159). In more detail, for the intra-band contributions we have (cid:88) k (Γ αα, k − r αα, k ρ α, k ) = 0 , (163)while the interband terms fulfill (cid:88) k Γ vc , k = (cid:88) k Γ cv , k , (cid:88) k r cv , k ρ c , k = (cid:88) k r vc , k ρ v , k . (164)To obtain compact two-level Bloch equations, we as-sume a k independent dipole matrix element d cv , k = d in Eq. (157b) and ignore the renormalization contribu-tion, yielding the usual definition Ω k = (cid:126) − d E . Sum-ming Eq. (156b) over k yields with Eqs. (159), (162b),(163) and (164) and P cv = V − (cid:80) k ρ cv , k ∂ t n cv = i (Ω P ∗ cv − Ω ∗ P cv ) + Γ − γ n cv . (165)For electrical pumping, the injection rate Γ = (cid:80) k Γ βα, k (with α = c , β = v or α = v , β = c ) canbe modeled as Γ = ηV − I/e , where η and I denote theinjection efficiency and current, respectively. The recom-bination rate γ , which includes nonradiative and spon-taneous transitions, is obtained by averaging over thecarrier distribution, γ = (cid:80) k r αβ, k ρ α, k /n cv . Proceedingin a similar manner for Eq. (156a) by neglecting the k dependence of ω cv , k is not feasible, due to the problemsarising from the summation over ( ρ c , k + ρ v , k − . Var-ious strategies have been developed to circumvent thisproblem. Here we follow the approach by Yao etal., formulated in the framework of the RWA. Thus,we start by introducing the slowly varying field envelope E and the transformed off-diagonal elements defined inEqs. (41) and (42), respectively. Furthermore assuming a3 k independent Rabi frequency Ω = (cid:126) − d cv E , Eq. (156a)yields with Eq. (162a) ∂ t η cv , k = (i∆ k − γ cv , k ) η cv , k − i2 Ω ( ρ c , k + ρ v , k − , (166)where ∆ k = ω c − ω cv , k with ω cv , k given by Eq. (157a).In the framework of semi-phenomenological macroscopicMB equation models, the renormalization term is oftenneglected. The k dependence can be modeled withEq. (158) or a more sophisticated description. DividingEq. (166) by (i∆ k − γ cv , k ) , summing over k and defining p cv = V − (cid:80) k η cv , k , we obtain ∂ t p cv = − i2 τ τ Ω n cv − τ − p cv , (167)where we have furthermore used Eq. (159) and introducedthe complex parameters τ and τ with V − (cid:88) k ρ c , k + ρ v , k − γ cv , k − i∆ k = τ n cv ,V − (cid:88) k η cv , k γ cv , k − i∆ k = τ p cv . (168)Equation (167) and Eq. (165) in RWA, ∂ t n cv = i2 (Ω p ∗ cv − Ω ∗ p cv ) + Γ − γ n cv , (169)constitute macroscopic Bloch equations for interbandtransitions in bulk semiconductor and quantum wellsystems. In order to obtain a Maxwell-Bloch model,Eqs. (167) and (169) can be coupled to the optical prop-agation equation in SVAA, Eq. (62), where the RWA po-larization term is obtained from Eq. (75) as P = 2 d cv p cv .The parameters τ and τ introduced in Eq. (168) canfor example be evaluated numerically, or by fitting toexperimental data. In this context, it has been foundthat τ and τ can be treated as independent of the opti-cal intensity, but that especially τ shows a pronounceddependence on n cv which should be considered in themodel, and also allows a phenomenological reintroduc-tion of renormalization effects. B. Quantum Well Intersubband Devices
Intersubband devices, such as QCLs, quantum cas-cade detectors and quantum well infrared photodetec-tors (QWIPs), commonly utilize optical intersubbandtransitions between quantized energy levels in the con-duction band Γ valley of a multiple quantum well struc-ture. The Maxwell-Bloch model has been extensively ap-plied to such devices, especially for the dynamic mod-eling of QCLs. The quantized states | ψ i , k (cid:105) , also re-ferred to as subbands, are characterized by their wave-function ψ i ( x ) where i is the subband index, and in-plane wavevector k = [ k y , k z ] T . These states are com-monly found by solving the one-dimensional effective mass Schrödinger equation, obtained from inserting theansatz Eq. (19) into Eq. (20), for the quantum well po-tential V ( z ) . As mentioned in Section III C 1, bandbending due to space charge effects is usually consid-ered by solving the coupled Schrödinger-Poisson equa-tion system, and also nonparabolicity effects, whichplay a role especially in mid-infrared devices, can beincluded. These calculations yield the eigenenergies E i and wavefunctions ψ i , and thus the transition frequen-cies ω ij = (cid:126) − ( E i − E j ) and the dipole matrix elements,Eq. (25), as input for the Bloch equations. We againchoose the semiconductor Bloch equations as a startingpoint, with a form analogous to Eq. (156). Due to thetypically low doping levels of QCLs, the Hartree–Fockrenormalization effects in Eq. (157) have been found tobe relatively small, and also Pauli blocking only playsa secondary role. Furthermore, we can assume k inde-pendent transition frequencies ω ij at least for terahertzQCLs, where the energetic level spacings are smaller thanin mid-infrared QCLs and thus the subbands have nearlyparallel dispersion relationships. Under these assump-tions, summing over k yields Bloch equations of the formEq. (17), where we have used that the dipole matrix ele-ment is k conserving, and introduced intersubband scat-tering rates r j → i which are related to the generally k de-pendent rates r j, q → i, k by r j → i ρ jj = (cid:80) k , q r j, q → i, k ρ jj, q .Assuming either moderate temporal variations of theintrasubband electron distribution ρ jj, k or a moderate k dependence of the rates, the r j → i are approximatelygiven by r j → i = (cid:88) k , q r j, q → i, k ρ jj, q / (cid:88) k ρ jj, k . (170)Here, ρ jj, k describes the steady state electron distri-bution in the subband. Notably, the intrasubbandelectron distributions in QCLs can often be reason-ably well approximated by Fermi–Dirac or Maxwell-Boltzmann distributions, parametrized by subband elec-tron temperatures which can significantly exceed the lat-tice temperatures.
By contrast, the off-diagonaldensity matrix elements generally vary strongly withtime, and no clearly defined concept exists how the k averaging should be performed to obtain an effectivedephasing rate γ ij from a relaxation term of the formEq. (162a), [ ∂ t ρ ij, k ] col = − γ ij, k ρ ij, k . This especially mat-ters if the ratio ρ ii, k /ρ jj, k has a strong k dependence,as is the case for significantly different subband electrontemperatures or highly non-thermal distributions. Of-ten, an average over the population inversion of the in-volved subbands is applied, γ ij = (cid:88) k γ ij, k (cid:12)(cid:12) ρ ii, k − ρ jj, k (cid:12)(cid:12) / (cid:88) k (cid:12)(cid:12) ρ ii, k − ρ jj, k (cid:12)(cid:12) , (171)and a comparison to an alternative way of averaginghas yielded similar results for terahertz QCLs. Apartfrom very few exceptions based on the full-wave MBequations, usually the MB equations in RWA4and SVAA, Eqs. (124) and (43), are used for the model-ing of QCLs and related devices. Also, apart from somecases including multiple subbands, typi-cally a two-level model is employed. In the case of mid-infrared QCLs, where nonparabolicity effects play a moreimportant role, an approach similar to Eq. (167) can beenvisioned. The transition and dephasing rates are usu-ally empirically chosen, or extracted from fits to experi-mental data. Alternatively, they can be calculated fromEqs. (170) and (171) based on the Hamiltonians of therelevant scattering mechanisms, such as electron-electroninteractions, scattering with acoustic and longitudinaloptical phonons, as well as impurity, interface rough-ness and alloy scattering.
Here, the use of dissipa-tion rates derived from steady-state models is consistentwith the Markovian and time-homogeneous character ofthe Lindblad dissipator, which provides the basis for theBloch equations. The corresponding scattering rates aretypically evaluated based on Fermi’s golden rule, whilethe associated pure dephasing rates can be obtained fromAndo’s model.
For QCLs, the MB equations have primarily beenused to study ultrashort pulse generation by mode-locking and the closely relatedformation of coherent instabilities, as well asthe generation of frequency combs.
In de-tail, it has been found that coherent multimode instabili-ties result in the emergence of sidebands around the orig-inal longitudinal mode, giving rise to broadband multi-mode operation.
Furthermore, active mode-lockinghas been investigated where short pulses are generatedby modulating the laser current at the cavity roundtripfrequency, yielding good agreement between simulationsand measurements.
Also the possibility of real-izing passive mode-locking in QCLs has been theoret-ically explored.
Here, pulse formation isobtained by adding saturable absorption regions, whereSIT mode-locking, discussed in Section VI A 2, consti-tutes a special variant. Besides, frequency comb opera-tion has been studied, where an equidistant line spectrumis generated, which serves as a ruler in the frequency do-main for spectroscopic and sensing applications. Here,a perturbative treatment of the MB equations imposinga comb-like spectrum has been employed, as wellas full numerical simulations.
In most of aboveworks, spatial hole burning has been considered based onEqs. (123) and (124), as it considerably affects the QCLdynamics. In addition, various other effects have beenimplemented which can play an important role for mode-locked and frequency comb operation in QCLs, such astunneling across thick barriers and group velocitydispersion due to the waveguide and bulk semiconductormaterial.
For optical excitation on veryshort timescales, memory effects become important andthe presuppositions of the Lindblad approach are too re-strictive, requiring the use of more complex models suchas the quantum-kinetic schemes.
As mentioned above, MB simulations have for example P o w e r s pe c t r u m [ a r b . u .] FIG. 20. Optical comb spectrum of a terahertz QCL, as ob-tained from (a) simulation and (b) experiment. P o w e r [ a r b . u .] Time (normalized to roundtrip time)(a)(b)
High freq. lobe Low freq. lobe
FIG. 21. Instantaneous optical power in the high and lowfrequency lobe as obtained from (a) simulation and (b) ex-periment. been used to model frequency comb operation of QCLs,identifying four-wave mixing as the primary comb form-ing mechanism and explaining experimentally observedfeatures.
In Fig. 20, a comparison betweensimulation and experiment is presented for the powerspectrum, generated by a THz QCL for frequency combgeneration. Good agreement is found; in particular, asplitting of the comb spectrum into a high and a lowfrequency lobe is observed in both simulation and ex-periment. Also the simulated temporal dynamics agreeswell with experiment. In Fig. 21, the simulated and mea-sured instantaneous optical power in the high and low fre-quency lobe of the comb is shown. Again, good agree-ment between theory and experiment is obtained, con-firming the validity of the MB model. In particular, thetemporal switching behavior between the two lobes is re-produced in the simulation.5 (a)(b)
ContactSubstrateWetting layer Light
FIG. 22. Schematic illustration of (a) wetting layer with QDsand (b) QD laser.
C. Quantum Dot Devices
Due to the strong carrier localization and discrete en-ergy spectrum resulting from the three-dimensional con-finement in QDs, they enable lasers and laser ampli-fiers with excellent gain, threshold, temperature, and dy-namic characteristics.
While these devices rely on in-terband optical transitions, also the possibility has beenstudied to exploit intraband transitions similarly as inQCLs to obtain lasing in the mid-infrared or terahertzregime.
Furthermore, intraband transitions be-tween bound electron or hole states (or from bound tocontinuum states) have been employed for quantum dotinfrared photodetectors.
In contrast to bulk semiconductors and quantum wellor wire structures which feature a continuum of statesdue to the free carrier motion in at least one dimen-sion, the QD possesses a discrete set of energy eigen-states. Thus, the application of phenomenological modelsbased on generic discrete-level MB equations appears tobe especially justified for QD systems. Phenomenologicaltwo-level MB equations have been employed for a largerange of applications based on QDs. This includes stud-ies of the spatiotemporal dynamics and SIT mode-locking in QD lasers, FDTD-based MB simulationsof QD photonic-crystal-cavity lasers, and QDs coupledto a nanoparticle or cavity. Furthermore, three-level MB equations have been used, for example to studyEIT, soliton propagation or all-optical switching in QD structures, and also four-level models have beendeveloped.
Optoelectronic applications employing large ensemblesof QDs are often fabricated utilizing self-assembly of theQDs on top of an initial quasi-two-dimensional semicon-
Wetting layerWetting layer Q uan t u m do t r c,1 r c,2 r v,1 r v,2 B andgapene r g y w c w v w cv , w cv , J p FIG. 23. Schematic energy diagram of QD and wetting layer.The big arrows represent coherent light-matter interaction,the dotted arrows indicate nonradiative intraband transitions,and the wavy arrows represent spontaneous emission. ductor layer, which is referred to as wetting layer, assketched in Fig. 22(a). The resulting structure is subse-quently covered by another layer of suitable semiconduc-tor material. The wetting layer effectively forms a quan-tum well, which serves as a reservoir for the carriers. Thethus obtained QD layer forms the basis of various devicessuch as QD lasers [see Fig. 22(b)], where commonly mul-tiple layers are stacked on top of each other to increasethe optical gain.In Fig. 23, a schematic energy diagram of the wet-ting layer and a QD is shown. A description of theQD dynamics based on the semiconductor Bloch equa-tions, Eqs. (156) and (157), features renormalized tran-sition and Rabi frequencies due to many-body Coulombinteractions.
These renormalization effects are of-ten neglected so that the conventional Bloch equations,Eq. (17), can be used as a starting point, which are fre-quently supplemented by a detailed model for Coulombscattering and other scattering mechanisms in the dissi-pation term.
In the following, the QD conduction( α = c ) and valence ( α = v ) band states are labeled byan index i . Furthermore, as described in Section VIII B,variations in QD size result in distributed resonance fre-quencies, and the associated inhomogeneous broaden-ing is included by dividing the QDs in correspondingsubensembles s containing a fraction p s of QDs. TheBloch equations, Eq. (17), thus have to be adapted byreplacing the density matrix elements ρ ij with ρ sαα (cid:48) ,ij forpairs of states | α, i (cid:105) and | α (cid:48) , j (cid:105) , where we write for com-pactness ρ sαα,ij = ρ sα,ij , ρ sαα (cid:48) ,ii = ρ sαα (cid:48) ,i and ρ sαα,ii = ρ sα,i .For the dipole matrix element vectors d ij , frequencies ω ij and dephasing rates γ ij , we proceed analogously. Simi-larly as for Eq. (156), the ρ s v ,i are taken as the hole oc-6cupation probability of the i th QD valence band level,i.e., matrix elements ρ ii in Eq. (17) referring to the elec-tron occupation probabilities of QD valence band stateshave to be substituted by (cid:0) − ρ s v ,i (cid:1) . Apart from the co-herent light-matter interaction, incoherent carrier transi-tions in QD systems mainly occur due to carrier-carrierscattering which gives rise to Auger-type processes, aswell as carrier-phonon interactions and spontaneous pho-ton emission. Rate equation terms of the formEq. (9) with phenomenologically chosen parameters arefrequently used to model incoherent transitions in QDsystems.
For a more detailed model-ing, it must be taken into account that important dissi-pative processes in QDs depend on the occupations of twoor more states, and that Pauli blocking is not includedin Eq. (9). This can be addressed by using an empiri-cal nonlinear rate equation model, or based on amicroscopic treatment.
For high pump cur-rents, Auger processes, where two carriers scatter fromtheir respective initial to final levels, involving QD andwetting layer states, constitute the dominant scatteringprocess.
This includes scattering of two electrons orholes, as well as mixed processes involving an electrontransition in the conduction band and a hole transition inthe valence band. The associated change of the occupa-tion ρ sα,i is in the following represented by a generic intra-band collision term (cid:2) ∂ t ρ sα,i (cid:3) intra , which can be gener-alized to also include other scattering-induced intrabandcarrier transitions, e.g., due to electron-phonon interac-tions. Additionally, spontaneous electron-hole recombi-nation is typically taken into account as an importantinterband process, which depends on the occupations ofthe initial and the final state. Within this model, the dis-sipation terms in Eq. (17b) are substituted by the moregeneral ansatz for incoherent processes (cid:2) ∂ t ρ sα,i (cid:3) inc = (cid:2) ∂ t ρ sα,i (cid:3) intra − (cid:88) j A s cv ,ij ρ sα,i ρ sβ (cid:54) = α,j , (172)with the spontaneous recombination coefficient A s cv ,ij .The Lindblad dephasing rate approach of the formEq. (12), as also used in Eq. (17a), has been arguedto be generally well suited to model dephasing inQDs. The dephasing rates can be calculated basedon microscopic models for carrier-carrier and carrier-phonon scattering, or are phenomenologicallychosen.
Considering that (cid:2) ∂ t ρ sα,i (cid:3) intra contains intra-QD tran-sitions as well as carrier exchange between the wettinglayers and QD states, this term not only depends on theoccupations of the QD states involved, but also on thecarrier densities in the wetting layers. Thus, for a closedcarrier transport model, the Bloch equations have to beextended by equations for the wetting layers, which canbe modeled by ∂ t w α = J p e − A cv w c w v − n (cid:88) i,s p s (cid:2) ∂ t ρ sα,i (cid:3) intra . (173) Here, w α denotes the overall carrier sheet densities inthe wetting layers, i.e., w c is the total number of con-duction band electrons in all wetting layers divided bythe area S of a wetting layer, and w v is defined analo-gously for the valence band holes. J p denotes the electricpump current density. Furthermore, n is the overallQD sheet density, and the factor accounts for the spindegeneracy of the QD states. A cv in Eq. (173) is the ratecoefficient for spontaneous band-band recombination inthe wetting layers. Sometimes the carrier injection fromthe bulk to the quantum well wetting layers is modeledby additional equations. For self-assembled quantumdash structures, the wetting layers can be considered inan analogous manner.
The extended Bloch equations, Eqs. (17), (172) and(173), are then coupled to Maxwell’s equations, Eq. (44),by the polarization current density for inhomogeneouslybroadened media given in Eq. (154), J q = n d g (cid:88) α,α (cid:48) (cid:88) i,j d α (cid:48) α,ji (cid:88) s p s (1 − δ α v δ α (cid:48) v δ ij ) ∂ t ρ sαα (cid:48) ,ij , (174)with the thickness of the gain medium d g . Here, the term (1 − ... ) compensates for the fact that the electron occu-pation probabilities of QD valence band states are in ourdensity matrix convention given by (cid:0) − ρ s v ,i (cid:1) . QD lasersand amplifiers usually operate in TE mode, due to thecharacter of the eigenstates for the QD shapes and strainsobtained with the widely employed Stranski–Krastanovgrowth mode. A corresponding one-dimensional MBmodel can be obtained by combining the extended Blochequations with Maxwell-type equations for TE operation,Eq. (92), where the finite overlap of the QD active regionwith the mode profile is considered by the field confine-ment factor, Eq. (88). The contribution of the sponta-neous emission processes in Eq. (172) to the optical fieldis in most cases neglected, but can be considered as dis-cussed in Section VIII C.For interband QD devices, optical intersubband tran-sitions can be neglected. The QD interband dipole ma-trix elements are given by Eq. (22). As discussed in Sec-tion III C 2, for the uppermost valence and lowest con-duction band states, in good approximation only opticalinterband transitions between states with equal quantumnumbers are allowed, and the corresponding envelopewavefunction overlap in Eq. (22) is (cid:104) ϕ i | ϕ j (cid:105) ≈ . Infact, for these states close to the band edge the index i is typically associated with a single quantum number, due to the typically small aspect ratio of QDs. Underabove assumptions, the Bloch equations simplify to ∂ t ρ s cv ,i = − i ω s cv ,i ρ s cv ,i + i (cid:126) − d cv (cid:0) − ρ s v ,i − ρ s c ,i (cid:1) E − γ cv ,i ρ s cv ,i ,∂ t ρ sα,i = 2 (cid:126) − (cid:61) (cid:8) d vc ρ s cv ,i (cid:9) E + (cid:2) ∂ t ρ sα,i (cid:3) intra − A s cv ,i ρ s c ,i ρ s v ,i , (175)where α = c , v , and d s cv ,i has been approximated by an s and i independent value d cv . For a closed description of7the carrier dynamics, Eq. (175) is again supplemented byEq. (173). The radiative and nonradiative transitionstaken into account in the resulting model are illustratedin Fig. 23. The RWA can be applied in the usual man-ner, as described in Section III E. The MB model hasbeen demonstrated to yield good agreement with exper-imental results for QD lasers and amplifiers, and to beinstrumental in interpreting the experimental findings.For example, the ultrafast gain dynamics in a QD am-plifier as well as the spatiotemporal dynamics and emis-sion characteristics of a QD laser were experimentallyand theoretically studied.
Furthermore, based on MBsimulations of ultrashort laser pulse propagation in a QDamplifier, it could be confirmed that the experimentallyobserved reshaping was in part due to coherent light–matter interaction.
As discussed in Section V C 2, longitudinal spatial holeburning, i.e., the formation of an inversion grating dueto the standing wave pattern in a Fabry-Pérot resonator,is automatically included in full-wave MB simulations.Assuming that tunneling between adjacent QDs can beneglected, the degradation of the inversion grating is gov-erned by carrier diffusion in the wetting layers, which canbe modeled by adding to Eq. (173) a diffusion term of theform Eq. (120).
X. CONCLUSION AND OUTLOOK
The goal of this review has been to discuss in detailthe underlying theoretical framework of the MB model,its extension and adaption to certain application areasand types of nanostructures, as well as special analyticalsolutions and suitable numerical methods. Apart fromthe intuitive appeal of the model and its adaptability, therelative compactness of the Bloch equations make themhighly suitable as an efficient quantum model for the ma-terial polarization in computational electrodynamics. Asshown in Section VII, their representation as a system ofordinary differential equations in time, where the posi-tion coordinates only enter as parameters, allows an effi-cient coupling to numerical schemes for Maxwell’s or re-lated propagation equations, such as the finite-differencetime-domain method. This compact form of the Blochequations is enabled by a mostly phenomenological treat-ment of dissipation based on the Lindblad formalism andrestriction to classical optical fields as well as discrete en-ergy levels. Fully microscopic descriptions of light-matterinteraction in a semiconductor, such as the semiconduc-tor MB equations shortly discussed in Section IX A, illustrate the limitations of semi-phenomenological Blochequations, and can serve as a starting point to developimproved compact Bloch equations. As an example, thisstrategy has been used to model the carrier dynamics ina semiconductor structure with a quasi-continuum of en-ergy levels in the conduction and valence band by macro-scopic discrete-level Bloch equations, as discussed inSection IX A. The main requirement for computational models isgenerally to combine numerical efficiency with accuracy,predictability and versatility. In this context, detailedmicroscopic theories can quickly become very compu-tationally demanding, which renders them impracticalfor applications such as device design.
Thus, a majorgoal is to further improve the quantitative accuracy andadaptability of the macroscopic MB equations by extend-ing the model accordingly, however without substantiallyincreasing its numerical complexity. This implies that itsgeneral form as a system of a few ordinary differentialequations should not be compromised.Probably the main limitation of the Bloch equationsis the phenomenological implementation of dissipationbased on the Lindblad formalism. As shortly discussedin Section IX C, an empirical treatment of certain pro-cesses, such as Pauli blocking or carrier-carrier scatter-ing, requires a generalization to nonlinear models. Here,special care must be taken to preserve the properties ofthe density matrix guaranteeing its physical character,which has for example been achieved by suitably extend-ing the Lindblad formalism. As mentioned in Section II C, the Lindblad model isonly realistic from a microscopic point of view if thememory decay of the environment occurs on a fastertimescale than the coherent system dynamics and relax-ation processes.
Although the macroscopic MB equa-tions often work surprisingly well on the verge of, or evenoutside, this microscopic validity range, advanced quan-titative modeling requires going beyond the Markovianapproximation in such cases. An ad hoc extension of theLindblad approach is obtained by replacing D k (ˆ ρ ) ( t ) inEq. (3) with (cid:82) t K k ( t − t (cid:48) ) D k (ˆ ρ ) ( t (cid:48) ) d t (cid:48) , where K k is thememory kernel. In certain cases, it is sufficient to treatthe populations in the usual manner and include mem-ory effects only for dephasing, which requires substitutingthe dephasing terms [d t ρ αβ ] relax = − γ αβ ρ αβ in Eq. (12)with − γ αβ (cid:82) t K αβ ( t − t (cid:48) ) ρ αβ ( t (cid:48) ) d t (cid:48) . The characteris-tic memory time and functional dependence of the mem-ory kernel, such as Gaussian or exponential, depend onthe underlying scattering mechanism.
Since theevaluation of convolution integrals is numerically expen-sive, a representation based on supplemental differentialequations is preferential. For exponential memory kernels K αβ = τ − αβ exp ( − t/τ αβ ) , such a differential equation iseasily derived, with d t s = − τ − αβ ( γ αβ ρ αβ + s ) where wehave introduced s = [d t ρ αβ ] relax for compactness. How-ever, since such modifications obviously do not preservethe Lindblad form of the dissipation terms, a physicalbehavior of the density matrix is not guaranteed, and infact highly nonphysical behavior can emerge. It should be pointed out that master equation mod-els with memory effects do not necessarily require con-volution integrals, and that memory kernel masterequations can even usually be cast into a time-localform. Thus, a promising approach towards a moregeneralized treatment of dissipation is to start with the8Lindblad equation in the form Eq. (4), and to general-ize the matrix D ijmn given in Eq. (6) for an arbitraryset of Lindblad operators. As already mentioned in Sec-tion II, time dependent Lindblad operators ˆ L k ( t ) , corre-sponding to time-varying dissipation rates in Eqs. (8) and(11), are unproblematic. Any further generalizationof D ijmn comes at the price of potentially unphysical re-sults. One example is the occurrence of temporarily neg-ative rates in Eqs. (8) or (11), which indeed introducesmemory effects into the Lindblad equation. Inthe construction of such a model, care should be takento avoid unphysical behavior, for example by adding cer-tain constraints.
Furthermore, a widely used modelof the form Eq. (4) is the Redfield equation, which is de-rived from microscopic considerations, i.e., a perturba-tive treatment of a quantum system weakly coupled tothe environment.
In this case D ijmn corresponds tothe generally time dependent Redfield tensor, which isdirectly related to the system-environment coupling andenvironment Hamiltonians.
The main advantages ofthe Redfield model are its strong connection to micro-scopic physics, and to some extent the inclusion of short-term memory effects.
However, in its commonlyused form, the Redfield equation does not guarantee pos-itivity of the density matrix which can lead to negativestate occupations. In practice, the emergence of this un-physical behavior appears to be a minor problem, and can also be cured.
Finally, applying the Lindblad formalism to a suitablyextended state space, a non-Markovian evolution witharbitrarily long memory times and strong initial correla-tions can be described.
This is achieved by represent-ing the reduced system density matrix ρ s with dimension N as a sum of a certain number M of positive matri-ces ρ i , i.e., ρ s = (cid:80) i ρ i where the traces of the ρ i mustadd up to one. Then, a big block diagonal density ma-trix ρ with dimension M N is constructed from the ρ i ,and the evolution of ρ is modeled by a Lindblad equa-tion for the extended system, where the operators arerequired to preserve the block diagonal form of ρ . Thisleads to M coupled evolution equations for the ρ i , wherethe dynamics is now defined by M arbitrary Hermitianoperators ˆ H i and M sets of arbitrary dissipation opera-tors (cid:110) ˆ L ij , . . . , ˆ L ijK (cid:111) . In this way, although the dynamicsof ρ is Markovian, the model can describe a highly non-Markovian evolution of ρ s , while intrinsically preservingthe physical properties of ρ s . An interesting subcase iswhen the evolution equations of the ρ i are decoupled, i.e., ˆ L ijk = ˆ L ijk δ ij . In this case, the evolution of each ρ i is described by an equation of the Lindblad form Eq. (3),but still a non-Markovian dynamics of ρ s is obtained. Ofcourse, for M = 1 , the standard Markovian Lindbladdynamics is recovered. CONFLICTS OF INTEREST
The authors declare no conflicts of interest.
SUPPORTING INFORMATION
Supporting information is available from the Wiley On-line Library or from the author.
ACKNOWLEDGMENTS
We thank Gabriela Slavcheva for valuable commentsand suggestions. The authors acknowledge financial sup-port by the European Union’s Horizon 2020 research andinnovation programme under grant agreement No 820419– Qombs Project ”Quantum simulation and entanglementengineering in quantum cascade laser frequency combs”(FET Flagship on Quantum Technologies), and by theGerman Research Foundation (DFG) within the Heisen-berg program (JI 115/4-2). N. Kirstaedter, N. Ledentsov, M. Grundmann, D. Bimberg,V. Ustinov, S. Ruvimov, M. Maximov, P. S. Kop’ev, Z. I.Alferov, U. Richter, P. Werner, U. Gosele, J. Heydenreich,
Elec-tron. Lett. , , 1416. N. Ledentsov, V. Ustinov, V. Shchukin, P. Kop’ev, Z. I. Alferov,D. Bimberg,
Semiconductors , , 343. D. Huffaker, G. Park, Z. Zou, O. Shchekin, D. Deppe,
Appl.Phys. Lett. , , 2564. J. P. Reithmaier, G. Eisenstein, A. Forchel,
Proc. IEEE , , 1779. X. Duan, Y. Huang, R. Agarwal, C. M. Lieber,
Nature , , 241. B. Mayer, A. Regler, S. Sterzl, T. Stettner, G. Koblmüller,M. Kaniber, B. Lingnau, K. Lüdge, J. Finley,
Nat. Commun. , , 15521. J. Faist, F. Capasso, D. L. Sivco, C. Sirtori, A. L. Hutchinson,A. Y. Cho,
Science , , 553. D. Hofstetter, M. Beck, J. Faist,
Appl. Phys. Lett. , ,2683. L. Gendron, M. Carras, A. Huynh, V. Ortiz, C. Koeniguer,V. Berger,
Appl. Phys. Lett. , , 2824. B. Levine,
J. Appl. Phys. , , R1. J. Phillips, K. Kamath, P. Bhattacharya,
Appl. Phys. Lett. , , 2020. S.-W. Lee, K. Hirakawa, Y. Shimada,
Appl. Phys. Lett. , , 1428. H. Liu, M. Gao, J. McCaffrey, Z. Wasilewski, S. Fafard,
Appl.Phys. Lett. , , 79. I. I. Rabi,
Phys. Rev. , , 652. A. Zrenner, E. Beham, S. Stufler, F. Findeis, M. Bichler, G. Ab-streiter,
Nature , , 612. S. Cundiff, A. Knorr, J. Feldmann, S. Koch, E. Göbel, H. Nickel,
Phys. Rev. Lett. , , 1178. A. Schülzgen, R. Binder, M. Donovan, M. Lindberg, K. Wundke,H. Gibbs, G. Khitrova, N. Peyghambarian,
Phys. Rev. Lett. , , 2346. H. Choi, V.-M. Gkortsas, L. Diehl, D. Bour, S. Corzine, J. Zhu,G. Höfler, F. Capasso, F. X. Kärtner, T. B. Norris,
Nature Pho-ton. , , 706. T. Stievater, X. Li, D. G. Steel, D. Gammon, D. Katzer, D. Park,C. Piermarocchi, L. Sham,
Phys. Rev. Lett. , , 133603. H. Kamada, H. Gotoh, J. Temmyo, T. Takagahara, H. Ando,
Phys. Rev. Lett. , , 246401. M. Kolarczik, N. Owschimikow, J. Korn, B. Lingnau, Y. Kap-tan, D. Bimberg, E. Schöll, K. Lüdge, U. Woggon,
Nat. Com-mun. , , 2953. O. Karni, A. Capua, G. Eisenstein, V. Sichkovskyi, V. Ivanov,J. P. Reithmaier,
Opt. Express , , 26786. A. Capua, O. Karni, G. Eisenstein, J. P. Reithmaier,
Phys. Rev.B , , 045305. S. L. McCall, E. L. Hahn,
Phys. Rev. Lett. , , 908. S. L. McCall, E. L. Hahn,
Phys. Rev. , , 457. S. Schneider, P. Borri, W. Langbein, U. Woggon, J. Förstner,A. Knorr, R. Sellin, D. Ouyang, D. Bimberg,
Appl. Phys. Lett. , , 3668. V. V. Kozlov,
Phys. Rev. A , , 1607. V. Kalosha, M. Müller, J. Herrmann,
J. Opt. Soc. Am. B , , 323. V. V. Kozlov, N. N. Rosanov, S. Wabnitz,
Phys. Rev. A , , 053810. C. R. Menyuk, M. A. Talukder,
Phys. Rev. Lett. , ,023903. R. Arkhipov, M. Arkhipov, I. Babushkin,
Opt. Commun. , , 73. L. V. Hau, S. E. Harris, Z. Dutton, C. H. Behroozi,
Nature , , 594. C. Liu, Z. Dutton, C. H. Behroozi, L. V. Hau,
Nature , , 490. D. Phillips, A. Fleischhauer, A. Mair, R. Walsworth, M. D.Lukin,
Phys. Rev. Lett. , , 783. J. B. Khurgin,
J. Opt. Soc. Am. B , , 1062. R. M. Camacho, C. J. Broadbent, I. Ali-Khan, J. C. Howell,
Phys. Rev. Lett. , , 043902. O. Firstenberg, M. Shuker, N. Davidson, A. Ron,
Phys. Rev.Lett. , , 043601. M. Lukin, S. Yelin, M. Fleischhauer,
Phys. Rev. Lett. , ,4232. A. Turukhin, V. Sudarshanam, M. Shahriar, J. Musser, B. Ham,P. Hemmer,
Phys. Rev. Lett. , , 023602. M. S. Bigelow, N. N. Lepeshkin, R. W. Boyd,
Phys. Rev. Lett. , , 113903. P. Ginzburg, M. Orenstein,
Opt. Express , , 12467. H. Borges, L. Sanz, J. Villas-Bôas, O. D. Neto, A. Alcalde,
Phys.Rev. B , , 115425. P. Tzenov, C. Jirauschek in
Proc. SPIE 10226, 19th Interna-tional Conference and School on Quantum Electronics: LaserPhysics and Applications , International Society for Optics andPhotonics, p. 1022603. J.-H. Wu, J.-Y. Gao, J.-H. Xu, L. Silvestri, M. Artoni,G. La Rocca, F. Bassani,
Phys. Rev. Lett. , , 057401. J. Gea-Banacloche, M. Mumba, M. Xiao,
Phys. Rev. B , , 165330. F. Bloch,
Phys. Rev. , , 460. R. P. Feynman, F. L. Vernon Jr, R. W. Hellwarth,
J. Appl.Phys. , , 49. I. Abella, N. Kurnit, S. Hartmann,
Phys. Rev. , , 391. L. Allen, J. H. Eberly,
Optical Resonance and Two-Level Atoms , Vol. 28 , Courier Corporation, . G. Lindblad,
Commun. Math. Phys. , , 119. V. Gorini, A. Kossakowski, E. C. G. Sudarshan,
J. Math. Phys. , , 821. A. Taflove, S. C. Hagness,
Computational Electrodynamics: TheFinite-Difference Time-Domain Method , Artech House, . G. Slavcheva, J. M. Arnold, I. Wallace, R. W. Ziolkowski,
Phys.Rev. A , , 063418. A. Klaedtke, O. Hess,
Opt. Express , , 2744. M. Sukharev, A. Nitzan,
Phys. Rev. A , , 043802. A. Pusch, S. Wuestner, J. M. Hamm, K. L. Tsakmakidis,O. Hess,
ACS Nano , , 2420. K. Lopata, D. Neuhauser,
J. Chem. Phys. , , 014701. H. Takeda, S. John,
Phys. Rev. A , , 053811. M. Dridi, G. C. Schatz,
J. Opt. Soc. Am. B , , 2791. W. Cartar, J. Mørk, S. Hughes,
Phys. Rev. A , , 023859. M. Riesch, P. Tzenov, C. Jirauschek in , IEEE, pp. 1–4. P. Tzenov, D. Burghoff, Q. Hu, C. Jirauschek,
Opt. Express , , 23232. C. M. Bowden, J. P. Dowling,
Phys. Rev. A , , 1247. G. Y. Slepyan, S. Maksimenko, A. Hoffmann, D. Bimberg,
Phys.Rev. A , , 063804. W. W. Chow, S. W. Koch, M. I. Sargent,
Semiconductor-LaserPhysics , Springer Science & Business Media, . H. Haug, S. W. Koch,
Quantum Theory of the Optical and Elec-tronic Properties of Semiconductors , World Scientific Publish-ing Company, . C. Ning, R. Indik, J. Moloney,
IEEE J. Quantum Electron. , , 1543. J. Yao, G. P. Agrawal, P. Gallion, C. M. Bowden,
Opt. Com-mun. , , 246. S. Balle,
Opt. Commun. , , 227. R. Rosati, F. Rossi,
Phys. Rev. B , , 205415. D. Stepanenko, G. Burkard, G. Giedke, A. Imamoglu,
Phy. Rev.Lett. , , 136401. S. G. Schirmer, A. I. Solomon,
Phys. Rev. A , , 022107. R. Rosati, R. C. Iotti, F. Dolcini, F. Rossi,
Phys. Rev. B , , 125140. A. Steinhoff, P. Gartner, M. Florian, F. Jahnke,
Phys. Rev. B , , 205144. H.-P. Breuer, F. Petruccione,
The Theory of Open QuantumSystems , Oxford University Press, . J. Preskill,
Lecture Notes for Physics 229: Quantum Informa-tion and Computation , California Institute of Technology, . K. Kraus,
Ann. Phys. (N. Y.) , , 311. H.-P. Breuer,
Phys. Rev. A , , 012106. C. M. Kropf, C. Gneiting, A. Buchleitner,
Phys. Rev. X , , 031023. E. B. Davies,
Commun. Math. Phys. , , 91. M. Le Bellac,
Quantum Physics , Cambridge University Press, . M. Kantner, U. Bandelow, T. Koprucki, H.-J. Wünsche in
Nu-merical Simulation of Optoelectronic Devices (NUSOD), 2015International Conference on , IEEE, pp. 151–152. D. Chruściński, A. Kossakowski,
Phys. Rev. Lett. , ,070406. M. Grifoni, P. Hänggi,
Phys. Rep. , , 229. D. K. Oi, S. G. Schirmer,
Phys. Rev. A , , 012121. P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, A. Aspuru-Guzik,
New J. Phys. , , 033003. S. Fathololoumi, E. Dupont, C. Chan, Z. Wasilewski, S. Lafram-boise, D. Ban, A. Mátyás, C. Jirauschek, Q. Hu, H. Liu,
Opt.Express , , 3866. T. Dinh, A. Valavanis, L. Lever, Z. Ikonić, R. Kelsall,
Phys.Rev. B , , 235427. E. Dupont, S. Fathololoumi, H. C. Liu,
Phys. Rev. B , ,205311. B. A. Burnett, B. S. Williams,
Phys. Rev. B , , 155309. H. Callebaut, Q. Hu,
J. Appl. Phys. , , 104505. G. Pfanner, M. Seliger, U. Hohenester,
Phys. Rev. B , ,195410. B. Palmieri, D. Abramavicius, S. Mukamel,
J. Chem. Phys. , , 204512. G. Kiršanskas, M. Franckié, A. Wacker,
Phys. Rev. B , ,035432. S. Kumar, Q. Hu,
Phys. Rev. B , , 245316. R. Terazzi, J. Faist,
New J. Phys. , , 033045. P. Tzenov, D. Burghoff, Q. Hu, C. Jirauschek,
IEEE Trans.Terahertz Sci. Technol. , , 351. C. Jirauschek,
J. Appl. Phys. , , 133105. A. Gordon, D. Majer,
Phys. Rev. B , , 195317. Y. Dubi, M. Di Ventra,
Nano Lett. , , 97. W. Freeman,
Phys. Rev. B , , 205301. R. W. Ziolkowski, J. M. Arnold, D. M. Gogny,
Phys. Rev. A , , 3082. S. Hughes,
Phys. Rev. Lett. , , 3363. V. Kalosha, J. Herrmann,
Phys. Rev. Lett. , , 544. O. Mücke, T. Tritschler, M. Wegener, U. Morgner, F. Kärtner,
Phys. Rev. Lett. , , 057401. J. R. Freeman, J. Maysonnave, S. Khanna, E. H. Linfield, A. G.Davies, S. S. Dhillon, J. Tignon,
Phys. Rev. a , , 063817. S. Mukamel,
Annu. Rev. Phys. Chem. , , 647. N. Gisin, I. C. Percival,
Journal of Physics A: Mathematicaland General , , 5677. P. Meystre, M. Sargent,
Elements of Quantum Optics , SpringerScience & Business Media, . M. O. Scully, M. S. Zubairy,
Quantum Optics , Cambridge Uni-versity Press, . M. W. Walser, C. H. Keitel, A. Scrinzi, T. Brabec,
Phys. Rev.Lett. , , 5082. S. Stobbe, P. T. Kristensen, J. E. Mortensen, J. M. Hvam,J. Mørk, P. Lodahl,
Phys. Rev. B , , 085304. P. Lodahl, S. Mahmoodian, S. Stobbe,
Rev. Mod. Phys. , , 347. K. Rzązewski, R. W. Boyd,
J. Mod. Opt. , , 1137. D. Bauer, D. Milošević, W. Becker,
Phys. Rev. A , ,023415. T. Yoshie, A. Scherer, J. Hendrickson, G. Khitrova, H. Gibbs,G. Rupper, C. Ell, O. Shchekin, D. Deppe,
Nature , ,200. D. Englund, A. Majumdar, A. Faraon, M. Toishi, N. Stoltz,P. Petroff, J. Vučković,
Phys. Rev. Lett. , , 073904. O. Hess, T. Kuhn,
Phys. Rev. A , , 3347. C. Jirauschek, T. Kubis,
Appl. Phys. Rev. , , 011307. M. A. Talukder, C. R. Menyuk,
Opt. Express , , 15608. R. Arkhipov, M. Arkhipov, I. Babushkin, N. Rosanov,
Opt. Lett. , , 737. P. Tzenov, I. Babushkin, R. Arkhipov, M. Arkhipov,N. Rosanov, U. Morgner, C. Jirauschek,
New J. Phys. , , 053055. A.-B. Chen, A. Sher,
Phys. Rev. B , , 3886. G. Bastard,
Wave Mechanics Applied to Semiconductor Het-erostructures , Les Editons de Physique, Paris, . S. Datta,
Quantum Transport: Atom to Transistor , CambridgeUniversity Press, . C. Jirauschek,
IEEE J. Quantum Electron. , , 1059. U. Ekenberg,
Phys. Rev. B , , 7714. D. F. Nelson, R. C. Miller, D. A. Kleinman,
Phys. Rev. B , , 7770. E. O. Kane,
J. Phys. Chem. Solids , , 249. E. O. Kane,
Handbook on Semiconductors , , 193. J. M. Luttinger, W. Kohn,
Phys. Rev. , , 869. G. L. Bir, G. E. Pikus,
Symmetry and Strain-Induced Effects inSemiconductors , Wiley New York, . C. R. Pidgeon, R. Brown,
Phys. Rev. , , 575. O. Stier, M. Grundmann, D. Bimberg,
Phys. Rev. B , ,5688. G. Baraff, D. Gershoni,
Phys. Rev. B , , 4011. D. J. Paul,
Phys. Rev. B , , 155323. B. A. Foreman,
Phys. Rev. B , , 4964. U. Rössler,
Solid State Commun. , , 943. S. Ridene, K. Boujdaria, H. Bouchriha, G. Fishman,
Phys. Rev.B , , 085329. J. H. Davies,
The Physics of Low-dimensional Semiconductors ,Cambridge University Press, Cambridge, . M. Burt,
J. Phys. Condens. Matter , , 4091. B. D. Gerardot, D. Brunner, P. A. Dalgarno, P. Öhberg, S. Seidl,M. Kroner, K. Karrai, N. G. Stoltz, P. M. Petroff, R. J. War-burton,
Nature , , 441. K. Karlsson, V. Troncale, D. Oberli, A. Malko, E. Pelucchi,A. Rudra, E. Kapon,
Appl. Phys. Lett. , , 251113. Y.-M. Niquet, D. C. Mojica,
Phys. Rev. B , , 115316. S. Cortez, O. Krebs, P. Voisin, J. Gérard,
Phys. Rev. B , , 233306. D. Bimberg, M. Grundmann, N. N. Ledentsov,
Quantum DotHeterostructures , John Wiley & Sons, . D. Bimberg,
Semiconductors , , 951. U. Woggon,
Optical Properties of Semiconductor QuantumDots , Springer, . J. Finley in
Handbook of Self Assembled Semiconductor Nanos-tructures for Novel Devices in Photonics and Electronics , Else-vier, , pp. 476–504.
F. T. Hioe, J. H. Eberly,
Phys. Rev. Lett. , , 838. G. Slavcheva, J. M. Arnold, R. W. Ziolkowski,
IEEE J. Sel.Top. Quantum Electron. , , 929. G. Slavcheva,
Phys. Rev. B , , 115347. G. Slavcheva, P. Roussignol,
New J. Phys. , , 103004. R. Marskar, U. Österberg,
Opt. Express , , 16784. F. Bloch, A. Siegert,
Phys. Rev. , , 522. N. Bloembergen,
Nonlinear Optics , World Scientific, Singapore, . R. W. Boyd,
Nonlinear Optics , Academic, . A. Wacker,
Phys. Rep. , , 1. M. Wegener,
Extreme Nonlinear Optics: An Introduction ,Springer Science & Business Media, . N. Owschimikow, C. Gmachl, A. Belyanin, V. Kocharovsky,D. L. Sivco, R. Colombelli, F. Capasso, A. Y. Cho,
Phys. Rev.Lett. , , 043902. D. Englund, A. Faraon, I. Fushman, N. Stoltz, P. Petroff,J. Vučković,
Nature , , 857. K. Srinivasan, O. Painter,
Nature , , 862. M. Belkin, F. Capasso, F. Xie, A. Belyanin, M. Fischer,A. Wittmann, J. Faist,
Appl. Phys. Lett. , , 201101. W. Maineult, L. Ding, P. Gellie, P. Filloux, C. Sirtori, S. Bar-bieri, T. Akalin, J.-F. Lampin, I. Sagnes, H. E. Beere, D. A.Ritchie,
Appl. Phys. Lett. , , 021108. A. Calvar, M. Amanti, M. Renaudat St-Jean, S. Barbieri,A. Bismuto, E. Gini, M. Beck, J. Faist, C. Sirtori,
Appl. Phys.Lett. , , 181114. M. R. St-Jean, M. I. Amanti, A. Bernard, A. Calvar, A. Bis-muto, E. Gini, M. Beck, J. Faist, H. Liu, C. Sirtori,
Laser Pho-ton. Rev. , , 443. F. Wang, K. Maussang, S. Moumdji, R. Colombelli, J. R. Free-man, I. Kundu, L. Li, E. H. Linfield, A. G. Davies, J. Mangeney,J. Tignon, S. S. Dhillon,
Optica , , 944. J. Faist, G. Villares, G. Scalari, M. Rösch, C. Bonzon, A. Hugi,M. Beck,
Nanophotonics , , 272. G. Grau, W. Freude,
Optische Nachrichtentechnik - Eine Ein-führung , Springer, Berlin, . A. Yariv,
Quantum Electronics , John Wiley & Sons, New York, . M. Sargent III, M. O. Scully, J. W. E. Lamb,
Laser Physics ,Addison-Wesley, Reading, MA, . R. G. Brewer, R. Shoemaker,
Phys. Rev. A , , 2001. G. M. Slavcheva, J. M. Arnold, R. W. Ziolkowski,
IEEE J. Sel.Top. Quantum Electron. , , 1052. P. Siddons,
J. Phys. B , , 093001. G. Slavcheva, O. Hess,
Phys. Rev. A , , 053804. G. Slavcheva, M. Koleva, A. Rastelli,
Phys. Rev. B , ,115433. X. Song, S. Gong, R. Li, Z. Xu,
Phys. Rev. A , , 015802. N. Schulz, K. Bierwirth, F. Arndt, U. Koster,
IEEE Trans. Mi-crow. Theory Techn. , , 722. A. S. Sudbo,
Pure Appl. Opt. , , 381. K. Chiang,
Opt. Lett. , , 714. G. Agrawal,
Nonlinear Fiber Optics , Academic, New York, . E. Schrödinger,
Ann. Phys. (Berl.) , , 437. M. M. Sternheim, J. F. Walker,
Phys. Rev. C , , 114. T. Visser, H. Blok, B. Demeulenaere, D. Lenstra,
IEEE J.Quantum Electron. , , 1763. J. D. Jackson,
Classical Electrodynamics , Wiley & Sons, NewYork, . Q. Y. Lu, N. Bandyopadhyay, S. Slivken, Y. Bai, M. Razeghi,
Appl. Phys. Lett. , , 131106. M. Born, E. Wolf,
Principles of Optics: Electromagnetic The-ory of Propagation, Interference and Diffraction of Light , Cam-bridge University Press, Cambridge, . R. Bräuer, O. Bryngdahl,
Appl. Opt. , , 7875. A. Maimistov, E. Manykin,
Zh. Eksp. Teor. Fiz , , 1177. E. Doktorov, R. Vlasov,
Opt. Acta , , 223. M. Nakazawa, E. Yamada, H. Kubota,
Phys. Rev. Lett. , , 2625. R. Guo, H.-Q. Hao,
Ann. Phys. (N. Y.) , , 10. F. DeMartini, C. Townes, T. Gustafson, P. Kelley,
Phys. Rev. , , 312. F. M. Mitschke, L. F. Mollenauer,
Opt. Lett. , , 659. J. P. Gordon,
Opt. Lett. , , 662. K. Nakkeeran, K. Porsezian,
J. Phys. A , , 3817. S. Kohen, B. S. Williams, Q. Hu,
J. Appl. Phys. , ,053106. J. Butler, J. Zoroofchi,
IEEE J. Quantum Electron. , ,809. T. Ikegami,
IEEE J. Quantum Electron. , , 470. P. C. Kendall, D. A. Roberts, P. N. Robson, M. J. Adams, M. J.Robertson,
IEEE Photon. Technol. Lett. , , 148. F. K. Reinhart, I. Hayashi, M. B. Panish,
J. Appl. Phys. , , 4466. C. Y. Wang, L. Diehl, A. Gordon, C. Jirauschek, F. X. Kärt-ner, A. Belyanin, D. Bour, S. Corzine, G. Höfler, M. Troccoli,J. Faist, F. Capasso,
Phys. Rev. A , , 031802. O. Hess, T. Kuhn,
Phys. Rev. A , , 3360. H. C. Torrey,
Phys. Rev. , , 563. T. Ando, A. B. Fowler, F. Stern,
Rev. Mod. Phys. , ,437. L. V. Asryan, R. A. Suris,
IEEE J. Quantum Electron. , , 1151. A. Capua, O. Karni, G. Eisenstein,
IEEE J. Sel. Top. QuantumElectron. , , 1900410. A. Gordon, C. Y. Wang, L. Diehl, F. Kärtner, A. Belyanin,D. Bour, S. Corzine, G. Höfler, H. Liu, H. Schneider, T. Maier,M. Troccoli, J. Faist, F. Capasso,
Phys. Rev. A , ,053804. N. Vukovic, J. Radovanovic, V. Milanovic, D. Boiko,
Opt.Quant. Electron. , , 254. H. Torrey,
Phys. Rev. , , 1059. A. Muller, E. B. Flagg, P. Bianucci, X. Wang, D. G. Deppe,W. Ma, J. Zhang, G. Salamo, M. Xiao, C.-K. Shih,
Phys. Rev.Lett. , , 187402. X. Xu, B. Sun, P. R. Berman, D. G. Steel, A. S. Bracker,D. Gammon, L. J. Sham,
Science , , 929. M. Wagner, H. Schneider, D. Stehr, S. Winnerl, A. M. Andrews,S. Schartner, G. Strasser, M. Helm,
Phys. Rev. Lett. , ,167401. C. Henry,
IEEE J. Quantum Electron. , , 259. D. Bimberg, N. Kirstaedter, N. Ledentsov, Z. I. Alferov,P. Kop’ev, V. Ustinov,
IEEE J. Sel. Top. Quantum Electron. , , 196. C. Jirauschek,
Opt. Express , , 25922. K.-J. Boller, A. Imamoğlu, S. E. Harris,
Phys. Rev. Lett. , , 2593. M. Fleischhauer, A. Imamoglu, J. P. Marangos,
Rev. Mod. Phys. , , 633. A. Kasapi, M. Jain, G. Yin, S. E. Harris,
Phys. Rev. Lett. , , 2447. G. P. Agrawal,
J. Opt. Soc. Am. B , , 147. Y. Hu, M. Lindberg, S. Koch,
Phys. Rev. B , , 1713. M. Sugawara, H. Ebe, N. Hatori, M. Ishida, Y. Arakawa,T. Akiyama, K. Otsubo, Y. Nakata,
Phys. Rev. B , ,235332. J. Khurgin,
J. Opt. Soc. Am. B , , 1673. M. Belkin, F. Capasso, A. Belyanin, D. Sivco, A. Cho, D. Oak-ley, C. Vineis, G. Turner,
Nat. Photon. , , 288. B. A. Burnett, B. S. Williams,
Phys. Rev. Appl. , ,034013. C. Jirauschek, H. Okeil, P. Lugli,
Opt. Express , , 1670. A. Vizbaras, M. Anders, S. Katz, C. Grasse, G. Boehm,R. Meyer, M. A. Belkin, M.-C. Amann,
IEEE J. Quantum Elec-tron. , , 691. P. N. Butcher, D. Cotter,
The Elements of Nonlinear Optics , Vol. 9 , Cambridge University Press, . Y. Shen,
The Principles of Nonlinear Optics , Wiley-Interscience, New York, . G. Lamb,
Rev. Mod. Phys. , , 99. A. Maimistov, A. Basharov, S. Elyutin, Y. M. Sklyarov,
Phys.Rep. , , 1. M. A. Talukder, C. R. Menyuk,
Appl. Phys. Lett. , ,071109. M. A. Talukder, C. R. Menyuk,
Opt. Express , , 5639. B. R. Mollow,
Phys. Rev. , , 1969. R. Bullough, P. Jack, P. Kitchenside, R. Saunders,
Phys. Scr. , , 364. A. Kujawski, J. Mostowski,
J. Opt. Soc. Am. B , , 1700. J. Javaloyes, S. Balle,
Freetwm: a simulation tool for semi-conductor lasers , https://onl.uib.eu/Softwares/Freetwm/ , . Kintechlab,
Electromagnetic Template Library , http://fdtd.kintechlab.com/en/start , . A. Deinega, T. Seideman,
Phys. Rev. A , , 022501. A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel,J. Joannopoulos, S. G. Johnson,
Computer Physics Commu-nications , , 687. M. Riesch, C. Jirauschek, mbsolve: An open-source solvertool for the Maxwell-Bloch equations , https://github.com/mriesch-tum/mbsolve , . Quantopticon,
Quantillion Software , https://quantopticon.co.uk , . P. D. Lax, B. Wendroff,
Commun. Pure Appl. Math. , ,217. H. Risken, K. Nummedal,
J. Appl. Phys. , , 4662. E. Isaacson, H. B. Keller,
Analysis of Numerical Methods , JohnWiley & Sons, . A. Harten,
J. Comput. Phys. , , 357. S. K. Godunov,
Matematicheskii Sbornik , , 271. P. K. Nielsen, H. Thyrrestrup, J. Mørk, B. Tromborg,
Opt. Ex-press , , 6396. J. Xiong, M. Colice, F. Schlottau, K. Wagner, B. Fornberg,
Opt.Quant. Electron. , , 447. G. Demeter,
Comput. Phys. Commun. , , 1203. J. C. Butcher,
Numerical Methods for Ordinary DifferentailEquations , John Wiley, . A. Iserles,
A First Course in the Numerical Analysis of Differ-ential Equations , Cambridge University Press, . C. W. Gear,
Numerical Initial Value Problems in Ordinary Dif-ferential Equations , Prentice Hall, . P. Arve, P. Jänes, L. Thylén,
Phys. Rev. A , , 063809. B. Gross, J. T. Manassah,
Opt. Lett. , , 340. B. Bidégaray, A. Bourgeade, D. Reignier,
J. Comput. Phys. , , 603. B. Bidégaray,
Numer. Methods Partial Differ. Equ. , ,284. K. Yee,
IEEE Trans. Antennas. Propag. , , 302. S. Krishnamoorthy, M. Baskaran, U. Bondhugula, J. Ramanu-jam, A. Rountev, P. Sadayappan,
SIGPLAN Not. , ,235. M. Riesch, N. Tchipev, S. Senninger, H.-J. Bungartz, C. Ji-rauschek,
Opt. Quant. Electron. , , 112. F. Schlottau, M. Piket-May, K. Wagner,
Opt. Express , ,182. Q. H. Liu,
Microw. Opt. Technol. Lett. , , 158. O. Saut, A. Bourgeade,
J. Comput. Phys. , , 823. B. Garraway, P. Knight,
Phys. Rev. A , , 1266. E. Hairer, S. P. Nørsett, G. Wanner,
Solving Ordinary Differ-ential Equations I 2nd ed. , Springer-Verlag Berlin Heidelberg, . G. Strang,
SIAM J. Numer. Anal. , , 506. J. Gallier, D. Xu,
Int. J. Robot. Autom. , , 10. C. Moler, C. V. Loan,
SIAM Rev. , , 3. D. H. Hailu, https: // arxiv. org/ abs/ 1610. 05951 . C. Weninger, N. Rohringer,
Phys. Rev. A , , 053421. L. Guduff, A. J. Allami, C. van Heijenoort, J.-N. Dumez,I. Kuprov,
Phys. Chem. Chem. Phys. , , 17577. R. Kosloff,
Annu. Rev. Phys. Chem. , , 145. M. Riesch, C. Jirauschek,
J. Comput. Phys. , , 290. S.-L. Chua, Y. Chong, A. D. Stone, M. Soljačić, J. Bravo-Abad,
Opt. Express , , 1539. C. Leforestier, R. H. Bisseling, C. Cerjan, M. D. Feit, R. Fries-ner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero, R. Kosloff,
J. Comput. Phys. , , 59. U. Peskin, R. Kosloff, N. Moiseyev,
J. Chem. Phys. , ,8849. J. C. Tremblay, T. Carrington Jr,
J. Chem. Phys. , ,11535. A. Gordon, C. Jirauschek, F. X. Kärtner,
Phys. Rev. A , , 042505. S. Blanes, F. Casas, A. Murua,
J. Chem. Phys. , ,114109. L. Pierantoni, D. Mencarelli, T. Rozzi,
IEEE Trans. Microw.Theory Tech. , , 654. I. Ahmed, E. H. Khoo, E. Li, R. Mittra,
IEEE Antennas Wirel.Propag. Lett. , , 914. I. P. Christov, M. M. Murnane, H. C. Kapteyn,
Phys. Rev. A , , R2285. E. Lorin, S. Chelkowski, A. Bandrauk,
Comput. Phys. Com-mun. , , 908. Q. Chen, H. Qin, J. Liu, J. Xiao, R. Zhang, Y. He, Y. Wang,
J.Comput. Phys. , , 441. Y. P. Chen, W. E. Sha, L. Jiang, M. Meng, Y. M. Wu, W. C.Chew,
Comput. Phys. Commun. , , 63. D. Masiello, E. Deumens, Y. Öhrn,
Phys. Rev. A , ,032108. E. Schelew, R.-C. Ge, S. Hughes, J. Pond, J. F. Young,
Phys.Rev. A , , 063853. M. Wegener, D. Chemla, S. Schmitt-Rink, W. Schäfer,
Phys.Rev. A , , 5675. G. Y. Slepyan, A. Magyarov, S. Maksimenko, A. Hoffmann,D. Bimberg,
Phys. Rev. B , , 045320. K. Dolgaleva, R. W. Boyd,
Adv. Opt. Photonics , , 1. C. M. Bowden, J. P. Dowling,
Phys. Rev. A , , 1514. K. Xia, S. Gong, C. Liu, X. Song, Y. Niu,
Opt. Express , , 5913. M. E. Crenshaw, C. M. Bowden,
Phys. Rev. A , , 1139. M. E. Crenshaw, K. U. Sullivan, C. M. Bowden,
Opt. Express , , 152. J. P. Dowling, C. M. Bowden,
Phys. Rev. Lett. , , 1421. A. Afanas’ev, R. Vlasov, A. Cherstvyi,
J. Exp. Theor. Phys. , , 428. M. Crenshaw, M. Scalora, C. M. Bowden,
Phys. Rev. Lett. , , 911. A. Afanas’ev, R. Vlasov, O. K. Khasanov, T. Smirnova, O. Fe-dotova,
J. Opt. Soc. Am. B , , 911. D. V. Novitsky,
Phys. Rev. A , , 013817. E. Paspalakis, A. Kalini, A. F. Terzis,
Phys. Rev. B , ,073305. Y. Mitsumori, S. Watanabe, K. Asakura, K. Seki, K. Edamatsu,K. Akahane, N. Yamamoto,
Phys. Rev. B , , 235305. A. E. Siegman,
Lasers , University Science Books, . N. Majer, K. Lüdge, E. Schöll,
Phys. Rev. B , , 235301. D. Polder, M. Schuurmans, Q. Vrehen,
Phys. Rev. A , ,1192. K. Wodkiewicz,
Phys. Rev. A , , 1686. J. Andreasen, H. Cao,
J. Light. Technol. , , 4530. C. Gardiner, P. Zoller,
Quantum Noise: A Handbook of Marko-vian and Non-Markovian Quantum Stochastic Methods withApplications to Quantum Optics , Vol. 56 , Springer Science &Business Media, . M. Lax, W. Louisell,
Phys. Rev. , , 568. P. Drummond, M. Raymer,
Phys. Rev. A , , 2072. E. Gehrig, O. Hess,
Phys. Rev. A , , 033804. J. E. Kim, E. Malic, M. Richter, A. Wilms, A. Knorr,
IEEE J.Quantum Electron. , , 1115. S. Wilkinson, B. Lingnau, J. Korn, E. Schöll, K. Lüdge,
IEEEJ. Sel. Top. Quantum Electron. , , 1900106. H. F. Hofmann, O. Hess,
Phys. Rev. A , , 2342. Q. Vu, H. Haug, O. Mücke, T. Tritschler, M. Wegener,G. Khitrova, H. Gibbs,
Phys. Rev. Lett. , , 217403. B. Witzigmann, V. Laino, M. Luisier, U. Schwarz, G. Feicht,W. Wegscheider, K. Engl, M. Furitsch, A. Leber, A. Lell,V. Härle,
Appl. Phys. Lett. , , 021104. F. Rossi, E. Molinari,
Phys. Rev. Lett. , , 3642. D. H. Marti, M. Dupertuis, B. Deveaud,
IEEE J. QuantumElectron. , , 848. D. Golde, T. Meier, S. W. Koch,
Phys. Rev. B , , 075330. T. Stroucken, J. Grönqvist, S. Koch,
J. Opt. Soc. Am. B , , A86. M. Hirtschulz, F. Milde, E. Malić, S. Butscher, C. Thomsen,S. Reich, A. Knorr,
Phys. Rev. B , , 035403. I. Waldmueller, W. W. Chow, E. W. Young, M. C. Wanke,
IEEEJ. Quantum Electron. , , 292. C. Jirauschek, P. Tzenov,
Opt. Quant. Electron. , , 414. M. S. Vitiello, G. Scamarcio, V. Spagnolo, B. S. Williams, S. Ku-mar, Q. Hu, J. L. Reno,
Appl.Phys. Lett. , , 111115. R. Nelander, A. Wacker,
J. Appl. Phys. , , 063115. D. Dietze, A. Benz, G. Strasser, K. Unterrainer, J. Darmo,
Opt.Express , , 13700. Y. Wang, A. Belyanin,
Opt. Express , , 4173. T. Ando,
J. Phys. Soc. Jpn. , , 765. T. Unuma, M. Yoshita, T. Noda, H. Sakaki, H. Akiyama,
J.Appl. Phys. , , 1586. V.-M. Gkortsas, C. Wang, L. Kuznetsova, L. Diehl, A. Gordon,C. Jirauschek, M. A. Belkin, A. Belyanin, F. Capasso, F. X.Kärtner,
Opt. Express , , 13616. A. K. Wójcik, P. Malara, R. Blanchard, T. S. Mansuripur,F. Capasso, A. Belyanin,
Appl. Phys. Lett. , , 231102. D. Revin, M. Hemingway, Y. Wang, J. Cockburn, A. Belyanin,
Nat. Commun. , , 11440. L. Columbo, S. Barbieri, C. Sirtori, M. Brambilla,
Opt. Express , , 2829. N. N. Vuković, J. Radovanović, V. Milanović, D. L. Boiko,
IEEEJ. Sel. Top. Quantum Electron. , , 1200616. J. Khurgin, Y. Dikmelik, A. Hugi, J. Faist,
Appl. Phys. Lett. , , 081118. G. Villares, J. Faist,
Opt. Express , , 1651. J. Bai, H. Wang, Q. Wang, D. Zhou, K. Q. Le, B. Wang,
IEEEJ. Quantum Electron. , , 1. R. C. Iotti, F. Rossi,
EPL , , 67005. S. Butscher, J. Förstner, I. Waldmüller, A. Knorr,
Phys. Rev.B , , 045314. I. Savić, N. Vukmirović, Z. Ikonić, D. Indjin, R. W. Kelsall,P. Harrison, V. Milanović,
Phys. Rev. B , , 165310. D. Burghoff, T.-Y. Kao, N. Han, C. W. I. Chan, X. Cai, Y. Yang,D. J. Hayton, J.-R. Gao, J. L. Reno, Q. Hu,
Nature Photon. , , 462. N. S. Wingreen, C. A. Stafford,
IEEE J. Quantum Electron. , , 1170. E. Zibik, T. Grange, B. Carpenter, N. Porter, R. Ferreira,G. Bastard, D. Stehr, S. Winnerl, M. Helm, H. Liu, M. S. Skol-nick, L. R. Wilson,
Nat. Mater. , , 803. N. Zhuo, F. Q. Liu, J. C. Zhang, L. J. Wang, J. Q. Liu, S. Q.Zhai, Z. G. Wang,
Nanoscale Res. Lett. , , 144. E. Gehrig, O. Hess, C. Ribbat, R. Sellin, D. Bimberg,
Appl.Phys. Lett. , , 1650. C. Sailliot, V. Voignier, G. Huyet,
Opt. Commun. , ,353. J. Mukherjee, J. G. McInerney,
Phys. Rev. A , , 053813. I. E. Protsenko, A. V. Uskov, O. Zaimidoroga, V. Samoilov,E. P. O’Reilly,
Phys. Rev. A , , 063812. M. Kulkarni, O. Cotlet, H. E. Türeci,
Phys. Rev. B , ,125402. E. Waks, D. Sridharan,
Phys. Rev. A , , 043845. G. T. Adamashvili, C. Weber, A. Knorr, N. T. Adamashvili,
Phys. Rev. A , , 063808. L. Schneebeli, T. Feldtmann, M. Kira, S. W. Koch, N. Peygham-barian,
Phys. Rev. A , , 053852. W. W. Chow, H. Schneider, M. Phillips,
Phys. Rev. A , , 053802. B. Bidégaray-Fesquet, K. Keita,
J. Math. Phys. , ,021501. P. Bardella, L. L. Columbo, M. Gioannini,
Opt. Express , , 26234. T. R. Nielsen, P. Gartner, F. Jahnke,
Phys. Rev. B , ,235314. H. H. Nilsson, J.-Z. Zhang, I. Galbraith,
Phys. Rev. B , , 205331. T. Koprucki, A. Wilms, A. Knorr, U. Bandelow,
Opt. Quant.Electron. , , 777. B. Lingnau, K. Lüdge, W. W. Chow, E. Schöll,
Phys. Rev. E , , 065201. D. Hadass, A. Bilenca, R. Alizon, H. Dery, V. Mikhelashvili,G. Eisenstein, R. Schwertberger, A. Somers, J. P. Reithmaier,A. Forchel, M. Calligaro, S. Bansropun, M. Krakowski,
IEEE J.Sel. Top. Quant. , , 1015. N. Yasuoka, K. Kawaguchi, H. Ebe, T. Akiyama, M. Ekawa,S. Tanaka, K. Morito, A. Uetake, M. Sugawara, Y. Arakawa,
Appl. Phys. Lett. , , 101108. M. van der Poel, E. Gehrig, O. Hess, D. Birkedal, J. M. Hvam,
IEEE J. Quantum Electron. , , 1115. A. Valavanis, L. Lever, C. Evans, Z. Ikonić, R. Kelsall,
Phys.Rev. B , , 035420. S. M. Barnett, S. Stenholm,
Phys. Rev. A , , 033808. G. Adamashvili, D. Kaup, A. Knorr, C. Weber,
Phys. Rev. A , , 013840. X. Hu, W. Pötz in
Coherent Control in Atoms, Molecules, andSemiconductors , W. Pötz, W. Schroeder (Eds.), Springer, ,pp. 127–145.
E. Laine, K. Luoma, J. Piilo,
J. Phys. B , , 154004. E.-M. Laine, J. Piilo, H.-P. Breuer,
Phys. Rev. A , ,062115. X.-M. Lu, X. Wang, C. Sun,
Phys. Rev. A , , 042103. J.-S. Tang, C.-F. Li, Y.-L. Li, X.-B. Zou, G.-C. Guo, H.-P.Breuer, E.-M. Laine, J. Piilo,
EPL , , 10002. D. Chruściński, S. Maniscalco,
Phys. Rev. Lett. , ,120404. A. G. Redfield,
IBM J. Res. Dev. , , 19. R. K. Wangsness, F. Bloch,
Phys. Rev. , , 728. R. S. Whitney,
J. Phys. A , , 175304. D. Egorova, M. Thoss, W. Domcke, H. Wang,
J. Chem. Phys. , , 2761. I. Knezevic, B. Novakovic,
J. Comput. Electron. , , 363. C. Weber, A. Wacker, A. Knorr,
Phys. Rev. B , , 165322. E. Vaz, J. Kyriakidis,
J. Phys. Conf. Ser. , , 012012. A. Pan, B. A. Burnett, C. O. Chui, B. S. Williams,
Phys. Rev.B , , 085308. A. Suárez, R. Silbey, I. Oppenheim,
J. Chem. Phys. , ,5101. H.-P. Breuer,
Phys. Rev. A , , 022103. A. A. Budini,