Modeling of Fluctuations in Dynamical Optoelectronic Device Simulations within a Maxwell-Density Matrix Langevin Approach

,


I. INTRODUCTION
2][3] Typically, such combs are used for measurements of optical frequencies in metrology and sensing, having revolutionized these fields by providing unprecedented accuracy and enabling numerous innovative applications. 4,5][29][30][31][32] The active gain medium of the aforementioned lowdimensional SCLs provides a large third-order nonlinearity, which gives rise to a broadband four-wave mixing (FWM) process resulting in mode proliferation. 2,59][40][41][42][43][44][45] Stable and robust OFC operation requires a narrow beatnote, which is a measure for the amount of amplitude and phase noise of the a) johannes.popp@tum.de;https://www.ee.cit.tum.de/cph/home/b) jirauschek@tum.decomb lines.Noise accompanying carrier transport and spontaneous emission noise can therefore have a significant impact on the OFC formation and the performance of SCLs.
7][48][49][50][51] Recently, intensity correlations in QCL harmonic frequency combs (HFCs) have been experimentally investigated to develop a new generation of semiconductor devices for generating light with non-classical properties. 52Endowing commercial devices with outstanding quantum features would pave the way to practical and high-performance applications in the field of quantum networks, 53,54 including quantum computation, 55,56 quantum communication, 57 quantum metrology, [58][59][60] and quantum simulation. 61,62Notably, photonic systems are quite attractive for the investigation and employment of non-classical features, such as the generation of so-called quantum combs, [63][64][65] corresponding to non-classical states of light with multimode squeezed and/or entangled output.As the emergence of non-linear and non-classical features in SCLs is directly linked to the noise properties, 52,66 the development of low-noise SCL sources based on detailed simulations is an important prerequisite.
8][69] Here, the Bloch equations are used for simulating the evolution of the quantum system and its coherent light-matter interaction with the optical field in the active medium.Additionally, the optical field propagation is treated classically within Maxwell's equations, where the coupling with the quantum system arises from the macroscopic polarization term. 69The density matrix formalism can be extended and adapted by adding further quantized states in addition to the laser levels, and tunneling between states.
Within the Maxwell-density matrix equations, the dissipation processes can be modeled on a microscopic level, by including relevant mechanisms such as carriercarrier and carrier-phonon scattering.1][72] Typically, the dissipation in the Maxwell-density matrix equations is treated by phenomenological dephasing rates and electron lifetimes, which are either empirically chosen or estimated based on experimental data.Alternatively, advanced carrier transport models of different levels of complexity and numerical efficiency can be used to determine the scattering and dephasing rates self-consistently.][75][76] Since the semiclassical EMC model cannot cover quantum coherence effects, such as resonant tunneling across thick injection barriers, e.g., in THz QCLs, quantum corrections based on the density-matrix formalism have been incorporated. 77,780][81][82] Both are self-consistent 3D approaches not relying on empirical or fitting parameters as input and taking into account intrasubband processes between different kinetic energy states within a subband.In general, the NEGF approach contains the full quantum information of scattering and dephasing effects and is therefore superior to semiclassical methods such as the EMC method.However, the numerical simulation of the NEGF solution is expensive, if incoherent scattering within the self-consistent Born approximation is taken into account.To reduce the numerical costs, NEGF simulations are typically executed by neglecting electron-electron scattering beyond the Hartree approximation. 74ypically, the Maxwell-density matrix equations are treated in the so-called rotating-wave approximation where a coarser spatiotemporal grid can be used and thus the numerical load is greatly reduced.This approximation is however only valid for relatively narrowband (and not too strong) optical fields. 67,69,83On the other hand, for frequency comb sources a broadband spectrum with a high number of modes is desired, and the QCL even offers the potential for generating spectra extending over a full octave and beyond. 84In several studies, dynamical simulations of different nonlinear optical phenomena in nanostructured devices, e.g.QCLs 37,39,[85][86][87][88][89][90][91][92][93] and QD lasers, 35,[94][95][96][97] were conducted.[100] As pointed out above, this work aims to implement noise sources into mbsolve, enabling a realistic simulation of the noise characteristics in low-dimensional SCLs, which is for example important in the context of OFC operation.In a semiclassical framework, noise can generally be included by adding stochastic terms. 101,102hese are typically numerically implemented by using a pseudorandom number generator producing uncorrelated, Gaussian-distributed random numbers for every gridpoint. 103,104The resulting Maxwell-Bloch equations are then commonly solved numerically with the finite-difference time-domain approach. 103,104The magnitude of the stochastic terms can be systematically derived from the quantum Langevin equations, [105][106][107][108] which are then represented by equivalent stochastic c-number equations, 103,[109][110][111] i.e., evolution equations for operator expectation values with additional stochastic terms.3][114][115] Spontaneous emission obviously plays an important role in optoelectronic devices.While the resulting recombination can simply be included by nonlinear rate terms for the carrier occupations, 95,116 the noise contribution is not included in the Maxwell-Bloch model due to its semiclassical nature.This effect can however be considered in terms of a Gaussian white noise source in the optical propagation equation. 103,117,1180][121][122][123] By virtue of the fluctuation-dissipation theorem, a decay of populations, coherences, or the optical field is generally accom-panied by fluctuations, and a Maxwell-Bloch equation model which includes such decay-induced fluctuations has been presented. 104,124,125We extend this approach by including incoherent tunneling.From this, we can derive the semiclassical noise terms for our Maxwell-density matrix Langevin approach, which we then incorporate in our open-source tool mbsolve.Notably, an extension of the stochastic c-number approach to incorporate nonclassical effects has been discussed, 111,126 which might enable a direct investigation of nonclassical light generation in the presence of noise.
In detail, the paper is organized as follows: In Section II, we calculate the stochastic noise terms of the quantum Langevin equations and derive the generalized description of the Maxwell-density matrix Langevin equations for a quantum-optoelectronic structure such as a QCL.Our model is illustrated schematically in Fig. 1.
Here, the structure is described in terms of the optical field represented by the electric and magnetic field vectors E , H and the density operator ρ, which are coupled to each other by the interaction Hamiltonian ĤI .For the calculation of drift and diffusion operators in the quantum Langevin theory, we take into account the influence of various reservoirs in our structure.Regarding the quantum system, the reservoir interactions with the semiconductor host, which for example includes phonons associated with (longitudinal-and transverse-optical andacoustic) thermal lattice vibrations, lattice imperfections in the form of impurities (such as dopants), interface roughness or atomic disorder in alloys, as well as vacuum fluctuations arising from spontaneous emission are considered.For the optical field, the interaction with noise arising from thermal radiation (blackbody) entering the active waveguide from the cavity walls can be taken into account by external noise sources. 127In this paper, we rather focus on the fluctuations arising from the quantum system and dedicate the investigation of thermal noise influences, which can for example play a role in THz QCL active waveguides, to future works.In Section III, we give an overview of the simulation tool and the implementation of the noise terms and validate the model by presenting the simulation results for a superfluorescence setup. 104,128Furthermore, we present the simulation results for a THz QCL harmonic frequency comb and discuss the effects of noise contributions on the comb characteristics.The paper is concluded in Section IV.

II. THEORETICAL MODEL
In the following, we restrict ourselves to optoelectronic devices where the optical field can be well modeled using 1D Maxwell's equations.This for example applies to semiconductor lasers with longitudinally invariant waveguide geometries, where the 3D Maxwell's equations can be reduced to an effective 1D model. 69The focus lies on the inclusion of noise arising for example from spontaneous emission and fluctuations associated with the electron transport.First, we introduce the quantum Langevin equations using a simple three-level resonant tunneling QCL system.Here, the reservoir variables are eliminated and replaced by drift and fluctuation terms within the Heisenberg equation of motion.The quantum Langevin equations can be transformed into associated c-number Langevin equations, which are then used to derive the stochastic noise terms incorporated into the full-wave Maxwell-density matrix equation system.

A. The quantum Langevin equations
The quantum Langevin equations are introduced by using a simple three-level resonant tunneling QCL system as depicted in Fig. 2. 87,129 The QCL exploits optical transitions between quantized states in the conduction band of a quantum well heterostructure, where the properties can be controlled by quantum design rather than being determined by the bulk material.This not only applies to the gain and lasing wavelength, but also to the nonlinear optical properties such as FWM.Besides confinement provided by the quantum wells, another important quantum effect is tunneling through the separating barriers, which significantly influences carrier transport, in addition to the incoherent scattering-induced transitions due to phonons, crystal imperfections and electronelectron interactions. 74,77Regarding non-stationary QCL operation as is the case for OFC emission, coherent lightmatter interaction as a further quantum effect plays a significant role in the dynamic behavior, e.g., leading to Rabi flopping, 130 i.e. oscillations of the electron population between the upper and lower laser levels driven by the resonant optical field.Dephasing due to incoherent scattering has to be taken into account for a realistic description, as it greatly affects tunneling and coherent light-matter interaction.
For the structure shown in Fig. 2, the lasing transition occurs between the upper laser level |3⟩ and the lower laser level |2⟩.Depopulation takes place via level |1⟩ and electrons are injected from the depopulation level |1 ′ ⟩ of the adjacent period via resonant tunneling.The resonant tunneling across thick injection barriers in THz QCLs is treated within the tight-binding model. 74,76,78,129,131,132Here, the tunneling between a doublet of states at the thick injection barrier is described by the coupling strength Ω ij = −ℏ −1 ⟨i| Vext − Vtb |j⟩, with the extended conduction band potential Vext and the tight-binding potential Vtb .The coupling strengths Ω ij between the states |3⟩, |2⟩, |1⟩ within the active period are zero. 131n general, the QCL laser system is then described by the reduced system Hamiltonian 113,133,134

Period length
Lp where ĤF is the Hamiltonian of the optical field, Ĥ0 is the Hamiltonian of the quantum system with ∆ Vtb describing the coupling of electron states in two adjacent periods within the tight-binding model, and ĤI constitutes the interaction Hamiltonian between quantum system and optical field.Here, ℏ is the reduced Planck constant, ω 0 the single mode lasing angular frequency, â † (â) denotes the creation (annihilation) operator of the radiation field, ϵ i is the energy of level |i⟩ and ℏΩ 1 ′ 3 the anticrossing energy gap between levels |1 ′ ⟩ and |3⟩.The dipole coupling constant g can be written in terms of the dipole matrix element, µ z,23 = q⟨2|ẑ|3⟩, as 105,133 where ϵ r is the relative permittivity, ϵ 0 is the vacuum permittivity and V p is the volume of each quantum system associated with an active QCL period.The Heisenberg-Langevin equation of motion for an operator Âµ (t) reads as 105,107,134,135 Here, the drift operator Dµ (t) and fluctuation operator Fµ (t) account for the influence of the reservoirs on the system.[•, •] denotes the commutator [ X, Ŷ ] = X Ŷ − Ŷ X.For the drift operator Dµ we can under the Markovian approximation write 134,135 where w ± are the reservoir spectral densities and Qi is a function of system operators.For a detailed description and derivation of this theory together with the calculation examples for specific operators Âµ , we refer to Refs.134 and 135.
The reservoir average of the fluctuation operator vanishes, ⟨ F † µ ⟩ R = ⟨ Fµ ⟩ R = 0.The diffusion coefficient for a Markovian system is defined as and can be calculated by applying the fluctuationdissipation theorem.Here, the δ-function indicates the very short memory period of the reservoirs.The generalized Einstein relation for the calculation of the diffusion coefficient is given by 105,135,136 2⟨ Dµν (t From Eq. ( 3) together with Eqs. ( 1) and ( 4) the quantum Langevin equations for the three-level QCL quantum system can be derived.Therefore, we introduce the electron population operators σii = |i⟩⟨i| and the coherence operators σij = |i⟩⟨j|.The term σ32 â † describes the creation of a photon accompanied by an electron transition from the lower to the higher lying energy level, and σ23 â the annihilation of a photon accompanied by an electron transition from the higher to the lower lying energy level.At this point we drop these counterrotating energy non-conserving terms in the interaction Hamiltonian ĤI as in the commonly used rotating wave approximation. 105,133This simplifies the following calculations of the noise terms.A more complete calculation should also include more than one mode of the optical field in the system Hamiltonian.However, these concessions do not affect the form of the specific noise terms which are ultimately used in our simulations.The corresponding equations of motion are given by where κ is the cavity decay rate, ∆ ij denotes the energy separation between levels |i⟩ and |j⟩, τ −1 i = i̸ =j r ji is the inverse population lifetime, r ij,i̸ =j represents the scattering rate from level j to i and is the dephasing rate.Here, γ ij,p is the pure dephasing rate, which for QCLs mainly consists of elastic scattering contributions due to impurity and interface roughness. 77The equivalent equations for â † (t), σ32 (t), σ1 ′ 3 (t) and σ1 ′ 2 (t) are given by the Hermitian conjugates of Eqs. ( 7)(a)-(d).
Using Eq. ( 6), we can calculate the second-order correlation function relevant for the polarization operator as A detailed derivation of the diffusion coefficient 2⟨ D3223 (t)⟩ R can be found in Appendix A.
With the same procedure, we determine the other nonvanishing second-order correlation functions 134,135 is the number of thermal photons in the lasing mode at temperature T , where k B denotes the Boltzmann constant.

B. The c-number Langevin equations
In order to derive the stochastic noise terms for the semiclassical Maxwell-density matrix equations, the operator Langevin equations have to be converted into the associated c-number Langevin equations.
The quantum Langevin equation for the operator Âµ (t) in chosen order is given by where we make use of the commutation relation Â † to bring the equation into the chosen order.We use the superscript c to highlight that we have put the operators in chosen order.To explain this formulation in more detail, we use the fluctuation operator Fµ (t) as an example, but the following description holds for the other operators in the same way.For a chosen order Â1 , . . ., Âµ , we can write where the fluctuation operator F c µ in the chosen order is, of course, equal to the fluctuation operator Fµ in the original order.The associated c-number fluctuation term F c µ (A 1 , . . ., A µ ) is obtained using c-numbers A ν .By defining a linear chosen ordering operator Ĉ we can further indicate 134,137 where the operator Ĉ has the function of replacing each A ν by the corresponding operator Âν and bringing all terms into chosen order.
If we now convert the quantum Langevin equation into the equivalent c-number Langevin equation, we may write with L µ (t) being the coherent term corresponding to the commutation of Âµ (t) with the system Hamiltonian Ĥs , and D µ (t) denoting the drift term.Furthermore, by the use of Eq. ( 14) we obtain the c-number equation In analogy to the reservoir average in the operator case, we may write the c-number equation where we can make use of the following relation under the Markovian approximation: 105,134 2⟨D µν (t The diffusion coefficients in the c-number Langevin equations may differ from the ones in the quantum Langevin equations, as the c-numbers commute, whereas the operators do not.By requiring the equivalence of Eq. ( 16) and Eq. ( 6) in both c-number and quantum Langevin theory, it can be shown that in general By taking our chosen ordered operator representation Therefore, we derive the c-number second-order moments and obtain for the populations differing terms compared to the operator case: As an example, we provide a detailed derivation of the diffusion term from Eq.( 19)(a) in Appendix A. Additionally, we obtain diffusion coefficients absent in the quantum Langevin theory, e.g.
The complete diffusion matrix D(A, t) including all relevant cross-correlation terms of the threelevel QCL system with the c-number vector In literature, 138,139 it has been shown that a set of Ito stochastic differential equations (SDEs) can be derived for the given c-number vector and can serve as an efficient basis for numerical simulations.The equivalent Ito SDEs to the Langevin theory are given by (21)   where ξ(t) is a vector with real, independent Gaussian random numbers.Here, a semi-definite and symmetric diffusion matrix D(A, t) is required, which then can be factorized into the form 111,138,140 where the derived noise matrix B(A, t) is not necessarily symmetric.
To calculate the full noise matrix B(A, t) for the threelevel QCL system, we can divide the diffusion matrix D(A, t) into four different submatrices, where a correlation between the corresponding terms is identified.The given subvector A ν as well as the submatrices B ν (A ν , t) and D ν (A ν , t) are illustrated in Table I.Here, we include correlations between three states by taking into account a tunneling transition followed by an optical transition.This leads to a substantial extension of the initially derived quantum theory of propagation of nonclassical radiation in a two-level system 111 and is of essential importance for the description of quantum fluctuations in THz QCL systems, where electron transport across thick barriers is mediated by tunneling between closely aligned energy levels.A detailed symbolic derivation of the noise submatrices B ν (A ν , t) and the resulting noise matrix B(A, t) for the three-level QCL system can be found in the GitHub project mbsolve. 141y calculating the operator expectation value in the Schrödinger picture; we can demonstrate that the c-numbers representing the quantum system can be replaced by the density matrix elements ρ 23 , ρ 31 ′ , ρ 21 ′ , ρ 33 , ρ 22 , ρ 1 ′ 1 ′ , ρ 1 ′ 2 , ρ 1 ′ 3 , ρ 32 .The expectation value can be written as Furthermore, we can write the interaction Hamiltonian ĤI of the quantum system and the optical field as where the electrical field operator Êz is defined as Here, e z denotes the unit vector in z-direction.For devices in which the intraband transitions between quantized states occur within the conduction band, e.g. the QCL quantum well heterostructure, only the dipole matrix element µ z for the polarization in growth direction z is nonzero and relevant.
C. Generalized Maxwell-density matrix Langevin equations in 1D In the following, we derive the generalized Maxwelldensity matrix Langevin equations with additional microscopic fluctuation terms and characterize the influence of spontaneous emission noise on the optical field evolution.For the description of the coherent carrier dynamics and the incoherent relaxation processes, as well as the interaction with the classical optical field, the generalized full-wave Maxwell-density matrix equations constitute a compact semiclassical model.By combining it with the Langevin approach the microscopic noise characteristics can be fully taken into account.Here, the carrier dynamics in a SCL system are described in the density matrix formulation using the Lindblad equation (26)   which is coupled to Maxwell's equations in one dimension where D(ρ) is the dissipation superoperator, F(ρ) is an additional Langevin fluctuation superoperator and the other operators have their usual meanings.The permittivity is given by the product ε = ε 0 ε r , σ is the material conductivity, µ is the permeability, and the confinement factor Γ ∈ [0, 1] gives the spatial overlap of the transverse optical field mode with the quantum system.As we focus in this work on optoelectronic devices with invariant transverse field distribution, the reduction to a one-dimensional model for the optical propagation in the waveguide is justified. 69The Lindblad equation is the general form of a time-local and Markovian linear master equation for a quantum system, described by its completely positive trace-preserving density matrix, interacting with an environment.Obviously, the conventional Bloch equations, corresponding to a two-level system, describe the interaction of the laser levels with the optical field E z and constitute a special case of the Lindblad equation given in Eq. ( 26).The interaction with the environment is here modeled by scattering and dephasing rates, r ij and γ ij .Further levels can be considered in Eq. ( 26), and additional effects such as tunneling are included in the Hamiltonian.Moreover, quantum fluctuations are considered in the model given by Eq. ( 26) by adding a suitable Langevin fluctuation superoperator F. Maxwell's equations capture the optical propagation through the waveguide resonator, where the coupling with the quantum system is described by the macroscopic polarization P z,qm arising from the contributions of the dipole matrix elements.The expectation value of the dipole moment operator μz is calculated by averaging over a large ensemble of quantum systems within an adequate volume V p around the position z, and we can write for the macroscopic polarization where n 3D is the carrier number density.The two classical contributions, P z,class = ϵ 0 χE z and σE z , account for the polarization caused by bulk and waveguide dispersion as well as the material losses. 99inally, the update equations of the density matrix elements for the QCL laser system depicted in Fig. 2 can be written as Via the macroscopic polarization P z,qm , the quantum fluctuations added to the coherence term of Eq. (29a) have an influence on the evolution of the classical optical field.The quantum-mechanical fluctuations are based on the symbolic evaluations derived in Sec.II B and can be found in Appendix C. For the reduction to a twolevel system, we obtain similar noise terms as derived by Drummond and Raymer. 111However, unlike Drummond and Raymer we can assure the preservation of the physical properties of the density matrix, i.e., positive definiteness and unit trace.This is accomplished by a suitable choice of the submatrices B ν (A ν , t) depicted in Table I.

III. SIMULATION
The density matrix equations can be used to model the carrier dynamics in low-dimensional SCLs. 142In QCLs, for example, optical intersubband transitions occur between quantized energy levels in the conduction band of a multiple quantum well heterostructure.The quantized energy state with index i refers here to a subband and is characterized by a wavefunction ψ i .Due to the 1D-confinement in growth direction, the free electron motion is restricted to the in-plane directions and is in carrier transport simulations characterized by an in-plane wavevector k = [k x , k y ] T .132][143][144][145] In THz QCLs with small energy level spacings, we can typically neglect nonparabolicity effects and obtain k independent optical transition frequencies ω ij .In addition, Pauli blocking and Hartree-Fock renormalization effects, which occur in interband SCLs, play a subordinate role in QCLs due to the typically low doping levels. 69,146Since dynamic simulations of nonlinear optical effects, such as the generation of frequency combs, have to run over thousands of cavity roundtrips to achieve convergence, a complete 3D Maxwell-density matrix simulation is not feasible from a numerical point of view.If we take into account the k conservation of the dipole matrix elements, we can sum over k which yields density matrix equations of the form in Eq. (29).Under the assumption of moderate temporal variations of intersubband electron distributions ρ jj,k , averaged scattering rates r ij can be calculated by 69,76 where ρ 0 jj,k ′ is the steady-state electron distribution in the subband j.The k dependent scattering rates r jk ′ →ik can be extracted, e.g., from self-consistent DM-EMC simulations where they are calculated using Fermi's golden rule. 74In QCLs, the electron distributions within one subband can often be well described by Fermi-Dirac or Maxwell-Boltzmann distributions.Here, the characteristic subband electron temperatures can significantly exceed the lattice temperatures. 74,147The effective dephasing rates γ ij are typically obtained by averaging over the inversion between the subbands as 148,149 As mentioned above, the Maxwell-density matrix equations are commonly treated in the so-called rotatingwave/slowly varying amplitude approximation to reduce the numerical load, which is only valid for relatively narrowband (and not too strong) optical fields. 67,69,83However, as discussed in the introduction, in the framework of this paper a broadband frequency comb with many modes is highly beneficial, and QCLs even offer the potential for generating spectra extending over a full octave and beyond. 84None of the available open-source platforms are suitable for our purposes, mostly because they employ the rotating-wave approximation.Thus, we have developed our own open-source project mbsolve, allowing for numerically extensive simulations of multilevel systems based on the full-wave Maxwell-density matrix equations. 98,141In detail, the development of the codebase has been based on various principles.Here, the generalized Lindblad equation [Eq.( 26)] instead of the usual, quite restrictive two-level Bloch equation model is used.We have developed numerical methods that preserve physical properties, such as the complete positivity and trace preservation of the density matrix. 69,150This is especially important in the context of long-term simulations, as required for frequency comb modeling.Furthermore, a computational speedup is obtained by using parallelization techniques, 151 which is also especially important for long-term simulations of rather complex quantum systems, e.g.intracavity THz frequency comb difference frequency generation by mid-IR QCLs. 100 Our scientific software package mbsolve is developed following sustainable software engineering strategies and includes all common and essential best software engineering practices. 151,152It is based on C++ for performance reasons and features an easy-to-use Python interface facilitating the setup and active quantum system of the low-dimensional optoelectronic structures.The central part of the software is the mbsolve-lib base library, providing a framework for defining a simulation setup and the infrastructure to add solver and writer components.Importantly, mbsolve supports different numerical methods for solving the Lindblad equation, 150 as well as different parallelization techniques, e.g.OpenMP for shared memory systems.For a detailed package description, the reader is referred to Ref. 98.
As pointed out above, we extended our code base considering vacuum fluctuations due to spontaneous emission and fluctuations associated with the electronic transport. 103,104,124The solver class solver cpu fdtd has the method run, which executes the simulation loop by updating the magnetic field H y , electric field E z , density matrix ρ and polarization P z,qm for all spatial and temporal grid points in the Yee grid.This update procedure is explained in detail in Ref. 98.We add a new density matrix algorithm class, which is based on the Strang operator splitting method, 150,153 to account for fluctuations accompanying the electronic transport and vacuum fluctuations.The density algorithm class algo lindblad reg cayley qnoise contains the method propagate fluctuation, which calculates the fluctuations for an update step by adding the product of the fluctuation superoperator and the time interval ∆t • F(ρ) to the updated density matrix ρ n+1/2 .Here, we investigate active SCL gain media in lasing operation above threshold.In the simulations, we have to find a balance be-tween numerical efficiency and modeling accuracy.Most of the fluctuation terms given in Sec.II arise from the operator ordering when reducing the operator equations to c-number Langevin equations.As it was proven in literature, 104,122 these terms are negligible in the lasing regime above threshold with strong optical fields in the laser cavity.The fluctuation terms for a N -level system featuring diagonal elements F ii and off-diagonal elements F ij = F † ji can thus be significantly reduced for the numerical treatment and are described by where N cell is the number of carriers in one grid cell.The ξ 2,ij , ξ 2,ij and ξ 3,ij are real Gaussian random numbers and fulfill the correlation function For future applications, in which a more detailed fluctuation treatment would be beneficial, it might be necessary to extend our numerical model by additional noise terms derived in the previous section.Furthermore, we have implemented a class ic density random 2lvl, which represents random initial conditions for the common Maxwell-Bloch two-level system.As the dipole moment operators σ 12 , σ 21 and the atomic operators σ 11 , σ 22 do not commute, we have to take into account a non-vanishing initial stochastic value for the polarization term following the uncertainty principle. 104,154The tipping angle θ is obtained by drawing a random number from a Gaussian distribution with a standard deviation σ = 2N −1/2 cell and the angle ϕ in xy-plane is obtained by drawing a random number from a uniform distribution. 124

A. Superfluorescence and amplified spontaneous emission
The system was tested with a superfluorescence (SF) setup in a two-level configuration. 104,128This setup describes the spontaneous build-up of a macroscopic coherent dipole moment in an initially inverted system, resulting in a collective emission of a superfluorescent pulse.This behavior can be reproduced numerically within our mbsolve framework by simulating an ensemble of excited ions and using a dephasing time T number per cell N cell = 3 × 10 4 , transition frequency f = 477 THz, dipole length d = 6.875 × 10 −2 nm and equilibrium inversion w 0 = −1 are specified.We further investigate a device with length L = 7 mm using a grid discretization ∆x = 70 nm.The simulated SF pulse is illustrated in Fig. 3(a), and compares well with previous numerical and experimental findings. 104,128,155By increasing the collisional dephasing rate within the system, the SF pulse is significantly disturbed and gets broadened until the spontaneous buildup of the coherent dipole moment is prevented.For a dephasing time T 2 = 14.3 ps below the critical point, the SF pulse is replaced by amplified spontaneous emission (ASE).The increased noise amplitude accompanying the smaller dephasing time is crucial for the modeling of ASE, which cannot be reproduced otherwise.The ASE simulation results are presented in Fig. 3(b).
Furthermore, the degree of decoherence is studied using the quantity ρ 3 /ρ B , where is the length of the Bloch vector.When the dephasing time T 2 is high [Fig.3(a)], the population inversion ρ 3 is quickly depleted through the spontaneous buildup of the macroscopic dipole moment and the SF emission, which clearly surpasses the decay of ρ 1 and ρ 2 and results in a rapid drop of ρ 3 /ρ B .In the second case [Fig.3(b)] we have used a smaller dephasing time T 2 , which prevents the macroscopic dipole moment build-up and limits the radiative decay.This decoherence state indicates a very slow decay of ρ 3 /ρ B , which stays close to one.

B. Noise characteristics in THz QCL harmonic frequency comb emission
Concerning the experimental investigation of intensity correlations in QCLs, 52 we aim to characterize the noise properties of a self-starting THz QCL HFC setup.Here, the THz QCL active region is based on a homogeneous four-quantum well design with a diagonal transition. 156he charge carrier transport in the active gain medium at a bias of 50 kV/period is analyzed using our in-house Monaco framework, consisting of a Schrödinger-Poisson solver and a density matrix-ensemble Monte Carlo modeling tool. 74,76,77,157,158For an appropriate description of the physical properties, we consider five wavefunctions in the active quantum well heterostructure.Furthermore, one incoherent tunneling transition from the injector state into the upper laser level and one optical transition are specified for the quantum mechanical description of the QCL system in the dynamical simulation.The Python script forrer 2021 50mVperperiod.py with the simulation setup to start the mbsolve simulation is given in Listing 1 of Appendix D and is furthermore included in the GitHub repository. 141Here, all input parameters for the full description of the quantum system are extracted from self-consistent DM-EMC simulations.
In the following, we present simulation results for a 4 mm long double-metal THz QCL with a free spectral range (FSR) of 9.94 GHz.The intensity spectrum of the THZ HFC at 3.5 THz with a mode spacing of 5 FSR is illustrated in Fig. 4(a).The THz QCL emits a broadband HFC with a cavity repetition rate of 49.7 GHz.In Fig. 4(b), the temporal evolution of the intensity at the facet and the calculated instantaneous frequency are depicted.We can identify a regular field pattern, which shows a periodic repetition with five times the roundtrip time.Here, only the three strongest modes are involved in the temporal evolution of the instantaneous frequency, as their intensities are of similar magnitude and contribute most to the overall comb emission power.
To specify the degree of coherence of the obtained HFC and for comparison with the experimental findings, we investigate the radio frequency (RF) spectrum using an observation time window of 2 µs.The obtained simulation results are shown in Fig. 5(a), and the clear appearance of the harmonic beatnote proves the purity of the harmonic state.The linewidth is substantially below the numerical frequency resolution of 500 kHz, which is confirmed by the zoom on the extremely narrow harmonic beatnote in the inset of Fig. 5(a).In addition, we can identify sub-beatnotes, which arise due to the beating of the center mode with the sub-comb lines.These subcomb lines are generated by FWM processes, where the strong harmonic sidemodes act as pump modes and generate weak sidebands with a frequency spacing of 1 FSR from the corresponding pump modes.As can be seen in Fig. 4(a), the intensities of the sub-modes are at least ∼ 5 orders of magnitude smaller than those of the pump modes.
To further analyze the noise characteristics of the THz QCL HFC setup, we calculate the relative intensity noise (RIN) for the total output power and for the power of the five harmonic comb lines contributing most to the HFC emission.Here, the RIN spectrum can be calculated by where T denotes the simulation time, and P i is either the power of a specific mode i or the total power P all .By numerically filtering the electric field at the facet E facet (t) using a filter with a 3 dB bandwidth of 20 GHz, we can extract the temporal electric field components E i (t) of the individual modes.The RIN results are depicted in Fig. 5(b) for the total power P all and the power of the five central harmonic modes P i with indices i = 1 . . . 5. The total power RIN is around −180 dBc/Hz, while for the three central harmonic modes having a similar power a RIN around −155 dBc/Hz is calculated.For the remaining two weaker modes 1 and 5 a higher RIN is obtained.This is in very good agreement with the experimental findings of a three-mode mid-IR HFC QCL setup. 52For increasing power, the RIN of the sidemodes decreases to that of the central mode, while sidebands closer to threshold exhibit a noisier behavior.Furthermore, we identify an overlapping RIN for sidemodes featuring a comparable power level, which indicates a comparable noise level.
A similar result could be retrieved from the mid-IR HFC RIN measurements. 52.CONCLUSION In this paper, we have theoretically derived a c-number Langevin approach for a three-level quantum system from the non-classical operator description within the Heisenberg-Langevin equations.Our approach is an extension of the well-known two-level quantum theory by Drummond and Raymer, 111 where we additionally take into account incoherent tunneling injection into the upper laser level.Within the generalized Maxwell-density matrix Langevin equations we can ensure the preservation of the physical properties of the density matrix, i.e., positive definiteness and unit trace.Furthermore, by including the derived noise terms into our open-source simulation tool mbsolve, we can model the fluctuations accompanying electronic transport and spontaneous emission in the dynamical simulations of light-matter interaction in multilevel quantum optoelectronic systems such as QCLs and QD lasers.The simulation approach is tested using a superfluorescence setup, where we prove the validity of our implementation by obtaining an excellent match with previous experimental and theoretical results. 104,124,128,155Additionally, we have characterized the noise properties of a coherent THz QCL HFC setup and obtained a good match with experimental findings. 52,156Our modeling approach based on the generalized Maxwell-density matrix Langevin equations shows great potential for the theoretical investigation of intermodal intensity correlations in photonic devices and the development of low-noise integrated light emitters, also with regard to the generation of non-classical light.and by the ESA (Discovery EISI) Project 4000142337: "Simulation toolbox for unconditionally secure on-chip satellite quantum communication networks operating in the telecom wavelength range".

DATA AVAILABILITY
The code and all scripts for reproducing the data presented in this paper are available on GitHub.The diffusion coefficient in Eq. ( 9) is calculated by using the generalized Einstein relation from Eq. ( 6).The detailed calculation is given by In addition, we prove the difference in diffusion coefficients, which arises through the transition from operator to c-number Langevin equations.Therefore, we calculate the diffusion coefficient D 3333 (t) based on the theoretical description presented in Sec.II B. By the use of Eqs. ( 6) and ( 7)(e), we obtain Here, the terms which are underlined in red are not in the chosen order defined in Sec.II B. The commutation relations are used to bring these terms into chosen order, and by exploiting the level orthogonality similar to Eq. (A1) we derive With this, we can restructure Eq. (A3) as follows: Here, the additional terms resulting from the operator ordering are underlined in green.With the use of Eqs. ( 14) and (15) we can derive the corresponding c-number equation If we now require the equivalence of the left-hand sides of Eqs.(A5) and (A6), we end up with the diffusion coefficient from Eq.( 19)(a).
Appendix B: Complete diffusion matrix for a three-level THz QCL system The complete diffusion matrix for a three-level THz QCL system with incoherent tunneling injection is derived within the framework of the c-number Langevin theory and is given by .
Appendix C: Quantum mechanical fluctuation terms within the generalized Maxwell-density matrix Langevin equations The quantum-mechanical fluctuation terms for the three-level QCL quantum system are derived within the framework of the Langevin theory.In this paper, we have calculated the full diffusion matrix resulting from the c-number Langevin equations.Exploiting the positive semi-definiteness of the diffusion matrix, one can show that there exists a set of Ito stochastic differential equations equivalent to the Langevin equations.We can factorize the diffusion matrix to obtain a noise matrix that can be directly integrated into the Maxwelldensity matrix approach for numerical modeling of fluctuations in dynamical optoelectronic devices.With a suitable choice of the noise matrix, one can guarantee a completely positive trace-preserving update map for long-term simulations.For the three-level QCL system, the fluctuation terms will fully account for the influence of the reservoirs and the properties of the nonlinear coupling between QCL system and optical field, including the incoherent tunneling transition, and can be represented as follows:

Figure 1 .
Figure1.Schematic illustration of the coupling of a SCL quantum system and field system and the interaction with their associated reservoirs.

Figure 3 .
Figure 3. Simulation results for a superfluorescence test setup 104 in an initially inverted two-level system using the Maxwell-density matrix Langevin equations.(a) Cooperative emission characteristic of superfluorescence for the dephasing time T2 = 100 ps.(b) Amplified spontaneous emission pulse for the dephasing time T2 = 14.3 ps.

Figure 4 .
Figure 4. Maxwell-density matrix Langevin simulation results of HFC with a mode spacing of 5 FSR in a 4 mm long THz QCL device with a metal-metal waveguide at 80 K and for V = 50 mV/period.(a) Intensity spectra of the optical radiation at the facet.(b) Simulated instantaneous intensity at the facet and calculated instantaneous frequency from the Hilbert transform of the simulated electric field over a single roundtrip time (RT).

Figure 5 .
Figure 5. (a) Simulated RF spectrum of the THz QCL HFC setup with a clear beatnote signal at 49.7 GHz.Inset, zoom on the harmonic beatnote, indicating a narrow linewidth below the numerical frequency resolution (500 kHz).(b) Calculated RIN spectra associated with the total power P all (blue) and the modal power Pi of each of the five harmonic modes contributing most the HFC emission.The colors of the individual RIN spectra correspond to those of the individual comb lines in Fig. 4(a).

Table I .
Division of the diffusion matrix into submatrices Dν (Aν , t) and the corresponding c-number subvectors Aν and noise submatrices Bν (Aν , t).In order to preserve the physical properties of the quantum system description, we have to interpret noise matrices Bν (Aν , t) differently for occupation and coherence terms.The differing matrix expressions for the coherence terms are highlighted here in red.