Theory of Hybrid Microwave–Photonic Quantum Devices

The Maxwell–Lindblad transmission line (MLT) equations are introduced as a model for hybrid photonic and microwave‐electronic quantum devices. The hybrid concept has long enabled high‐bandwidth photodetectors and modulators by using waveguide structures supporting microwave–optical co‐propagation. Extending this technique to quantum‐engineered lasers provides a route for greatly improved short‐pulse and frequency comb operation, as impressively demonstrated for quantum cascade lasers. Even more, the self‐organized formation of coupled microwave and optical oscillation patterns has recently enabled photonics‐based ultralow‐noise microwave generation. The MLT model provides a suitable theoretical description for hybrid quantum devices by combining optical and microwave propagation equations with a Lindblad‐type approach for the quantum active region dynamics. A stable numerical scheme is presented, enabling realistic device simulations in excellent agreement with experimental data. Based on numerical results and analytical considerations, it is demonstrated that the functionality of the investigated hybrid quantum devices does not only rely on local coupling effects, but rather benefits from microwave‐optical co‐propagation, settling an ongoing debate. The MLT equations provide the basis for systematic device development and extension of the concept to other quantum‐engineered hybrid devices based on, for example, quantum dot or interband cascade structures, as well as for simplified analytical models providing additional intuitive insight.


Introduction
The ability to generate periodic optical waveforms [1] with comblike spectra [2] directly from lasers has revolutionized metrology and sensing, and has also impacted many other fields such as materials processing and communications.This is achieved DOI: 10.1002/lpor.202300461 by active or passive mode-locking, resulting in periodically amplitude-modulated (AM) and/or frequency-modulated (FM) optical fields. [3,4]From an application point of view, electrically pumped semiconductor lasers have many significant advantages such as compactness, robustness, and cost-effectiveness, and are available for a broad range of frequencies.[8][9] This hybrid device concept has recently been taken to another level by applying it to quantumengineered mode-locked sources in the form of quantum cascade lasers (QCLs).These are a special type of nanostructured semiconductor laser which employs the quantized electronic states in a multiple quantum well structure as laser levels. [10,11]Thus the lasing wavelength and gain profile, as well as the nonlinear optical and microwave properties, can be custom-tailored.This enables a wide range of functionalities in the mid-infrared (MIR) and terahertz (THz) range such as broadband frequency comb generation, [12,13] while for conventional semiconductor lasers comb operation remains challenging. [14]ince QCLs typically have cavity lengths in the mm range, the repetition rate of the optical waveform, and hence the spacing between two adjacent frequency comb lines, is on the order of 10 GHz, which is well-suited for the implementation of hybrid optical and microwave functionalities. [15]The key lies in combining the QCL active region with a metal-metal waveguide structure, formed by contact layers on the top and bottom of the active region. [16]In this way, simultaneously an optical waveguide and a microwave transmission line is obtained, which can be designed to support co-propagation of the optical field and microwave-electronic signal. [17]he device structure is schematically illustrated in Figure 1.In such a setup, the generation of ultrashort mode-locked THz pulses, which is highly relevant for numerous application fields, [18][19][20] was demonstrated by exploiting co-propagation of the externally applied microwave modulation signal and the optical wave. [21,22]A similar strategy was exploited for the generation of broadband and low-noise THz frequency combs, [23] which are essential in several application domains. [24]Notably, this approach is fundamentally different from previous injectionlocked short-pulse [25] and frequency comb [12] setups, where the applied microwave signal only affects the injection region of the waveguide. [21,26]A particular appeal of the hybrid concept lies in the potential to exploit self-seeded microwave signals in a freerunning device, thus eliminating the dedicated external circuitry and electronics required for active modulation.This leads to a new class of devices where both the optical and microwave dynamics can unfold in a self-organized manner.The microwaveoptical coupling is here provided by the quantum-engineered active region, featuring both optical and low-energy transitions. [27]here has been experimental evidence that microwave selfseeding in microstrip QCLs contributes to comb stabilization, [23] and may even enable mode-locking at the second or higher-order harmonics in free-running QCLs. [28]Based on the same concept, a hybrid terahertz-microwave source has been demonstrated, offering a new route for the generation of extremely low-noise microwave signals of up to hundreds of GHz. [27]This paradigm of hybrid device design has recently culminated in an approach for independently tailoring the terahertz and microwave properties based on a planarized waveguide geometry. [15]ince the above-demonstrated hybrid device concept relies on the coupled propagation of the optical radiation and the injected or self-generated microwave signal, a moderate microwave loss across the waveguide length is crucial, and also the role of microwave-optical phase matching has been emphasized. [17,21,27]hus, the proper choice of waveguide type and careful microwave design of the overall device are of critical importance.While microstrips in the form of metal-metal waveguides have been identified as the most suitable type for hybrid QCL devices, [17] it has been pointed out that the microwave loss may still be a limiting factor. [29][32][33] Thus, there is an ongoing debate if the observed phenomena are really due to microwave-optical co-propagation or can be explained by conventional modulation effects, with the measured microwave beat-notes along the waveguide arising purely from local stimulation by the laser radiation. [26]Clarification of the underlying mechanisms requires a detailed understanding of the processes involved, and a comprehensive description of the hybrid opticalmicrowave dynamics is all the more important for the targeted design of mode-locked hybrid devices.In this context, it has long been pointed out that a suitable modeling approach is still missing. [8,34]Ultimately, the significance of a comprehensive theoretical model goes far beyond QCLs, since further types of quantum-engineered lasers, such as interband cascade, [35] quantum dot/dash (QD), [36,37] and not least the more conventional quantum well lasers, [38] have been shown to possess favorable properties for AM or FM mode-locking. [39]This makes them highly suitable candidates for extending above discussed hybrid concepts, such as the self-organized coupled microwave-optical dynamics in free-running devices, to other spectral regimes.
The lack of a comprehensive hybrid model is resolved in this paper.A common approach for modeling the semiconductor laser dynamics are the Maxwell-Bloch (MB) equations.Combining Maxwell's equations for the optical field propagation with a density matrix (DM) approach describing the coherent electron-light coupling, this semiclassical theory constitutes a basic and very widely used model for the coupled dynamics of an optical field interacting with an ensemble of two-level quantum systems. [40,41]44][45] On the other hand, by extending the original two-level model to multiple quantum states and introducing a generalized system Hamiltonian, the resulting framework consisting of coupled Maxwell and DM equations constitutes a powerful tool for the detailed dynamic modeling of quantum-engineered optoelectronic devices. [46]We refer to this approach as Maxwell-Lindblad (ML) equations since the DM dynamics is here described using the Lindblad formalism. [46,47][50] Modeling of steady-state mode-locking requires long-term simulations of the laser dynamics.Full-wave 3D simulations are too computationally expensive, and thus the model equations are commonly reduced to the waveguide propagation direction. [46]The microwave propagation is not contained in the MB/ML equations, but can typically be regarded as quasi-transverse electromagnetic for microwave-photonic structures and thus also be reduced to a single spatial dimension in the form of a distributed equivalent-circuit transmission-line representation. [51]Co-propagation effects in traveling-wave devices have been modeled by coupling this description to basic optical propagation equations, [52,53] and it was also proposed as early as 2004 to include the transmission line dynamics into the theoretical description of mode-locked QCLs. [54]Along those lines, we introduce an adequate model for hybrid quantum devices by combining the transmission line equations with a ML-type approach featuring a suitable Hamiltonian, and refer to the resulting model as Maxwell-Lindblad transmission line (MLT) equations.Although such an approach may seem straightforward, several issues must be considered which may have impeded the development of a self-consistent model including the microwave dynamics: i) The strongly frequency dependent microwave loss Laser Photonics Rev. 2023, 17, 2300461 of the transmission line must be adequately considered to obtain a realistic description; ii) the coupling between quantum active region and transmission line has to be self-consistently included in the model in form of a hybrid interaction Hamiltonian; iii) for the numerical evaluation of the resulting partial differential equation set, a stable numerical scheme is required.
In detail, the paper addresses the following issues: In Section 2, we introduce a self-consistent model for the coupled optical, quantum electronic, and microwave dynamics in the form of the MLT equations.Section 3 contains an analytical discussion of microwave propagation in transmission line structures.In Section 4, we present a stable and efficient numerical scheme for the MLT equations, enabling the direct simulation of experimentally realized devices.Based on the obtained results, the presence and significance of microwave-optical co-propagation effects is demonstrated.The paper is concluded in Section 5.

Model
The modeling of mode-locked laser operation requires longterm simulations over many hundreds or thousands of cavity roundtrips. [48]To minimize the numerical load, the model is usually reduced from three spatial dimensions to a single one, involving only the propagation coordinate x and time t. [46]As opposed to fully 3D approaches, [55] it is then not practical to describe both the optical and microwave field by a single set of Maxwell's equations for several reasons: In particular, the optical and microwave frequencies differ by more than two orders of magnitude even for THz QCLs, and also the transverse mode profiles deviate from each other (all the more so if the optical and microwave confinement are designed separately), [15,17] resulting in quite different effective parameters in the 1D model.Furthermore, the optical field is typically treated in the slowly varying envelope approximation (SVEA) in the propagation equation, along with the rotating wave approximation (RWA) in the density matrix equation, to reduce the numerical effort associated with the fast oscillations. [46]On the other hand, the microwave propagation in transmission line structures is usually described in the framework of a distributed-element model in terms of voltage and current rather than field quantities, establishing a direct connection to related characterization techniques and measurement results [17] and readily allowing for adaption to different transmission line configurations. [52,53,56]Thus, we model the microwave and optical dynamics by separate, but coupled propagation equations.

Microwave Propagation
The microwave propagation on the waveguide structure is described in terms of the voltage u(x, t) and current i(x, t) (see Figure 1).Here, it is advantageous to write u(x, t) = u 0 + u Δ (x, t), where u 0 denotes the DC component of the input voltage at the position of the wire bond, defining the operating point of the device, while u Δ contains the modulation.The transmission line equations can then be written as R ′ , L ′ , and C ′ have their usual meaning as transmission line resistance, inductance, and capacitance per unit length in the direction of the transmission line (x coordinate in Figure 1), while J(x, t) denotes the current density flowing through the active region of width w.The inclusion of the waveguide microwave losses, represented in Equation ( 1) by R ′ , is crucial for realistic modeling. [26]Furthermore, R ′ has a strong frequency dependence which must be taken into account, [17] leading to additional time derivative terms in Equation ( 1). [57]This model can be readily adapted to more complex transmission line configurations by adding further distributed components. [52,53,56]In the following, we assume that u represents the voltage drop across the active region.The boundary condition at x = 0 is given by the applied pump current or voltage, which can be constant or modulated.Assuming a waveguide of length L with open-circuit termination, the boundary condition at x = L is i = 0.

Optical Propagation
In the MB equations, typically the SVEA is invoked to reduce the numerical load. [46]The optical field is then written as with the complex slowly varying envelopes E ± (x, t) of the forward and backward propagating resonator field components, optical carrier frequency  c , and propagation constant of the guided mode .E represents the optical field strength in polarization direction, for example, for QCLs parallel to the growth direction z of the heterostructure since the other field components do not interact with the quantum active region. [46]The field propagation equation is in the SVEA given by [46] v −1 Here, a denotes the waveguide power loss coefficient, v g is the group velocity, and p ± contains the polarization induced by the quantum system (see Equation (A6) in the Appendix).Also the background group velocity dispersion (GVD), characterized by the coefficient  2 , is included in Equation ( 3), since it can considerably affect both ultrashort pulse generation and frequency comb operation. [35,42,43,58,59]S sp ± represents spontaneous emission, which is in the simulation used to start the laser action, and furthermore adds fluctuations to the optical field.This process is numerically implemented as spatially distributed random noise source. [60]The boundary conditions for reflection coefficients r 1,2 at the left and right facets are given by E + = r 1 E − at x = 0 and E − = r 2 E + at x = L. [46]

Quantum System Dynamics
In the MLT model, the dynamics of the quantum active region is described by the DM formalism.Here, we do not restrict ourselves to the conventional two-level approximation which is often too restrictive for realistic device modeling, but rather use a Lindblad-type multilevel approach, which may be supplemented  by additional equations for the incoherent scattering-induced transport. [46]For inclusion of the coupling to the microwave transmission line, our model has to be formulated in dependence of bias voltage, [49,61] and yield the current density as an input for the transmission line equations [48,61] (see Figure 2).The temporal evolution of the density operator ρ is described by the evolution equation where ℏ denotes the reduced Planck constant.Within the 1D model, ρ(x, t) denotes the density operator of a representative quantum system at the macroscopic position x in the laser waveguide.The Hamiltonian for the quantum active region, for example, a representative QCL period or QD featuring N relevant states |n⟩ with energies E n and associated transition frequencies The unperturbed Hamiltonian is given by Ĥ0 = ∑ n E n |n⟩⟨n|.In general, the system Hamiltonian may also contain off-diagonal elements, which are included in ĤT .For example, many QCL designs feature thick injection and/or extraction barriers, across which electron transport is mediated by tunneling between closely aligned energy levels, | mn | ≪  c .[64] Furthermore, ĤI = −E(d mn |m⟩⟨n| + d nm |n⟩⟨m|) is the light-matter interaction Hamiltonian in dipole approximation for an optical transition between two levels m and n.Here, d mn represents the dipole matrix element.It is given by d mn = −e⟨m|ẑ|n⟩ for QCL intersubband transitions where e denotes the elementary charge, while it depends on the conduction and valence band Bloch functions for interband transitions. [46]If the system contains several relevant optical transitions, each of them is represented by a corresponding Hamiltonian.For optical transitions close to resonance, | mn | ≈  c , the RWA is commonly employed by inserting Equation ( 2) and discarding the rapidly oscillating terms.Likewise, the interaction with the microwave field can be represented by a Hamilto-nian of the form ĤI , where the modulation field strength across the active region of thickness d is given by E = u Δ ∕d.Due to the comparably slow dynamics in the microwave range, the RWA is here not applied, and the relevant transitions have frequencies correspond in first-order perturbation theory to the change in eigenenergies ΔE n caused by a bias field u Δ .Writing the system Hamiltonian as ĤS = Ĥ + ĤRWA , Ĥ describes the coherent dynamics of the system due to Ĥ0 , ĤT and the interaction with the microwave field.Thus, the matrix elements of Ĥ are given by Here, E n,0 are the level energies at the operating point u 0 .In a tight-binding framework, the tunneling transitions have practically vanishing dipole moments, that is, d mn ≈ 0 for Ω mn ≠ 0. [65] ĤRWA is the light-matter interaction Hamiltonian treated in RWA.][68] For the resulting DM equations in RWA, see Equation (A1)-(A5) in Appendix A. Dissipation effects are in Equation ( 4) considered by the term  t ρ| inc .Here, we use the usual Lindblad-type model, [46,47] where incoherent transitions from a state |m⟩ to |n⟩ with a scattering rate r mn result in terms  t  nn | inc = − t  mm | inc = r mn  mm , and analogously for the scattering in the reverse direction with rate r nm .In this way, incoherent transport between the energy levels of the quantum system can be taken into account.The associated dephasing gives rise to terms  t  mn | inc = − mn  mn ,  t  nm | inc = − nm  nm , where the dephasing rate is modeled as  mn = (r m + r n )∕2 +  p,mn , with the total outscattering rate r m = ∑ s≠m r ms and pure dephasing contribution  p,mn =  p,nm .In principle, bias-dependent scattering and dephasing rates r ij (u Δ ) and  ij (u Δ ) may be used in Equation (A1)-(A5). [49]This effectively introduces a time dependence into the Lindblad operators, which is compatible with the Lindblad approach. [46,69]The inclusion of bias-dependent scattering rates may, for example, be necessary for modeling scattering-assisted QCL designs, where the electron injection/extraction is mediated by scattering rather than tunneling.In this case, the current transport, and thus the optical gain, strongly depends on the scattering cross-sections, [70] while the behavior of resonant-tunneling designs is governed by tunneling, already contained in the bias-dependent Hamiltonian of Equation (6).Additional terms beyond the Lindblad framework, such as nonlinear rate terms, may be required to apply the model, for example, to QD lasers (see also Appendix A).
To obtain a self-consistent model which only requires the device specifications and material parameters rather than using fitting or empirical parameters, we couple the density matrix equations of Appendix A to carrier transport simulations at the operating point, providing the scattering and dephasing rates. [49]he carrier transport simulations contain the full dependence on the electron in-plane wave vector k, thus also providing the subband electron distributions, [71] while this dependence has been integrated out in the dynamical DM equations to reduce the numerical load. [49]This approach enables long-term simulations and seems especially appropriate for operating conditions where the bias and current density modulations are only a few percent of the DC values; thus, it is well-suited for the examples discussed in Section 4. Its validity (without the inclusion of transmission line effects) has been confirmed for various QCL structures by comparison to experimental results. [49,50,72][77] While this adds versatility to the model, the numerical cost increases significantly, thus impeding long-term dynamical simulations.

Coupling
The density matrix, optical propagation, and transmission line equations are coupled as illustrated in Figure 2, resulting in a closed model.The coupling induced by the optical field and the polarization in the ML (i.e., optical propagation and DM) equations is obtained in the usual manner via the polarization term Equation (A6) contained in Equation (3), and the light-matter interaction terms in Equation (A1)-(A5) resulting from the Hamiltonian ĤRWA .The quantum system and the microwave transmission line are coupled via u Δ in Equation ( 6) and the current density J in Equation ( 1), which is computed from the density matrix. [46]For example, Equation (B1) in the Appendix can be used under certain assumptions for QCLs, and a model for the current density in QD media can be found in literature. [61]While SHB in the active region in principle causes short-scale spatial modulations of J, these do not transfer to u Δ since they are strongly damped by the transmission line and, thus, are here neglected for both J and u Δ .

Analytical Considerations
As discussed in Section 2.4, the transmission line couples to the quantum active region via the voltage, and modulation thereof in-troduces a variation of the optical gain.Thus, for actively modelocked semiconductor lasers, microwave-optical co-propagation effects can only be expected if a considerable voltage modulation is obtained along the whole waveguide structure.For QCLs, the fulfillment of this premise has been debated in literature [26] and can be investigated based on simple analytical considerations.As already discussed, for a realistic description of transmission lines, it is vital to consider the dependence of R ′ () (and possibly other parameters) on the frequency  = 2f , [17] resulting in an operator function R ′ (i t ) in Equation ( 1).In the following, we assume that the DC resistance is negligible, R ′ ( = 0) ≈ 0. [17] Besides using u(x, t) = u 0 + u Δ (x, t), we insert the ansatz i(x, t) = i 0 (x) + i Δ (x, t) into Equation ( 1), where u 0 and i 0 denote the DC components defining the operating point of the device, while u Δ and i Δ contain the (zero time-average) modulation signal.Assuming a small modulation voltage u Δ and neglecting nonlinearities as well as dynamical effects in the quantum active region, we can express Jw in Equation ( 1) as Jw ≈ J 0 w + G ′ u Δ , where J 0 and G ′ are the active region current density and differential conductance per unit length at u 0 .With the boundary condition i 0 (L) = 0 as in Section 2.1, the DC current is then given by i 0 (x) = J 0 w(L − x).For the time varying component, we correspondingly have i Δ (L) = 0. Furthermore, we assume at x = 0 a pump voltage modulation u Δ = u  cos(t).The transmission line then forms a lossy half-wavelength resonator with the fundamental resonance frequency  r = c∕(Ln = r ), as used in experiment to match the optical mode spacing. [17,21]Here, n  is the microwave effective index at a frequency , and c denotes the vacuum speed of light.The analytical solution for a given modulation frequency  can be obtained by inserting u Δ = ∑ ± ℜ{u ± exp(±i  x − it)} (where u + and u − are the amplitudes of the forward and backward propagating components) into Equation (1), along with a corresponding ansatz for i Δ .This yields with Z ′ () = R ′ () − iL ′ and Y ′ () = G ′ − iC ′ the standard results for the propagation constant   = i(Z ′ Y ′ ) 1∕2 and characteristic impedance Z  = ±u ± ∕i ± = (Z ′ ∕Y ′ ) 1∕2 .The modulated voltage is obtained as Separating the propagation constant into real and imaginary parts yields   = n  ∕c + ia  with the damping coefficient a  .Exemplarily, we determine the transmission line parameters for the THz ultrashort-pulse QCL structure [21] as follows: We use R ′ = 4.5 × 10 −2 (f ∕Hz) 1∕2 Ω m −1 as measured for a similar waveguide structure, [17] and G ′ = 80 S m −1 extracted from measurements. [21]The Hammerstad formula yields L ′ = 144 nH m −1 for the given microstrip geometry. [78]Finally, C ′ = 1.13 nF m −1 is determined such that n  ≈ 3.9 is obtained at f = 12 GHz. [21]The ultrashort-pulse QCL is driven close to  r . [21]The ratio between the modulation amplitude at x = L and x = 0 is then obtained from Equation (7) as 1∕ cosh(a  L) = 0.26 for a  = 657 m −1 and L = 3.1 mm, that is, a significant voltage modulation is achieved across the device despite the high microwave loss.
A more detailed assessment of the loss-induced damping is obtained by evaluating the roundtrip-averaged voltage modulation.For (fundamental) mode-locked operation, the optical roundtrip time is locked to the modulation period.We again assume a small modulation u Δ as above (which is realistic considering the strong microwave attenuation in the external circuitry [79] ) and fast gain dynamics.The latter applies, for example, to QCLs where the gain recovery time is typically much shorter than the roundtrip time. [21,80]The modulated part of the gain coefficient can then approximately be written as g Δ (x, t) = g u u Δ (x, t), with some proportionality factor g u .To evaluate the modulation in the vicinity of the optical pulse, a co-moving reference frame is introduced via the retarded time coordinate t ′ = t − xn  ∕c.The modulated part of the roundtrip gain is then G Δ (t ′ ) = exp[2Lg u ūΔ (t ′ )], with the roundtrip-averaged modulation voltage Inserting Equation (7), we obtain at This voltage waveform gives rise to a net gain window centered at t ′ ≈ 0, as required for active mode-locking. [3]For a  = 657 m −1 , the roundtrip-averaged modulation amplitude is reduced by less than 50% as compared to its maximum value u  ∕2 at a  = 0.These basic considerations already show that despite considerable microwave damping in this THz QCL waveguide, the modulation signal not only acts locally at the injection point, rather, propagation along the waveguide induces strong interaction with the co-propagating optical pulse, increasing the effective modulation depth, and thus enabling devices with improved modelocking [21] and modulation bandwidth. [17]he conventional scenario, that is, modulation of a waveguide section with length L m < L, [25] can also be investigated with above model.[83] The modulation bias is then given by u Δ = u  cos(t) for 0 ≤ x ≤ L m , and u Δ = 0 otherwise.Inserting u Δ into Equation (8), we obtain at  =  r the roundtrip-averaged modulation voltage For small modulation section lengths, L m ≪ L, the effective modulation strength increases linearly with L m as intuitively expected, ūΔ ≈ (L m ∕L)u  cos(t ′ ).Furthermore, we obtain ūΔ → 0 for L m approaching L, indicating that short-pulse generation by active mode-locking is not possible in this case, as also discussed in previous theoretical works. [82,83]However, as also pointed out in this context, [82] neglecting the spatial dependence is not realistic for long modulation sections, explaining the discrepancy between this result and the full solution, Equation (9), which is in line with the experimental findings. [21]n a simplified model, the back-action of the optical field onto the microwave transmission line signal can be estimated in the following way: The net number of generated photons per unit volume and time at a given cavity position is ṅp = (g − a)I∕(ℏ c ). Neglecting SHB, I = I + + I − with I ± =  0 cn 0 |E ± | 2 ∕2 denotes the optical intensity, where  0 and n 0 are the vacuum permittivity and optical refractive index.Furthermore, g is the power gain coefficient.Assuming each net photon emission event contributes an electron to the current flow, we obtain a light-induced change in current density where L p denotes the length of a single QCL period.Taking into account the voltage dependence of the gain, we have g = (g 0 + g u u Δ ).The simplest way to include the fast gain dynamics is in terms of an instantaneous saturation model g = (g 0 + g u u Δ )∕(1 + I∕I s ) with the saturation intensity I s ; furthermore, gain filtering due to the finite spectral gain bandwidth must be considered. [3]owever, the full light-matter interaction dynamics can only be captured by a quantum mechanical model, such as the DM approach described in Section 2.3.By considering the contribution of Equation ( 11) to the current density in Equation ( 1), its impact on u Δ can be evaluated, which in turn again affects the gain g.

Numerical Model
While the analytical model in the section above provides valuable insights, it does not self-consistently include the back-action of the optical dynamics on u Δ , which is at best justified for active modulation close to threshold.This limitation can be overcome by a comprehensive numerical model, fully accounting for the microwave-optical co-propagation dynamics and resulting effects, such as the self-organized formation of oscillation patterns.
We have developed a stable and efficient numerical framework for solving the Maxwell-Lindblad-transmission line model, consisting of Equation (3), Equations (A1)-(A5), and Equation (1) along with the adequate boundary conditions.The optical propagation equation, Equation (3), is discretized using the Risken-Nummedal (RN) finite differences scheme, [84] taking advantage of its favorable numerical properties. [46]For the DM part, Equations (A1)-(A5), we use an explicit Adams-Bashforth three-step method, [46] featuring good efficiency and stability.The transmission line equations, Equation (1), are solved using a special variant of the finite difference time domain (FDTD) method which incorporates frequency dependent elements. [57]For this purpose, L ′ and R ′ are represented in terms of a Debye rational approximation, the parameters of which are here extracted using Bayesian optimization.We have extended this approach to incorporate the current density J(x, t) contained in Equation (1).J, u, E ± , and the DM elements are defined on the same spatial grid.Using a grid spacing Δ x , the time step of RN is given by Δ t = Δ x ∕v g , [46,84] while FDTD requires a smaller interval due to the Courant stability criterion.This is here achieved by performing two FDTD update steps during Δ t .FDTD requires a staggered grid, that is, i is defined on a spatially and temporally shifted grid relative to u, and the J grid is temporally shifted.Based on this grid, the coupling of the equations via E ± , p ± and u (see Figure 2) can be straightforwardly implemented, while the values of J at the required time steps are determined by extrapolating the data obtained from Equation (B1).In the following, we present exemplary simulation results for an actively mode-locked short pulse [21] and a freerunning frequency comb [13] structure.

Example Actively Mode-Locked QCL Structure
While active mode-locking of QCLs can be described within the conventional MB/ML equations, [81] modeling of microwaveoptical co-propagation effects requires the MLT framework.We exemplarily simulate a microstrip QCL featuring a 3.1 mm long and 10 m thick metal-metal waveguide, which is used for the generation of ultrashort THz pulses in the ps range. [21]Our active region model, illustrated in Figure 3, is adapted from a related structure [85] which this three-well design is based on.In our model, scattering between all states is taken into account.The periodic arrangement of quantum well stacks is implemented in the usual way by taking a representative period and adding suitable boundary conditions (see Appendix B).The states of the representative period are in Figure 3 arranged such that coherent transport due to tunneling and optical transitions only takes place inside the period; that is, interperiod transport is exclusively mediated by scattering.The active region parameters, such as eigenenergies, dipole moments, coupling energies, and scattering as well as dephasing rates, are extracted from self-consistent simulations based on the density matrix Monte Carlo approach coupled to the Schrödinger-Poisson equation system. [71]In the tight-binding picture used for our model in Figure 3, all the relevant microwave transitions between closely aligned states coincide with tunneling transitions. [65]The bias dependence enters our model via E ′ n (see Equation ( 6)), which causes an energetic detuning of the tunneling transitions.The resulting stationary current-voltage characteristics features a differential resistance close to the measured value, [21] confirming the validity of our model.The transmission line parameters L ′ , R ′ , and C ′ have been determined in Section 3.
In the following, we compare our simulations to available measurements [21] for actively mode-locked operation (cf. Figure 4a in ref. [21]).Like in experiment, [21] operation close to threshold is assumed for our simulation.We apply direct bias modulation at the left microstrip end (see also Section 3) in order to have a controlled modulation amplitude, which is not affected by frequency dependent losses in the feed microwave circuitry and imperfect impedance matching.We assume a bias modulation of 0.8 kV cm −1 corresponding to ≈ 26 mW, which is realistic considering the experimentally injected microwave Bias [arb.u.] power of 450 mW [21] and the power attenuation of ≈ 12 dB in the circuitry, measured for a comparable setup at 13 GHz. [79]We use the MLT equations to investigate the degree of mode-locking in dependence of the modulation frequency by evaluating the power and phase noise quantifiers. [86]Similar to experimental observations, [21] the most stable mode-locked operation with clearly separated periodic pulses is not obtained for modulation at exactly the free-running repetition rate (13.05GHz in the simulation), but rather for a slightly lower modulation frequency of 12.86 GHz.In Figure 4, the resulting microwave and optical propagation dynamics in the waveguide cavity is illustrated for a single cavity roundtrip (for the animated movie, see Movie S1 in the Supporting Information).In this case, the pulse precedes the voltage bump, and thus the gain peak, for the most part of the roundtrip, as also expected from standard active modelocking theory for the averaged pulse dynamics. [87]The pulse outcoupled at the right facet (light blue central vertical line in Figure 4) shows an asymmetry similar to experiment [21] with a steep leading edge, and its full width at half-maximum (FWHM) duration of 16 ps is comparable to the experimental value of 11 ps (cf. Figure 4a in ref. [21]).We note that at somewhat lower modulation frequencies, the simulation yields even shorter pulses (e.g., 11 ps at 12.7 GHz), but then the noise quantifiers do not quite fulfill the mode-locking condition. [86]The optical pulse generates a photon-induced current component in the gain medium, resulting in a local discharging of the transmission line capacitance and thus producing a voltage dip on the microwave signal.The modulation bias at the open end of the transmission line oscillates between −0.192 and 0.163 kV cm −1 .The minimum value is in accordance with the analytically predicted reduction of the modulation amplitude to 26% (see Section 3), while the maximum value is somewhat smaller due to the voltage dip.As noted above, this back-action of the optical dynamics on the microwave signal is not self-consistently included in the analytical theory presented in Section 3.This effect becomes especially important for high optical intensities and in the context of self-seeded microwave modulation, as discussed in the next example.
To demonstrate the beneficial effect of the electronic modulation co-propagating with the optical pulse, we have repeated the simulations with roughly doubled microwave damping a  , achieved by quadrupling the resistance to R ′ = 0.18(f ∕Hz) 1∕2 Ω m −1 .While at the left microstrip end, the modulation remains unchanged, it experiences stronger damping along the transmission line, and the roundtrip-averaged modulation amplitude defined in Equation ( 9) reduces to 63% of its previous value.Consequently, the simulated pulse duration increases by ≈ 30%, and the noise quantifiers indicate only partial mode-locking.This clearly confirms the beneficial effect of exploiting microwave-optical co-propagation to enhance the modulation depth, as also suggested by experimental evidence. [17,21]

Example Free-Running Frequency Comb QCL Structure
As a second example, the MLT equations are applied to a freerunning THz frequency comb design. [13]This structure has previously been simulated using ML-type equations, [48,49,59] that is, without considering microwave propagation.It has been experimentally and theoretically confirmed that synchronized microwave modulation supports the comb formation process. [23,83][83] Alternatively, the experimentally observed emergence of radio frequency (RF) beatnotes on the metal-metal waveguide at the roundtrip frequency and multiples thereof [13,23,27] can potentially be exploited in a free-running device to generate a self-seeded microwave signal, resulting in the self-organized formation of microwave-optical oscillation patterns.This approach is highly attractive since it does not require the dedicated external circuitry and electronics needed for active modulation.In this context, the use of an adapted waveguide for a targeted enhancement of the microwave signal has been proposed. [23]This self-seeding effect is not contained in the MB/ML formalism; rather, the MLT equations are required for proper modeling.Especially, this model can clarify to what extent the waveguide design matters, that is, if transmission line effects beyond purely local oscillations play a role.
The four-well active region of the frequency comb structure is modeled with a tight-binding scheme and periodic boundary conditions similarly as above, employing the parameters used in previous simulations. [48,49]We set  2 = 0 in Equation (3), corresponding to the not correctly compensated case [48] since residual GVD caused by the quantum system is still present.In contrast to active mode-locking investigated in Section 4.1, we do here not impose a bias modulation at the left microstrip end, but rather assume a constant current injection while self-sustained bias oscil- lations may occur.Since here, as for active mode-locking, metalmetal waveguides appear to be especially suited, [23] we tentatively assume the same active region dimensions as above, and use the same L ′ and R ′ as in Section 4.1.C ′ is treated as a degree of freedom to explore the influence of the microwave effective index on device operation, and in particular to investigate if a favorable parameter range corresponding to the above-mentioned adapted waveguide exists.We have performed MLT simulations over 15 000 roundtrips for a large range of C ′ between 0.3 and 2.1 nF m −1 .The last 5000 roundtrips were evaluated, resulting in a two-lobed discrete spectrum for all values of C ′ , as also obtained in experiment and previous MB simulations not including microwave propagation effects. [13,48,49]However, the degree of self-locking, characterized again in terms of the power and phase noise quantifiers, greatly depends on C ′ .Between ≈ 1.2 and 1.3 nF m −1 , stable self-locking is obtained, as also indicated by the temporal periodicity of both the optical and microwave field, as well as the narrow spectral widths of the individual comblines and the RF beatnote.For the other cases, irregular or multiperiodic behavior arises.Also the spectral comb shape depends on C ′ , albeit to a lesser extent, while the average outcoupled intensity varies only slightly between 6.55 and 6.78 kW cm −2 .
In the following, we exemplarily discuss simulation results for C ′ = 1.22 nF m −1 where stable self-locking is obtained, and for C ′ = 0.78 nF m −1 .The first case represents a design with an adapted waveguide, while the second case produces an unstable comb, with similar shape and power/phase noise quantifier values as obtained from conventional MB/ML simulations of this not fully dispersion compensated setup.In  shown for both cases, as well as the associated optical and microwave spectra.For the adapted waveguide (Figure 5a,b), the intensity and bias time traces are highly periodic, while for C ′ = 0.78 nF m −1 (Figure 5c,d), irregularities appear.The associated optical and microwave power spectral densities are shown in Figures 5e and 5f,g, respectively.The total optical comb power is almost identical for both cases as mentioned above, but the overall comb shapes differ somewhat, and more importantly, the individual spectral lines are considerably broadened and thus have a much lower peak value for the not adapted waveguide (see Figure 5e).This broadening is also reflected by the fundamental RF beatnote at f R ≈ 13.1 GHz, shown in Figure 5f: For the not adapted microstrip, the FWHM linewidth is 10.2 MHz, while in the adapted case, it drops well below our numerical resolution limit of ≈ 2.7 MHz, and a similar behavior is obtained for the intensity beatnote (not shown).The full RF spectrum consists of regularly spaced lines and extends to several 100 GHz, as shown in Figure 5g and experimentally observed for a similar device. [27]Also here, for the not adapted waveguide, the lines are much broader, and thus have a significantly lower peak value, while the overall spectral power is comparable in both cases.In Figure 6, the microwave and optical propagation dynamics over a cavity roundtrip are illustrated for the structure with the adapted waveguide (for the animated movie, see Movie S2 in the Supporting Information).The microwave bias oscillation has a similar fundamental frequency as in Figure 4 due to the identical waveguide parameters (apart from the slightly different C ′ ).Also the arising modulation pattern looks somewhat similar but is less smooth, giving rise to the higher harmonics visible in Figure 5g.The intensity, on the other hand, shows a much weaker modulation than in Figure 4, as is typical for QCL frequency combs. [12]e spatiotemporal variation of the current density for the setup with the adapted waveguide, here characterized by its standard deviation over a roundtrip, is ΔJ I = 0.0194 kA cm −2 .This value can directly be compared to the analytical theory in Section 3: Using the standard deviation of the gain-intensity product Δ[(g − a)I] = 59.1 kW cm −3 extracted from the simulation, Equation ( 11) yields with L p =54.7 nm and  c ∕(2) = 3.7 THz an estimate of ΔJ I = 0.0211 kA cm −2 , in good agreement with the numerical result.
Besides assuming a constant current at the left microstrip end, simulations of this structure have also been performed imposing a constant voltage with or without an additional 3.1 mm long feed line, yielding stable self-locking around a similar value of C ′ as above.Interestingly, with the feed line self-locking is obtained over a broader range of C ′ (≈ 1.0..1.3nF m −1 ), indicating that the electric boundary conditions imposed by the external feed circuitry should also be considered in the design process of freerunning hybrid photonic and microwave-electronic QCL devices.
As already mentioned, the analytical transmission line model of Section 3 does not self-consistently include the impact of the optical dynamics on the microwave signal and, thus, is especially suited for external modulation and operation close to threshold.By contrast, the active region of the free-running structure investigated here cannot be represented by a simple conductance G ′ as in previous work [17] and Section 3, since it features a strong photon-induced current component driving the microwave dynamics.This clearly demonstrates the need for a comprehensive model including the full optical, microwave, and active region dynamics for an in-depth analysis and targeted design of hybrid quantum devices, as presented in this paper.We have also performed simulations with increased total microwave damping a  L across the waveguide, achieved either by doubling a  as for the actively mode-locked structure or by extending the length to L = 5 mm.No self-locking is obtained since the noise quantifiers significantly exceed the mode-locking threshold, [86] although for L = 5 mm, they do get noticeably reduced around a certain C ′ .
Similarly as for the actively mode-locked structure, the obtained results suggest that at moderate losses, microwave propagation effects beyond purely local oscillations can build up along the transmission line, enabling the observed self-locked operation for adapted values of C ′ .

Conclusion
In conclusion, a theoretical model for hybrid microwaveelectronic and photonic quantum devices has been introduced in the form of the MLT equations, along with a suitable numerical framework.Simulations of both directly modulated and freerunning hybrid QCL structures yield excellent agreement with available experimental results and demonstrate the importance of microwave-optical co-propagation effects.As analytically and numerically shown in accord with experiment, the hybrid concept enables more efficient device modulation, and thus improved active mode-locking.Moreover, our simulations of free-running devices reveal the emergence of broadband hybrid terahertz-microwave generation, as recently observed in experiment.Notably, for structures with suitably engineered waveguides, we numerically demonstrate the self-organized formation of regular microwave-optical oscillation patterns, potentially Laser Photonics Rev. 2023, 17, 2300461 enabling self-stabilized comb operation.[37][38][39] Furthermore, in analogy to simplified models for phase-locked operation of quantum optoelectronic devices obtained from the MB equations, [42][43][44][45] the MLT framework may serve as the starting point for the development of (semi-)analytical models for the hybrid operating regime.

Appendix A: Density Matrix Equations
We consider a quantum system with one or (in case of multiple laser levels) several relevant optical transitions (d mn ≠ 0) in near-resonance (| mn | ≈  c ).Furthermore, the system can contain tunnel coupling (Ω mn ≠ 0) and relevant microwave transitions (d mn ≠ 0) between closely aligned states (| mn | ≪  c ).We assume a Fabry-Pérot-type cavity where counter-propagating waves with amplitudes E + and E − form a standing wave pattern, which results in SHB.Diffusion counteracts the formation of an inversion grating and should thus be included in the model. [46]hile the diffusion coefficient is often assumed to be identical for all levels, this is for example not true for quantum dot media where diffusion mainly takes place in the reservoir. [67]Furthermore, there is a debate if diffusion should also be considered for the coherences. [88]Thus, we use a generalized approach, adding a term  t  mn = ⋯ + D mn  2 x  mn with diffusion coefficient D mn to the evolution equation of each DM element.We start with Equation ( 4), taking the Hamiltonian ĤS = Ĥ + ĤRWA .The matrix elements of Ĥ given by Equation (6) imply that also the transition frequencies  mn = Ω mm − Ω nn are bias dependent.For employing the RWA, we substitute Equation (2) into Equation (4).Furthermore, for all near-resonant optical transitions with d mn ≠ 0, we substitute  mn = ∑ ±  ± mn exp[±ix − i c sgn( mn )t] for the corresponding off-diagonal DM element.As for SHB, the inversion grating's periodicity corresponds to that of the optical power.Thus, we make for all the remaining DM elements the ansatz  mn =  0 mn + ∑ ±  2± mn exp(±2ix).The equations in RWA are then obtained by discarding terms oscillating rapidly in time.For a more compact notation, we introduce r m = ∑ n≠m r mn ,  mn =  2 D mn and f ± with f + = E + , f − = E * − .For the diagonal DM elements, we obtain and  2− nn = ( 2+ nn ) * .For the off-diagonal DM elements  ± mn of optical transitions (d mn ≠ 0) with  mn > 0, we get  t  ± mn = i( c −  mn ) ± mn + i and  2− nm = ( 2+ mn ) * .The remaining off-diagonal DM elements  mn are set to 0 in our RWA model since carrier transport between the corresponding levels m and n is then exclusively mediated by incoherent transport, considered in Equation (A1) and (A2) by the rates r mn .The polarization term in Equation ( 3) is obtained from the  ± mn in Equation (A3) as d nm  ± mn , (A6) with the electron number density n 3D and overlap factor Γ. This model can straightforwardly be adapted to QD lasers by adding suitable nonlinear rate terms describing carrier exchange among the QD states and with the wetting layer, and by considering both the electron and hole dynamics which is most easily done using the excitonic approximation. [68]For compatibility with the MLT model, carrier injection must be included in dependence of the bias voltage rather than the current density. [61]

Appendix B: Treatment of Periodic Structures
QCLs consist of many equivalent periods of length L p , corresponding to coupled, ideally identical quantum systems.Usually, periodicity-breaking effects, such as the influence of the electric contacts, the transverse dependence of the optical field and possibly domain formation, are neglected, and the QCL active region is modeled by a representative period, described by a quantum system containing the N relevant states of a single period with suitable periodic boundary conditions. [89]As discussed in Section 2.3, tunneling across thick barriers (such as the injection and/or extraction barriers in many QCL designs) can be included by using tightbinding states rather than an extended basis. [62,90]The tight-binding basis set is obtained by dividing the quantum structure into modules separated by the thick barriers, and solving the Schrödinger-Poisson equation for each module separately. [62,90]An alternative method for obtaining a suitable localized basis set without having to identify thick barriers are the EZ states, obtained from a basis transformation of the extended energy eigenstates. [65]Apart from the optical and tunneling transitions, the transport in QCLs is largely incoherent. [71,90,91]Thus, the states of a representative period can usually be chosen in the simulation such that tunneling as well as optical transitions happen only inside the period and interperiod transport is incoherent, that is,  mn = 0 if the levels m and n are in two different periods.This enables an implementation of periodic boundary conditions in the same way as for rate equation models of QCLs, [92,93] using that  mn =  m±N,n±N .Taking into account only scattering to adjacent periods, periodic boundary conditions can then be implemented into our model by replacing r mn in Equations (A1) and (A2) with the effective rates r ′ mn = r mn + r m,n+N + r m,n−N , and making sure that the total outscattering rate r m from level m also contains scattering to adjacent periods.In this case, also the evaluation of the current density is possible in a simplified manner, yielding

Figure 1 .
Figure 1.Schematic representation of the microstrip QCL as an example of a hybrid optoelectronic and microwave quantum device.The microstrip is represented by its equivalent circuit for an infinitesimal segment, coupling to the quantum active region.The parameters and variables included in the figure correspond to those introduced in Section 2. In the lower right corner, the coordinate system used in our model is shown.

Figure 2 .
Figure 2. Schematic illustration of the coupling between the DM, optical propagation, and transmission line equations.The current density J and polarization p ± are computed from the DM elements  ij .More specifically, in our model, p ± is obtained from the relevant off-diagonal elements  ij in the rotating reference frame.

Figure 3 .
Figure 3. Conduction band diagram of the THz active region in the tightbinding picture.The probability densities for the relevant states belonging to the representative period (full lines) and to adjacent periods (dashed lines) are shown.The arrows indicate the dominant transport mechanisms, namely scattering (block arrows) and coherent transport due to optical (wavy arrow) or tunneling (circular arrows) transitions.

Figure 4 .
Figure 4. Simulated optical power and modulation bias voltage along the waveguide for a single intracavity roundtrip in the unfolded cavity representation, where the right half of the figure shows the backward propagating optical pulse.The time t is expressed in units of roundtrip time T R .The central vertical line corresponds to the right facet, and the discontinuity of the optical power is due to facet outcoupling.

Figure 5 .
Figure 5. Simulation results for free-running QCL comb device.a-d) Outcoupled intensity and AC bias voltage at the right facet for an adapted (a,b) and a not adapted (c,d) microstrip as a function of time (in units of roundtrip time).e-g) Power spectral density of the optical field (e) and the bias voltage (f,g) for the cases shown in (b) and (d).

Figure 5 ,
the timedependent intensity and modulation bias at the right facet are Laser Photonics Rev. 2023, 17, 2300461

Figure 6 .
Figure 6.Simulated optical power and modulation bias voltage along the waveguide for C ′ = 1.22 nF m −1 , represented in analogy to Figure 4.
0 mi + E ∓  2± mi )− ( mn +  mn ) ± mn , (A3)and use  ± nm = ( ∓ mn ) * to obtain the remaining elements.The condition | in | ≪  c in above summation implies that we only include terms with slowly oscillating elements  0 in and  2± in since the other terms are discarded in the RWA, and analogously for | mi | ≪  c .For these DM elements, we obtain the evolution equations t  0 mn = i ∑ i ( Ω in  0 mi − Ω mi  0 in ) −  mn  0 mn + i 2ℏ