Chaotic Behavior of Quantum Cascade Lasers at Ignition

The ignition of Quantum Cascade Lasers can occur from a state of oscillating fie ld domains. Here, the interplay between lasing and the kinetics of traveling dom ain boundaries provides complex oscillation scenarios. We analyze our numerical findings in detail for a device operating at terahertz frequencies and manifest chaotic evolution by positive Lyapunov exponents. This shows that these importan t devices can exhibit chaotic behavior even without periodic driving, which need s to be taken into account in their design.


Introduction
Negative differential conductivity (NDC), i.e., the decrease of current with increasing electric field, is a common source of instabilities in semiconductor devices [1,2,3]. In extended systems like the Gunn diode [4], it leads to the formation of spatial domains with different electric fields.
The Quantum Cascade Laser (QCL) [11,12] is currently the most important device for mid and far-infrared radiation. QCLs are based on carefully designed semiconductor heterostructures, see, e.g., Fig. 1(b). These guide the electron flow by tunneling and scattering to establish electronic inversion for a pair of quantum levels (the laser levels) at a specific electric field, the nominal operating point (NOP). In order to increase the total gain, a module of several layers including the laser levels is repeated several times, so that the electrons traverse the total structure like water in a cascade. Thus, the field distribution in QCLs exhibits domain formation if driven in an NDC region [13,14,15]. As the devices are most efficient, if all modules contribute equally to the gain, it is a common strategy to avoid NDC around the NOP in the QCL-design. On the other hand, resonant tunneling is prone to provide NDC above alignment [16,17,18]. Therefore, NDC is ubiquitous in layered structures such as QCLs and instabilities close to threshold [19,20] are not always avoidable.
Here we focus on device V812 from [19], a QCL operating with good performance at terahertz frequencies, where the NOP is actually in the NDC region. Recently, some of us showed [21], that ignition occurs in the state of oscillating field domains and that the arising lasing field afterwards stabilizes the behavior, see Fig. 1(a). In the transition region, where lasing starts and coexists with domain formation, see Fig. 1(c), our numerical simulations provide interesting complex dynamics including chaos, which we analyse in detail here. While chaotic behavior had been recently found in QCLs under external periodic driving [22,23], we note, that our system is autonomous.
Our article is organised as follows: In Sec. 2 we briefly repeat our model detailed in [21].
Here we focus on the differential equations describing the time evolution of electric fields F k in the modules and the occupations of the relevant lasing modes N i ph , which are our main variables. Section 3 presents a detailed analysis of the spatio-temporal evolution in the QCL just after ignition.
Subsequently, we show in Sec. 4, that the irregular behavior observed exhibits positive Lyapunov exponents, which proves that we observe chaos.

Model
In order to study the formation of field domain formation, we consider the dynamical evolution of the (average) electric field F m in module m of the QCL containing N = 222 modules of thickness d = 45.07 nm, see Fig. 2(b). In full analogy to superlattices [7,8] and earlier QCL studies [24], we where r = 12.9 is the average relative permittivity, C s is the capacitance of the QCL structure and C p is a parasitic capacitance in parallel to the device, see Ref. [25] for a derivation. In the following we assume C p = C s /4. The current J(t) is fed via the circuit shown in Fig. 2(a) resulting in is the total bias drop over the QCL (including a boundary region with field F 0 ). Here, U 0 is the external bias applied to the device, which is the control parameter for our system. From the experimental setup [21] we extract a load resistance R L = 41.2 Ω and a probe resistance R p = 1050.4Ω. Finally, A = 0.15 mm 2 is the cross section of the QCL and V B = 0.8 V reflects the bias drop due to a Schottky barrier at the metal-semiconductor contact.
The currents J k→k+1 between the modules are determined as follows: For k = 1 . . . N − 1, we use an expression J k→k+1 (F k , n k , n k+1 , {N i ph }) based on our non-equilibrium Green's function (NEGF) scheme [26] (using 7 states per module). The homogeneous results in Fig. 1 for the areal doping density of n D = 3 × 10 10 /cm 2 (located in the center of the largest wells). For the inhomogeneous case, the areal electron density in module k is given by Details are given in Ref. [21]. The boundary currents at the beginning and the end of QCL structure are estimated by a phenomenological conductivity σ = 0.15 A/Vcm using J 0→1 = σF 0 and J N →N +1 = n N σF N /n D in analogy to Refs. [7,24]. We use a lattice temperature of 77 K throughout this work.
The occupations N i ph of the cavity electromagnetic modes i with frequency ω i 0 are changing due to the interplay by gain G from the QCL medium and waveguide/mirror losses (quantified by the threshold gain g th = 20/cm) as where n g = 3.6 is the group refractive index assumed to be constant here. We also considered spontaneous emission with a time τ i sp = 3 ms, where n U LS k is the areal carrier density in the upper laser level in module k. The gain G(ω) is the sum of the gain contributions G(F k , ω, {N i ph }) for all modules k, which are extracted from our NEGF calculations, see Ref. [21] for details.
Our model provides a closed system of equations for the fields F k with k = 0, 1, . . . N and the photon occupations N i ph , where we have 38 relevant modes with frequencies between 10 and 16 meV in the cavity. We tacitly assumed that the internal electron dynamics inside the QCL is instantaneously adapting to the actual fields and mode occupations. This is probably a good approximation, as typical scattering times are shorter than 1 ps. In comparison, the photon lifetime is n g /(g th c) = 6 ps and the dielectric relaxation time is at least r 0 (dF/dJ) = 5 ps, based on the maximal slopes in Fig. 1(a).  detailed in Ref. [21]. This shows the validity of our simulations, which have no fit parameters except for the contact conductivity and assuming a higher lattice temperature than in the experiment, which mimics heating of the phonon distribution [27,28].
The central input to the model are the functions for gain and current, which show a wide variation with the system parameters as displayed in Fig. 3. The detailed fitting process for the gain and the current is discussed in [21]. The mode frequencies chosen for these plots span the whole range and demonstrates that the data are highly frequency-dependent. For the linear response, the current density is largely unchanged. For higher intensity, gain saturation is observed together with an increase in the current due to the stimulated intersubband transitions.

Oscillating Field Domains and Chaos
In this section, we analyse the dynamical behaviour in the region shown in Fig. 1(c) in detail.
When U 0 becomes larger than 54.2 V, the operation point reaches the NDC region for a homogeneous field distribution. This causes the formation of field domains with boundaries travelling through the device. This provides oscillations in current and bias and therefore the corresponding timeaverages have been plotted in Fig. 1(c). Without lasing, Fig. 1(a) shows that the condition of  Fig. 1(c). This lasing field strongly modifies the current and gain as demonstrated in Fig. 3 which results in the complex behaviour we observed. Here, we identified 5 distinct regions with essentially different behaviour as shown in Fig. 4.
In region I, Fig. 4(a) shows the characteristic oscillations due to travelling field domain boundaries. The homogeneous field distribution becomes unstable (e.g. at t = 12 ns) and splits up in a high-and low field domain. Afterwards the electron accumulation layer separating both domains travels to the positive contact and the scenario repeats, after the field in the low-field region has increased to maintain the bias. Just after its formation, the high-field domain can be large enough to provide sufficient gain to compensate losses, but this holds only for a short time, so that the lasing intensities are never sufficient to effect the behaviour in this region.
With increasing bias U 0 the high-field domain becomes more extended and the lasing field becomes stronger, so that it significantly enhances the current around 63 mV, see Fig. 3(d). In order to keep the current density, the field in the high-field domain needs to diminish as can be seen Fig. 4(b), which is characteristic for region II, where the average current increases stronger with bias than in region I. However, with shrinking length of the high-field domain, the gain drops and the original bias per module of 63 mV is restored in the high-field domain before a new instability appears in the low-field domain associated with a peak in current.
In region III, as the external drive increases, the lasing intensity becomes more strong and covers a wider range of each domain cycle, as seen in Fig. 4(c). However, lasing always stops before a new domain boundary forms. In the local U qcl − J relation (see Fig. 1(c)), U qcl stabilizes while the intensity continues increasing.
In region IV, lasing starts to persist most of the time and some high field domains form around the NOP and the average bias drops as seen in Fig. 4(d). Finally, in region V, lasing persists all the time and all the high field domains form around the NOP. U qcl almost stabilizes as the current increases slightly with the intensity. (see Fig. 4(e) and Fig. 1(c)).
As discussed in Fig. 4, the system shows fundamentally different behavior in the five regions.
In order to check for periodicity, we plot the local maxima of the U qcl (t) signal as a function of the  Oscillations observed in GHz range.
control parameter in Fig. 5. Such a diagram is well-known for identifying routes to chaos [29,22].
In order to identify irregular behavior, the time series U qcl (t) is divided into three equal interval in times marked with colors pink, orange and black after removing initial behavior in the first 160 ns. The black dots are from the latest time interval and cover all earlier dots with the same value.
The orange dots are from the middle interval and their persistence indicates non-periodic solutions.
Finally, pink dots from the first period in time indicate that periodicity had not been reached within the first 160 ns, and their presence without yellow dots indicates long transients rather than chaos [30].
In region I, consistently with Fig. 4(a), oscillations are regular. A couple of points in region II exhibit regular period two and three behaviour. In region III, there is a point U 0 at 56.1V that shows some irregularity which we did not analyse further. In general U qcl solutions continue to be regular. Shortly after the transition between region III to IV, we have strong indications of chaos among the region IV. The region V again has a regular structure. In the following we focus on the transition range which is highlighted as cyan color background. In Fig. 6, we show the U qcl time series for three U 0 points from the cyan highlighted region in In contrast irregular behaviour is observed for slightly larger biases U 0 = 56.95 V and U 0 = 57 V as shown in panels (b) and (c). These irregular patterns are the first indications of chaos [31].
In Fig. 6, Fast Fourier Transforms (FFT) of the corresponding three U 0 points are shown in the insets locate on the right of the panels with oscillations in GHz range. Here, in the transition between regular to irregular oscillations, the structure clearly becomes continuous rather than discrete. This is not a proof of chaos but an indication consistent with the time series.
With constructing the time series and FFT of the three U qcl data discussed in Fig. 6, we will show these rigorous numerical results providing chaos. The most common methods of proving chaos are usually constructing the phase spaces and deriving the Lyapunov exponents. To understand the dynamical evolution of the system, phase diagram is one of the main steps to investigate the chaos [32]. To analyze the results we have chosen U qcl and the time derivative of U qcl as two variables out of our complex system with 223 variables to define a new reduced phase diagram. For simplicity, this is referred to as phase diagram in the following. In Fig. 7, phase diagrams of the three U 0 data from Fig. 6 are shown. In Fig. 7(a) and (d), one can clearly see that the regular data is following the same path in each period without any deviation. Therefore, small perturbations of the data do not change any dynamics of the system. However, In Fig. 7(b), (c), (e) and (f) one can clearly observe the deviations of the path spread around the phase diagram. While the trajectories come close to each other at some places, tiny differences grow to qualitative different behaviour in the course of the time-evolution, as characteristic for chaos. There are some similar phase diagram constructions done in other works including tangential junctions with stable and unstable manifolds (see Ref. [33,34,35,36,37]).
In Fig. 8 we analyse the behavior for the periodic oscillation at U 0 = 56.90V in detail. The phase diagram in panel (a) shows, that two subsequent periods (marked by red plus and green cross symbols) lie on top of each other. The lasing intensity [red dashed line in panel (b)] is essentially dropping to zero around symbol 5 before the new domain boundary forms, which is associated with a sharp peak in the U qcl (t) signal at symbol 1. A smaller bias peak arises close to symbol 2 just after the domain formation. This behavior is typical for region III as discussed in Fig. 1(c) and to appear as seen in Fig. 9. Here the same method is used to plot the variables as in Fig. 8 but with three major peaks of bias marked with 1. Three markers (red plus, green cross and magenta circle) now represent the trajectories after the respective bias peak in the phase diagram . We see, that the trajectories deviate significantly in each cycle. These deviations bring some new scenarios in the trajectories. While the red and magenta markers follow a similar path with a significant deviation, the green trajectory stays on a whole different path. Eventually all the markers get together around 9.75V in Fig. 9(a). It appears this point of junction is the place where the system decides how to evolve. Also, as seen in Fig. 9(b) and (c) lasing stays persistent even while new domain boundaries form in the middle part with green cross markers. In contrast, at the major bias peaks, the domain boundaries form a state with vanishing intensity.

Lyapunov Exponents
As we seek to quantify the presence of chaos in the QCL we consider the largest exponent from the Lyapunov spectrum of the system [38,39]. Let d(t) be the distance between two closely lying state vectors in a bounded phase space at time t then the largest Lyapunov exponent λ describes the time-evolution of the distance When λ is positive, small deviations in the state vector lead to exponentially growing deviations in phase space, i.e. chaos (provided the phase space is bounded).
To estimate λ from the time series we employ a procedure similar to the ones in Refs. [40,41] albeit with the difference that we have access to the true phase space of the system, i.e. the state of the modules, hence, we can ignore the question of how to reconstruct the attractor. Consider the field of the m'th module with time series F m (t) and define its normalized variable with F m and std(F m ) being, respectively, the mean and standard deviation of F m (t) over time.
This yields the normalized state vector at time t: v(t) = (v 1 (t), v 2 (t), · · · , v N (t)). From the normalization of the components of v(t) the variation of every variable is included equally in the analysis.
For two state vectors at times t i and t j , respectively, we define the Euclidean distance after the evolution time t evo as When estimating the Lyapunov exponents we traverse a fiduciary trajectory of v(t). At time t i we search the time series for another time t j for which d ij (0) is minimal. Here, we exclude any state vectors that are within the interval (t i − ∆; t i + ∆) to avoid trivial pairing between the state vector and its negligibly shifted self.
In a system of limited phase space volume the distance between any two state vectors is bounded.
Hence, the deviation between the pairs will saturate and the picture of exponentially increasing distances between any two close pairs only holds for small initial deviations. This is ensured by restricting the sampling to pairs with initial distances below a certain discriminator d ij (0) < D.
From the identified pairs t i and t j we estimate the contribution to the Lyapunov exponent at t evo as The average contribution of all pairs (i, j) yields the estimated Lyapunov exponentλ(t evo ) = λ ij (t evo ) .
The procedure above neglects the phase orientation information between one pair and the next when determining the sequence of pairs which is otherwise included in the algorithm of Wolf et al. [42]. Yet, as we are only interested in the largest Lyapunov exponent such phase information is unnecessary [41]. Considering first Fig. 10(a) with bias 56.95 V, the curves show that the estimated Lyapunov exponents (solid blue and dashed orange lines) exhibit a plateau for evolution times around 2 ns, the plateau resembles the stationary of the Lyapunov exponents for certain evolution times observed in in Refs. [43,42,41,23]. The plateau yields a positive estimate for the Lyapunov exponent of the order of 6 ns −1 > 0 which corroborates our interpretation of the irregular pattern in Fig. 9 as chaotic behaviour. Beyond the plateau the estimate decreases according to 1/t evo -behaviour (the black dotted line), as expected from (9) with saturation bound on the distance. Repeating the analysis for Fig. 10(b) with bias 56.90 V we find that the estimates for both discriminators are two orders of magnitude smaller than for Fig. 10(a) and hover around zero, in agreement with the regular behaviour observed in Fig. 8.

Conclusion and Outlook
In this work, we analyzed chaotic behavior in a QCL without external driving. This occurred for a particular QCL operating in its NDC region, where traveling field domains form. Here the highfield domain exhibited gain, resulting in pulses of light during periods, where this high-field domain is sufficiently large. With increasing driving, the pulses get more pronounced and can modify the subsequent formation of a new domain boundary. Here, two different formation scenarios exist and their succession becomes irregular, resulting in the chaotic signal. For even higher driving, the lasing never stops and we recover an ordinary laser operation.
In phase diagrams, we could identify the points in which the trajectories deviate from each other. Here we observe the sensitive dependence on initial conditions as demonstrated by positive Lyapunov exponents. Thus, QCLs form a further autonomous system of technological relevance showing chaos.
It would be interesting to study this experimentally for undriven QCLs. Indeed the device studied here showed subharmonics of period three [21] at a single operation point, which can be seen as an indication for chaos [44]. More data are very welcome to reconstruct the phase space in similar devices which could unambiguously demonstrate chaos in these technologically important systems.