Multi-Domain Modeling of Free-Running Harmonic Frequency Comb Formation in Terahertz Quantum Cascade Lasers

Optical frequency comb (OFC) emission in quantum cascade lasers (QCLs) is highly attractive for applications in metrology and sensing. Recently, coherent OFC mode-locking with large intermodal spacing was demonstrated in QCLs. These self-starting harmonic frequency combs (HFCs) show highly phase-stable operation and promise interesting perspectives towards optical or even quantum communication. Here, we present a detailed multi-domain modeling approach for the numerical study of QCLs. Our theoretical characterization is divided into stationary carrier transport simulations, based on the ensemble Monte Carlo method, and dynamical simulations of the light-matter interaction, based on multi-level Maxwell-density matrix equations. We investigate the influence of the chosen eigenstate basis on the gain spectrum and present self-consistent simulation results of stable HFC operation in a double metal terahertz QCL. In our simulations, the studied QCL gain medium shows self-starting harmonic mode-locking for different bias and waveguide configurations, resulting in a mode spacing of up to twelve times the cavity round trip frequency. Furthermore, we characterize the spectral time evolution of the coherent HFC formation process, yielding the spontaneous build-up of a dense multimode state which is gradually transferred into a broad and clear harmonic OFC state.

Fig. 1.Schematic of a double metal QCL waveguide with THz OFC emission through the laser active region facet.Here, the difference between a fundamental OFC and a HFC is illustrated by the varying mode spacings within the emitted comb spectrum.For the fundamental comb the mode spacing corresponds to the roundtrip rate f rt , while for the HFC spectrum, we obtain a mode spacing of n × f rt , where n is an integer multiple of one.
Recently, OFC emission in a novel QCL operating regime was detected, featuring harmonic modes separated by multiples of the cavity roundtrip (rt) frequency.In both the mid-IR and THz spectrum, harmonic frequency comb (HFC) states of varying orders were obtained for different bias points and waveguide geometries [16], [17], [18], [19], [20], [21].Potential applications of HFC QCL setups arise in the area of wireless terahertz communication networks [17], THz imaging [19] and quantum optics experiments [11], [21].Free-running HFC formation in THz QCLs is mainly based on double-metal waveguide configurations.A schematic of such a THz HFC QCL setup is illustrated in Fig. 1, where the active gain medium is sandwiched into a double metal waveguide, and the outcoupled coherent THz light is either specified by a fundamental or a harmonic OFC spectrum depending on the order of the intermodal spacing.Recently, stable self-starting second-order HFC emission could be demonstrated for a single-plasmon THz QCL with a cavity length of 15 mm [22].Effective coherent emission of harmonic combs is further achieved by active mode-locking [23], [24], [25], or can be controlled and manipulated through external forcing, e.g.macroscopic defects [26], [27], [28], external cavities [29] or optical seeding [30].
In this work, we provide a substantial numerical study of THz QCLs based on a diagonal transition design, which have shown self-starting fundamental and harmonic mode-locking [16], [61], [62], [63].Our results show the self-assembling of highorder HFCs in a free-running device, where the mode-locking results from the internal laser dynamics.In contrast, previous fully numerical studies have observed HFCs with large combline separation (more than 2 × f rt [44]) only in the presence of external forcing (macroscopic defects [26], [27], [28], active modulation [24], or external cavities [29]), to the best of our knowledge.On the other hand, analytical models have been derived and predict gain for largely detuned side modes [20], [58], or provide an explanation based on mean-field treatment [64], [65].For the characterization of the active gain medium, we use our in-house open-source monacoQC framework, featuring a versatile wavefunction solver [35], [66], [67] and a stationary carrier transport model, which is based on the DM-EMC method [35].For the investigation of the dynamical behavior in the THz QCL, we use the open-source solver tool mbsolve for the full-wave generalized Maxwell-DM equation system, which is composed for the modeling of light-matter interaction in multi-level quantum systems without invoking the rotating wave approximation (RWA) [40], [68].All input parameters for the dynamical simulations can be self-consistently extracted from the quantum cascade device modeling tool monacoQC.In Section II a detailed description of the used multi-domain simulation approach for THz HFC QCLs is presented.As pointed out in several theoretical analyses of QCL gain media, the chosen wavefunction basis has a significant impact on the gain characteristics [35], [36], [67].In order to analyze the influence of the eigenstate basis set on the laser output state, we calculate wavefunctions in extended [66] and localized [35], [36], [67] basis, for a THz QCL design that showed HFC generation in experiment [16].The solutions of the Schrödinger-Poisson (SP) equation system are depicted in Fig. 2, at a bias of 50 mV/period, showing a significant difference in terms of wavefunction localization.The influence of the basis set on the gain characteristics and the harmonic mode-locking behavior is discussed in Section III-A.Additionally, we illustrate frequencyresolved results of the unsaturated gain profile, obtained by Fig. 3. Overview of the monacoQC framework [69].The project consists of five modules: optimizer, setup, mbsolve, solver and results.The module ext.solver is not included in the GitHub repository.
pump-probe dynamical simulations, and long-term simulations of the four quantum well QCL in harmonic comb operation.In Section III-B, we investigate the influence of the applied bias and the waveguide geometry on the harmonic ordering and characterize the spectral time evolution of the self-starting harmonic mode-locking mechanism.The paper is concluded in Section IV.

II. MULTI-DOMAIN SIMULATION APPROACH
For spatiotemporal simulations of self-starting THz HFC QCL emission, we use a self-consistent multi-domain approach.To model the light-matter interaction in a THz HFC QCL, our open-source Maxwell-DM tool mbsolve is used [49], [50], [68].All input parameters required for a complete description of the quantum system are extracted from our in-house open-source monacoQC framework [50], [69], a quantum cascade (QC) device simulation tool [70], [71].In the following, the monacoQC framework with the base library for QC device specifications and the charge carrier transport simulation results are explained in more detail.Furthermore, the mbsolve framework for dynamic full-wave Maxwell-DM simulations of coherent light-matter interaction in quantum optoelectronic structures is introduced.

A. Device Description and Carrier Transport
The monacoQC framework is an open-source tool, consisting of an object-oriented Matlab-based device engineering tool for QC structures and including different modules [70], [71].A schematic overview is depicted in Fig. 3.A setup library for the description of the QC active region and the simulation parameters is provided.The class device includes all important information about the active gain medium and the waveguide.Here, the active region is formed by one period consisting of the layer sequence with barrier and well materials.An object of class layer includes the information about the layer width, the doping, interface orientation and an object of class material.
Here, the abstract class material is divided into three subclasses binary, ternary and quaternary providing important material information, e.g.bandgap energy E g and effective mass m eff of the III-V compound semiconductors and their ternary and quaternary alloys [72].The most relevant material systems for QCL device engineering are InGaAs/InAlAs in the mid-IR regime and AlGaAs/GaAs in the THz regime.All required parameters of the waveguide model, e.g. the cavity length l c , the linear loss term α 0 , the facet reflectance R and the confinement factor Γ, are summarized in the class waveguide.Other relevant environmental simulation parameters, e.g.applied bias V and temperature T , are specified in the class scenario.
Taking into account all important parameters in our device description, the corresponding wavefunctions ψ i and system Hamiltonian Ĥs are calculated within a Schrödinger-Poisson solver library.The numerical approach used for solving the Schrödinger equation is based on the transfer matrix method [31], [66], which is implemented in the corresponding class tm_solver.Additionally, a class poisson_solver takes care of the mean-field treatment of the electron-electron interaction and determines the additional electrostatic potential.Accounting for an appropriate quantum mechanical description, different boundary conditions for the solution of extended states (boundary_condition_ext) [31], tight-binding states (boundary_condition_tb) [35], [36], or EZ-states (boundary_condition_ez) [67] can be specified.The solver method solve solves iteratively the Schrödinger and Poisson equation until convergence is reached and returns an object of class eigenstates and of class conduction_band, respectively.
Here, eigenstates comprises the quantum mechanical description with wavefunctions, Hamiltonian and averaged effective masses of the subbands, while conduction_band contains the total potential consisting of conduction band profile and band bending effects due to space charges.Additionally, a Bayesian optimization module (BO_sim) based on Gaussian process regression is integrated into the monacoQC framework for developing high-performance QC devices in a structured manner [71], [73], [74].Here, the input space can be chosen arbitrarily and spans from layer width variations in the active QC heterostructure to changes in the material composition or the doping density.Selecting a suitable figure of merit (merit_fun), which is evaluated based on carrier transport simulation results with specific post-processing tools, results in an efficient and flexible optimization procedure.
Our charge carrier transport simulations are based on the EMC method, a self-consistent 3D approach taking into account all necessary intrasubband processes [31].The scattering-induced transport of electrons between quantized levels is modeled semiclassically by evaluating the Boltzmann equation.Here, the scattering rates for phonon scattering, electron-electron scattering, impurity scattering, interface roughness scattering and alloy scattering are considered.The scattering rates are calculated self-consistently at the beginning of the simulation based on Fermi's golden rule.Therefore, the quantized energy levels and corresponding wavefunctions have to be provided.An ensemble of electrons is distributed over the different states characterized by the subband index i and in-plane wavevector k.According to the Monte Carlo method, the QCL system will converge to the stationary solution by stochastically evaluating the scattering events.To overcome the lack of quantum coherence and quantum mechanical dephasing mechanisms, the EMC model was extended by essential quantum correction terms [35].The DM-EMC approach includes elements of the density matrix formalism to account for incoherent tunneling, which plays a significant role in THz QCLs.To account for the tunneling transport between tight-binding wavefunctions ψ i , ψ j in the DM-EMC approach across thick injection barriers, we calculate the Rabi frequencies Ω ij = ψ i |V − V tb |ψ j / based on the extended and tight-binding potential V and V tb [35], [36].Within the monacoQC framework EMC simulation results can be integrated and processed using the following classes: r scattering_rates: includes a map of all considered scatter- ing mechanisms and corresponding scattering rates.
r dephasing_rates: includes a map with all considered pure dephasing contributions and level broadenings with respect to the kinetic energy of each subband.
r carrier_distribution: contains the subband occupation probabilities and the electron distribution within each subband.The stationary carrier transport simulation results, as needed for the dynamical solver, can be summarized in the backend class mbsolve_sim.The resonant tunneling and required optical transitions for mid-IR (ω mid−IR ) and THz (ω THz ) frequencies have to be identified.Using various helper methods, a reduced model, which describes the quantum system adequately, is retrieved.Furthermore, the methods calc_gain, calc_suscept_2 and calc_gvd for calculations of the gain characteristics, the second order susceptibility |χ (2) | and the group velocity dispersion (GVD) are provided in the class mbsolve_sim.The class method generate composes a Python script of the quantum mechanical description in mbsolve syntax, which can be directly used as input script for a mbsolve simulation via its Python interface [68].The quantum mechanical description comprises the level occupations ρ ii , the system Hamiltonian Ĥs with eigenenergies E i and anticrossing energies Ω ij , the dipole moment operator dz , the dephasing rates Γ ij and the scattering rates r ij .For the one-dimensional dynamic Maxwell-DM simulations, the energy-resolved dephasing rates are simulated within the EMC approach and have to be averaged over the population inversions of the involved subbands, as it is described in detail in [35] and [75].The base library together with the presented carrier transport results in extended and EZ-configuration can be found on GitHub [69].

B. Generalized Maxwell-Density Matrix Equations in 1D
In our modeling tool mbsolve, the light-matter interaction is treated by the full-wave Maxwell-DM equations.In this semiclassical framework, a quantum mechanical description of the electron dynamics and a classical description of the optical field are used.The density matrix ρ and its interaction with the environment is described by the Lindblad equation, the general form of a time-local and Markovian linear master equation for quantum systems, which can be written as The interaction Hamiltonian − dz E z describes the interaction of the quantum system with the optical field E z .The dissipation superoperator D(r ij , Γ ij , ρ) displays the interaction with the environment and is here described by a function of the scattering rates r ij , the dephasing rates Γ ij and the density matrix operator ρ.Multiple levels can be considered in (1), and additional effects such as tunneling can be included in the system Hamiltonian Ĥs .The optical propagation through the waveguide resonator is described classically within Maxwell's equations.For optical field propagation in waveguides with invariant transverse field distribution, a one-dimensional model is justified and helps to reduce the numerical load [40].The time evolution of the magnetic field H y and the electric field E z is given by where ε 0 ε r,∞ is the product of the vacuum permittivity and the relative permittivity at high frequencies, σ is the material conductivity and i ∂ t P i z,class is the sum of polarization terms accounting for GVD in dispersive media [59], [76].Light-matter interaction is included by taking into account the macroscopic polarization P z,qm , which is obtained by summing over the expectation values of the microscopic dipole moments of the quantum systems, where n 3D is the charge carrier density.Equations ( 1) and ( 2) constitute a semiclassical model for the simulation of the spatiotemporal dynamics in quantum optoelectronic devices.Detailed information about the macroscopic Maxwell-DM equation system and the mbsolve tool can be found in [40], [68].

III. RESULTS AND DISCUSSION
In the following, we present simulation results for a selfstarting THz HFC QCL setup [16], [62].The device consists of a four quantum well active region, which is embedded in a doublemetal waveguide featuring high facet reflectance.A schematic of the THz HFC QCL setup is illustrated in Fig. 1.As described above, we investigate the influence of the chosen eigenstates basis on the gain spectrum and present self-consistent simulation results of stable HFC operation with a mode spacing of 4 × f rt .Furthermore, we investigate the influence of the applied bias and waveguide geometry on the harmonic ordering.Finally, we characterize the spectral time evolution of the self-starting harmonic mode-locking mechanism.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

A. Influence of Wavefunction Basis Sets on the Coherent Operation in THz HFC QCLs
For an adequate optical and electrical description of the THz QCL active gain medium, we take into account five wavefunctions in one active period.The QCL system is investigated using extended and localized wavefunction configurations within the SP solver.Electron injection into the upper laser level (ULL) is governed by resonant tunneling from the injector state (INJ) of the adjacent period [blue rectangle in Fig. 2(b)] and can be appropriately described within the tight-binding model [35], [36].Additionally, we take into account further coherences arising from closely aligned energetic levels by applying an EZ-transformation [67].Here, eigenstates separated by an energy of less than 5 meV are summarized within a multiplet of states.These subsets of eigenstates are diagonalized with respect to the dipole moment operator [67].A transformed triplet of states within the investigated THz QCL configuration is schematically illustrated in Fig. 2(b) by an orange rectangle.For the characterization of the full dynamical range extending from single-mode operation to fundamental and harmonic comb states, carrier transport simulations over an extended bias range have been conducted for the two configurations, respectively.In the following, we analyze the unsaturated gain behavior of the four quantum well QCL gain medium.
In Fig. 4, the peak gain and center frequency for the bias range 40 mV/period to 60 mV/period are illustrated.We can identify two operating regimes.In the lower bias range < 51 mV/period, both the peak gain and center frequency for the localized states tend towards higher values compared to the extended states.Due to the strong localization of ULL and lower laser level (LLL) in space, the transition exhibits a smaller lifetime broadening for the optical transition, resulting in narrower gain curves with increased peak values.The center frequency is also higher, as the EZ-transformation shifts the eigenenergies of the triplet of wavefunctions [LLL, depopulation level 1 (DP1), and 2 (DP2)] towards a collective energy level, which is below the LLL eigenenergy of the extended states.In the second bias interval > 51 mV/period, the peak gain value for the extended states collapses, while the optical frequency continues to increase successively.The wavefunctions depicted in Fig. 2 are obtained for an applied bias of 50 mV/period, which is close to the intersection of the two intervals.The energy position of the two wavefunctions INJ and LLL at the injection barrier becomes crucial within the second interval [green rectangle in Fig. 2(a)].By increasing the bias, the anticrossing energy gap is gradually reduced until the crossing point of eigenenergies is reached.This crossing process for the extended states is accompanied by a wavefunction extension for the INJ and ULL levels across the injection barrier.The increasing impurity dephasing rate due to the ionized donors in the injection well results in a broadened gain curve with a reduced peak value.In contrast, we have a strong spatial localization of these wavefunctions in the localized states configuration as illustrated in Fig. 2(b) (blue rectangle), leaving the gain characteristics unaffected by artificial broadening mechanisms.At 55 mV/period, we identify a small frequency drop for the localized states configuration, as the energy gap between the lower laser level LLL and the depopulation level DP1 exceeds 5 meV and the EZ-transformation for this subset of wavefunctions is no longer conducted.The wavefunctions of the LLL have a similar shape in both configurations.Besides, the peak gain value in the extended states configuration gradually grows for increasing bias values, as the two wavefunctions INJ and ULL separate again in energy.The presented carrier-transport results are publicly available and can be found in the monacoQC GitHub repository [69].From the stationary carrier transport simulation results, we can extract the quantum mechanical description of the QCL active gain medium and use it as input for the dynamical Maxwell-DM solver.To specify the influence of the tunneling transition INJ → ULL on the spectral gain profile, we perform a dynamical gain analysis within the Maxwell-DM framework.As incoherent tunneling is absent in the extended state basis, but this is the major injection mechanism to the ULL in the diagonal transition design, we do not consider this configuration in the following dynamical simulations.For the modeling of the gain characteristics within the mbsolve framework, we seed a weak Gaussian pulse on the left facet and record the electric field at the middle of the 4 mm waveguide to measure the light amplification within the active gain medium [46], [77].A linear field loss term α = 6.5 cm −1 extracted from COMSOL simulations for a copper-copper (Cu-Cu) double-metal waveguide [16] is included in the simulation setup and has to be taken into account for the interpretation of the dynamical simulation results.The obtained unsaturated gain spectrum is depicted in Fig. 5 for the localized states configuration.We obtain a strong gain lobe for the optical transition accompanied by a weak side gain lobe at the lower frequency, which results from the incorporation of the tunneling transition INJ → ULL in the system Hamiltonian.In recent years, the linewidth enhancement factor has become an important parameter for the description of QCL dynamics [43], [57].It gives a measure for the coupling between amplitude and phase fluctuations in semiconductor lasers with fast gain dynamics and was initially introduced in the linewidth theory of semiconductor laser, going beyond the Schalow-Townes limit [78].Microscopic models of varying complexity, e.g. based on EMC [79] and NEGF [80] simulation approaches, were used to calculate the LEF factor in QCLs.In our fullwave Maxwell-DM equations, the LEF factor is introduced by the gain asymmetry resulting from the tunnel coupling between the INJ and ULL states, while assuming that other contributions, e.g.due to nonparabolicity, play a secondary role in this THz QCL.To extract the LEF, we carry out pump-probe simulations, with a Gaussian pulse as seed, where we iteratively increase the pulse power until the gain clamping condition is reached.The LEF, also known as α-factor, can be computed at the center frequency ω c using Here, χ(ω) is the complex susceptibility of the QCL gain medium.By slightly changing the pulse power in the dynamical simulations, we can calculate a LEF in the QCL design.An α-factor of ∼ −0.1 is obtained for a bias of 50 mV/period.This is in good agreement with experimental and theoretical findings [29], [44], [80], [81], [82], [83].Furthermore, we conduct long-term simulations with an end time of t e = 1000 ns characterize the particular mode-locking behavior.We obtain harmonic comb emission with a mode spacing of 4 × f rt , located around the center of the gain maximum, determined by the pump-probe simulations.The intensity spectrum is depicted in Fig. 5, illustrating HFC emission with seven clear and narrow comb lines.The appearance of the associated harmonic beatnote, and absence of the fundamental beatnote and its intermediate harmonic beatnotes, prove the purity of the harmonic state.A linewidth substantially below the numerical frequency resolution of 1 MHz is detected in the inset at 35.7 GHz of Fig. 5.

B. Analysis of Harmonic Mode-Locking Regimes in THz QCL Waveguides
In this section, we proceed with dynamical simulations using the QCL system based on the localized states configuration.From the simulation results in Fig. 5, we find that the contribution of the gain lobe caused by injection tunneling to the HFC mode proliferation is negligible.This reflects in the small magnitude of the α-factor, as discussed above.Furthermore, all HFC modes lie within the gain lobe arising from the optical transition between ULL and LLL.Therefore, we omit the Rabi energies corresponding to injection tunneling in the quantum mechanical description and include additional scattering rates for the carrier injection.
In the following, we present a HFC analysis of different bias values and waveguide geometries.Furthermore, we analyze the temporal evolution of the HFC modes.
1) Applied Bias: Initially, we set the facet reflectance to 0.5, which arises from the high impedance mismatch in doublemetal waveguides [16], [46], [84].Here, the facet reflectance can vary between R = 0.5 − 0.9 depending on the waveguide dimensions relative to the wavelength [84], [85].The chosen reflectance value is thus at the lower limit.Various strategies for reducing the waveguide loss and attenuating the lasing of higher-order transverse modes in THz QCL setups are discussed in the literature [86].Harmonic comb emission in THz QCLs was recently demonstrated by the use of a low loss Cu-Cu waveguide [16].In contrast to gold-gold waveguides, where the HFC gradually builds up from a single mode, the HFC in Cu-Cu waveguides emerges spontaneously and alternates over the operating bias range with fundamental comb and high-noise states.The lower waveguide loss in Cu-Cu waveguides leads to a wider dynamic range and increased intracavity fields.The latter becomes important for the FWM processes, which are proportional to the cube of the electric field intensity.Harmonic comb states at bias points far above the laser threshold are detected.We conduct simulations for different bias values within the operation regime presented in Section III-A and obtain HFC emission spectra of varying harmonic orders.The distinct coherent comb regimes for three representative bias points are illustrated in Fig. 6.At smaller bias values of around 48 mV/period we obtain a dense/fundamental comb [Fig.6(a) top].A clear beatnote at 9.05 GHz is detected.The radiofrequency (RF) spectrum is shown in Fig 6(b), also containing higher beatnotes arising from the beating of wider-spaced modes.A regular field pattern consisting of a periodically repeating waveform at each roundtrip is retrieved, where a snippet of two roundtrips is depicted in Fig. 6(c).In recent years, the spectral and time domain properties of free-running THz QCLs have been extensively investigated both theoretically [42], [44], [46], [87], [88] and experimentally [8], [9], [86], [89], [90].Self-starting OFC emission in THz QCLs is characterized by a simultaneous frequency-and amplitude-modulated signal.In the upper panel of Fig. 6(c), parts of the electric field profile exhibit a flat behavior, which can be explained by the help of the time-varying instantaneous Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.frequency profile, also illustrated in Fig. 6(c).In the first time interval of ∼ 20 ps, the electric field is governed by a single frequency mode, which acts as the dominant one in the intensity spectrum in Fig. 6(a).By applying a higher bias of 50 mV/period a harmonic comb with a mode spacing of 4 × f rt ≈ 36.12GHz is obtained [see Fig. 6(a) middle], similar to the HFC results in Section III-A.The associated RF spectrum is presented in Fig. 6(b).In the time domain, a periodic waveform with four repetitions per roundtrip is visible [Fig.6(c)].As expected, the corresponding instantaneous frequency features the same periodicity as the amplitude.On closer inspection, it is clearly visible that only the strongest modes within a ∼ 10 dB power range of the optical spectrum contribute significantly to the temporal behavior of the instantaneous frequency.The quasilinear chirp is in good agreement with recently published experimental results for THz QCL OFCs [8], [86].By further increasing the bias towards regimes of decreasing unsaturated gain values, the separation of the locked modes reduces again [see Fig. 6(a Here, the linear chirp trend is even more pronounced as for the HFC of fourth order.By increasing the bias further, we obtain a noisier state, which however reveals tendencies towards higher order harmonics.The system is then apparently not stable enough to evolve into a clear HFC state.
The obtained results are in good agreement with experimental findings [16], [18], where OFC states of varying order are retrieved and the transition from a dense state to a harmonic comb state and back for bias sweeping in a given interval is documented.Furthermore, we can reproduce the bias-dependent alternating between different comb states, as it is documented for Cu-Cu waveguides [16].Obviously, a change in bias results in a change of center frequency and gain shape, and also the dynamical properties of the quantum system, e.g. the level lifetimes and dephasing rates, vary with the bias.It is worthwhile mentioning that in our simulations only one optical transition is considered.This implies that self-starting harmonic mode-locking can also emerge without the assistance of an asymmetric gain profile, which has been interpreted in earlier theoretical approaches to be essential [16].
2) Waveguide Geometries: In the experimental studies, the THz HFC QCL setup was characterized for different waveguide configurations with variations in the cavity length and width [16].In order to emulate the changed impedance mismatch for different waveguide widths in our simulations, we vary the facet reflectance R for the THz HFC QCL setup at an applied bias of 50 mV/period.In addition to the results of Fig. 6, we present simulation results for facet reflectances of 0.32 and 0.8, which represent a single-plasmon and a second double-metal waveguide with different geometries, respectively.For the reflectance of 0.32, a HFC emission spectrum with a mode repetition rate of 3 × f rt is obtained [see Fig. 7(a)].Interestingly, two twin lobes around the center frequency evolve.A similar shape of the optical spectrum can be derived from experiments in which a HFC state with a mode spacing of 5 × f rt is obtained at a drive current of 800 mA [16].In Fig. 7(b), the emission spectrum for a reflectance of 0.8 is presented, revealing a HFC state with eighth harmonic order.The changed intracavity dynamics appear to have an emphasized effect on the harmonic mode ordering.For increased reflectance values, we obtain higher intracavity fields, affecting both the FWM processes and the spatial hole burning (SHB), i.e., the arising spatial gain modulation in a linear cavity with standing waves.For a noisy operating regime with tendencies towards higher order harmonics, an increase of the field strength within the cavity could result in a stabilization to a coherent state.We therefore execute a simulation of the THz QCL at a bias voltage of 55 mV/period and facet reflectance of 0.8.A stable HFC with a mode spacing of 12 × f rt is retrieved.It is obvious that the free-running THz QCL system is very sensitive to the different environmental parameters, highly affecting the dynamic behavior of the quantum system and thus the optical output.
3) Spectral Time Evolution: To gain an intuitive understanding of the spectral time evolution of the harmonic mode-locking in the THz HFC QCL setup, we investigate the OFC dynamics in the first 50 ns of the Maxwell-DM simulations.Here, we present simulation results for the three QCL configurations with facet reflectances of 0.32, 0.5 and 0.8, (see Fig. 8).The optical field within the THz QCL cavity is triggered by spontaneous emission events and starts from dense multimode lasing.From that, it evolves into the HFC state.The temporal characteristics of the three configurations are similar, whereby small deviations are evident.In Fig. 8(a), the HFC formation for the QCL setup with a reflectance of 0.32 is illustrated.Here, a strong center mode only changes weakly, while the energy of the emerging harmonic sidemodes increases gradually.The transition from the dense multimode state to a broad and clear HFC state with a mode repetition rate of 3 × f rt is clearly recognizable.In the QCL configuration with a facet reflectance of 0.5, we observe a more inert behavior during the HFC generation [see Fig.  the center mode is rather weak, with two stronger side lobes, or the center mode is strong, with a strong decay of intensity in the side modes.
The QCL gain medium features a large third-order nonlinearity χ (3) [5], [6], which leads in combination with SHB to a broadband FWM process and strong mode proliferation.In the case of harmonic mode-locking, the interaction of strong harmonic pump modes reduces the gain in their close spectral environment and thus acts as parametric suppression of adjacent sidemodes [19].In addition, a parametric enhancement of widely detuned modes is present and can clearly be observed in Fig. 8(c), showing the spectral time evolution of the HFC with eighth harmonic order for the QCL setup with a facet reflectance of 0.8.During the temporal evolution of the HFC, the adjacent sidemodes still get populated by degenerate FWM processes [46].As long as the parametric suppression is too weak, the sidemodes can survive and no clear harmonic state evolves.In the case of the HFC state with the three dominant modes being equally strong and neighboring [see Fig. 6 second row and Fig. 8(b)], the parametric suppression competes with the seeding of the sidemodes, which results in a more inert behavior.A supplementary movie showing the evolution of the twelfthorder HFC state at a bias of 55 mV/period and for a reflectance value of 0.8 is provided.Here the above-mentioned dynamics in the initial phase of HFC evolution can be clearly seen.

C. Discussion
With our multidomain simulation approach, we can replicate the experimental results of stable self-starting HFC generation of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
higher orders in THz QCLs.Furthermore, we provide a detailed analysis of the active gain medium and investigate the influence of the applied bias and waveguide geometries on the HFC formation.Harmonic states of different orders are retrieved for the different configurations and the experimentally documented behavior of alternating fundamental and harmonic comb states over the operating bias range is replicated in the simulations.To better categorize the results, we give a brief overview of existing simulation results of HFC emission in QCLs based on the MB equation system and compare them with our approach.Furthermore, we discuss ideas for the next steps to shed light on the physical origins of harmonic mode coupling in QCLs.
Numerical simulation results of QCL harmonic comb generation using a three-level Maxwell-density matrix model were presented in [58].For the generation of HFCs of varying order, their simulation was seeded with an electric field consisting of a weak central mode and very weak sidemodes at higher orders of the intermodal spacing.Here, they could demonstrate, that the harmonic states are sustained after switching off the electric field seeding and can describe the typical QCL dynamics when interacting with a strong coherent laser field.The seeded HFC states were found to be only stable for a few thousand roundtrips, however, degrading to dense/fundamental states in the long run.More recently, simulation results of fundamental and harmonic comb generation in QCLs embedded in a Fabry-Perot cavity were published based on the ESMB equation approach [42], [44].An extensive study on dense and harmonic regimes is provided and experimental features such as the simultaneous amplitude and frequency modulation behavior were reproduced.However, self-starting HFC operation could only be demonstrated in second order using the ESMB equations.
The simulation approaches used in the aforementioned publications for the simulation of self-starting HFC generation in QCLs are based on the common rotating wave/slowly varying envelope approximation, generally used to reduce the numerical effort associated with the fast field oscillations [37], [40].Optical interference effects such as SHB in Fabry-Perot resonators significantly influence the QCL dynamics and are taken into account by an inversion grating [40], [45], [52], [55], [57], [64], [91].Further extensions towards third-order polarization grating were proposed and considered in studies of Risken-Nummedal-Graham-Haken instabilities in QCLs [56].The ESMB equations additionally include a non-zero LEF to account for certain semiconductor specifics, such as asymmetric gain and dispersion profiles.Within our multidomain simulation approach, we utilize the full-wave Maxwell-DM equations and thus intrinsically account for all required effects of self-starting HFC formation, e.g.higher-order polarization and population gratings, nonlinearities and off-resonant dynamics.We have extended this approach by the inclusion of stochastic noise terms derived from the quantum Langevin equations to characterize the noise properties of the THz QCL HFC setup.The simulation results are illustrated in [50] and show good agreement with the experimental results.In fact, this realistic spontaneous emission noise is not amplified to macroscopic distortions, but intrinsically kept at low values.Therefore, the intermodal beatnote is still observed to be extremely narrow (at the numerical frequency resolution limit ≈ 1 MHz), indicating a high purity of the observed comb states.Higher-order nonlinearities and grating terms were omitted in earlier HFC simulations employing the RWA, as their importance for HFC formation was considered negligible.Due to the complexity of our simulation approach, it is not possible to abstract the model.The higher-order grating terms and nonlinearities, which we will refer to as beyond-RWA effects, appear to have an essential influence on self-starting HFC formation in QCLs.To provide insights into the governing mechanisms of self-starting harmonic mode-locking, the crucial beyond-RWA effects have to be identified and characterized using simplified models.Possibly, the adiabatic elimination of specific terms may provide additional insights.This approach has already been used for the prediction of purely frequencymodulated combs in mid-IR QCLs [57] or recently for the description of passive single-pulse mode-locking in THz QCLs with distributed saturable absorbers [52].
Our numerical approach also shows great potential for the modeling of harmonic comb operation in mid-IR QCLs.We have recently used this full-wave framework to model THz OFCs created by difference frequency generation in mid-IR QCLs, in good agreement with experimental data [70].Numerically reproducing the entire optical spectrum spanning from the THz to mid-IR regime, is only possible by considering full-wave Maxwell-DM simulations.Broadband mid-IR HFC spectra with large intermodal spacing up to 26 times the free spectral range of the laser cavity were experimentally demonstrated [17], [20], pushing towards the intrinsic bandwidth limitations of MB models invoking the RWA.For such scenarios, our here-discussed approach might not only be applicable, but even advantageous.

IV. CONCLUSION
In summary, we provide a substantial theoretical analysis of self-starting THz HFC emission in QCLs.The simulation results of harmonic combs featuring different orders are conducted using a self-consistent multi-domain modeling approach, coupling a DM-EMC carrier transport simulation tool to a dynamical Maxwell-DM solver.The investigation of different basis states reveals their influence on the gain characteristics and the formation of HFC states in THz QCLs.By comparing extended to localized states in the here simulated QCL structure, we identify an underestimation of the injection transition into the upper laser level in the extended configuration.The resulting small gain appears to be unsuitable for the formation of a HFC in this setup.In contrast, the more pronounced gain spectrum in the localized states basis with a dominant optical transition ULL → LLL alleviates the formation of a dense FC comb, which evolves over time into a harmonic comb.Additionally, we investigate the influence of the applied bias and waveguide geometry on the HFC formation and obtain different dynamical regimes with varying harmonic orders.Finally, we characterize the spectral time evolution of the self-starting harmonic mode-locking mechanism, providing new insights into the HFC creation process.In the future, we aim to pin down the physical mechanisms enabling harmonic mode-locking by reducing the model complexity.

Fig. 2 .
Fig. 2. Calculated conduction band profile and probability densities of the investigated THz QCL structure based on GaAs/Al 0.15 Ga 0.85 As for a temperature of 80 K.The layer sequence (in nm) of one active period is 4.1/18.3/5.8/10.6/2.9/11.5/3.9/9.2.Barriers are in boldface, for the underlined layers a doping density of 2.3 × 10 16 cm −3 (n-type) is assigned and the applied bias is 50 mV/period.Comparison of Schrödinger-Poisson solutions based on extended (a) and localized (b) states [66], [67].

Fig. 4 .
Fig. 4. Monte Carlo simulation results of the THz QCL gain medium using extended (orange diamonds) and localized states (blue squares).Simulated unsaturated peak gain values (a) and center frequencies (b) are shown as a function of the applied bias V .The bias window is divided into two operating regimes, which are highlighted by shades of green and blue.The specific characteristics of the different configurations within the two regimes are explained in the text.

Fig. 5 .
Fig. 5. Spectrum of the Maxwell-DM simulation setup at a bias of V = 50 mV/period.The individual comb lines follow well the unsaturated gain curve, proving the validity of the self-consistent multi-domain simulation approach.Inset: Zoom on the radiofrequency beatnote with a numerical frequency resolution of 1 MHz.

Fig. 6 .
Fig. 6.Comparison of three coherent regimes in a 4 mm long THz QCL device with a metal-metal waveguide at 80 K [16].Results of a fundamental frequency comb (top panel, 48 mV/period), harmonic combs of 4 × f rt (middle panel, 50 mV/period) and with a mode spacing of 2 × f rt (bottom panel, 54 mV/period) are presented.(a) Intensity spectra of the optical radiation for the three regimes exhibiting clear mode spacings.(b) RF spectra with beatnotes at the fundamental (9.05 GHz), fourth harmonic (36.12 GHz), and second harmonic (17.83 GHz) roundtrip frequency f rt , respectively.(c) Time-resolved electric field amplitudes and calculated instantaneous frequency from the Hilbert transform of the simulated electric field.Arrows indicate the individual waveform period for the three regimes.
) bottom].At a bias value of 54 mV/period a beatnote at the second harmonic roundtrip frequency of 17.83 GHz is detected [see Fig. 6(b)].The harmonic time signal repeats twice per roundtrip and is together with the instantaneous frequency shown in Fig 6(c).

Fig. 7 .
Fig. 7. Intensity spectra of two harmonic states in a THz QCL at an applied bias of 50 mV/period and a temperature of 80 K. Simulation results for a facet reflectance of 0.32 (a) and 0.8 (b) are presented.
8(b)].The sub-comb lines with a mode spacing of 1 × f rt from the strong harmonic pump modes only slowly fade out over time.This operation regime corresponds to the middle panel of Fig. 6(a), where we identify three strong central modes with similar intensity.In contrast, the two configurations with R = 0.32 and R = 0.8 (Figs. 7, and 8(a), (c), respectively), do not show similarly strong modes at the center frequency.Either

Fig. 8 .
Fig. 8. Time evolution of the the spectral intensity I for self-starting harmonic mode-locking in three THz HFC QCL waveguide configurations at an applied bias of 50 mV/period and a temperature of 80 K. Simulation results for a facet reflectance of 0.32 (a), 0.5 (b) and 0.8 (c) are presented.